Optimal Penalization and Nonlinear Solver Convergence for a DG-Based Richards’ Equation Model of Variably Saturated Flows

Camille Poussel¹, Mehmet Ersoy¹, Yolhan Mannes¹, Aimed Ajroud², Frédéric Golay¹

¹Institut de Mathématiques de Toulon (IMATH), Université de Toulon, La Garde, France.
²École d’Ingénieurs SeaTech, Université de Toulon, La Garde, France.

DOI: 10.4236/jamp.2025.1311227                

Abstract

In this article, we study the convergence of an IIPG (Incomplete Interior Penalty Galerkin) Discontinuous Galerkin numerical method for the Richards equation. The Richards equation is a degenerate parabolic nonlinear equation for modeling flows in porous media with variable saturation. The numerical solution of this equation is known to be difficult to calculate numerically, due to the abrupt displacement of the wetting front, mainly as a result of highly nonlinear hydraulic properties. As time scales are slow, implicit numerical methods are required, and the convergence of nonlinear solvers is very sensitive. We propose an original method to ensure convergence of the numerical solution to the exact Richards solution, using a technique of auto-calibration of the penalty parameters derived from the Galerkin Discontinuous method. The method is constructed using nonlinear 1D and 2D general elliptic problems. We show that the numerical solution converges toward the unique solution of the continuous problem under certain conditions on the penalty parameters. Then, we numerically demonstrate the efficiency and robustness of the method through test cases with analytical solutions, laboratory test cases, and large-scale simulations.

Keywords

Porous Media, Richards Equation, Discontinuous Galerkin, Backward Differentiation Method, Incomplete Interior Penalty Galerkin (IIPG), Broken Sobolev Space, Picard’s Fixed Point, Minimal Regularity Solution

1. Introduction

The behavior of flows in variably saturated porous media can be modeled by the Richards’ Equation (RE). One of the key advantages of RE is its ability to represent the porous medium, incorporating both saturated and unsaturated zones. While it doesn’t consider the air phase, RE effectively incorporates the effects of gravity and capillarity, enabling the modeling of complex processes across various scales. Notably, RE is a nonlinear parabolic equation that can transform into an elliptic equation under complete saturation conditions.

The history of RE begins with Darcy’s law, which was formulated experimentally by Darcy in 1856 [1] for saturated porous media. This result was later extended to multiphase flows by Buckingham in 1907 [2], resulting in the Darcy-Buckingham law, which serves as the cornerstone for the derivation of RE. The equation was first established by Richardson in 1922 [3], although it was later attributed solely to Richards, who independently published the equation in 1931 [4]. Initial attempts to numerically solve the RE date back to the late 1960s with the works of Rubin [5] and Cooley [6]. From the 1980s, RE was extensively studied from both theoretical and numerical perspectives.

In this paper, RE is introduced by providing its expression and constitutive laws. As the main objective of this work is to solve RE using Discontinuous Galerkin (DG) methods, the weak problem associated with RE is given and its discretization using the Incomplete Interior Penalty Galerkin (IIPG) formulation. Additionally, an overview of the penalization method is provided. The fully discrete IIPG formulation is derived through time integration using the implicit Backward Differentiation Formula (BDF) method. Due to the nonlinear nature of RE, its fully discretized nonlinear formulation is linearized using Picard’s fixed-point method. Theoretical results related to the solution of stationary nonlinear elliptic problem are produced, including existence, uniqueness, and convergence results. Furthermore, an automatic calibration method is obtained for penalization parameters. The solution of RE using the previously mentioned IIPG formulation is implemented in an in-house numerical code named RIVAGE, which is then validated against numerical benchmarks.

2. Governing Equation

RE is a classical nonlinear parabolic equation used to describe flow in both unsaturated and saturated zones of an aquifer (for a detailed derivation of the equation, please refer to Clément’s 2021 thesis [7]).

The so-called mixed formulation of the RE, commonly used in hydrology, is

$\theta \left(h-z\right)-\nabla \cdot \left(\mathbb{K}\left(h-z\right)\nabla h\right)=0$ (1)

where $h:=\psi +z$ is the hydraulic head with $\psi$ the pressure head, $z$ is the elevation, $\theta$ is the water content and $\mathbb{K}$ is the hydraulic conductivity tensor.

The tensor of hydraulic conductivity $\mathbb{K}$ is split, in general, into two parts: the intrinsic or saturated hydraulic conductivity tensor ${\mathbb{K}}_s$ and the relative hydraulic conductivity $K_r$ :

$\mathbb{K}\left(\psi \right)={\mathbb{K}}_s{K}_r\left(\psi \right).$ (2)

The intrinsic hydraulic conductivity tensor ${\mathbb{K}}_s$ depend on the material of the porous media.

The relative hydraulic conductivity is a function of the pressure head controlling the behavior of groundwater flow within the porous media and it is defined as

$K_r\left(\psi \right)=\{\begin{array}{ll}1\hfill & \text{if}\psi \ge {\psi }_e,\hfill \\ K_{e,\text{law}}\left(\psi \right)\hfill & \text{otherwise}\hfill \end{array}$ (3)

where $K_{e,\text{law}}$ is given by empirical laws, see Table 1 and Figure 1. The quantity ${\psi }_e$ , corresponding to the entry of the air pressure, the pressure head transition value between the saturated and unsaturated zones. The saturated zone corresponds to $\psi \ge {\psi }_e$ and the unsaturated zone to $\psi <{\psi }_e$ . The water table corresponds to $\psi ={\psi }_e$ by definition.

Table 1. Hydraulic relations for hydraulic conductivity and effective saturation.

Name	Expression	Parameters
Gardner-Irmay relations (1954) [8]	$\begin{array}{l}{S}_e={\text{e}}^{\frac{\alpha \psi }{m}}\\ K_r={\text{e}}^{\alpha \psi }\end{array}$	$\alpha$ : Pore-size distribution $m$ : Tortuosity
Vachaud’s relations (1971) [9]	$\begin{array}{l}{S}_e=\frac{C}{C+{\left|\psi \right|}^D}\\ K_r=\frac{A}{A+{\left|\psi \right|}^B}\end{array}$	$A,B$ : Empirical shape parameters $C,D$ : Empirical shape parameters
Van Genuchen-Mualen relations (1980) [10]	$\begin{array}{l}{S}_e={\left(1+{\left(\alpha \left|\psi \right|\right)}^n\right)}^{-m}\\ K_r=S_e^l{\left(1-{\left(1-S_e^{\frac{1}{m}}\right)}^m\right)}^2\end{array}$	$l=0.5$ : Pore connectivity $\alpha$ : Linked to air entry pressure inverse $n>1$ : Pore-size distribution $m=1-\frac{1}{n}$ : Pore-size distribution

Figure 1. Hydraulic laws for effective saturation and hydraulic conductivity.

The water content law is expressed in terms of the effective saturation $S_e$ :

$S_e\left(\psi \right)=\frac{\theta \left(\psi \right)-{\theta }_r}{{\theta }_s-{\theta }_r},$ (4)

where ${\theta }_r$ is the residual water content and ${\theta }_s$ is the saturated water content corresponding to the minimal and maximal saturation, respectively. The effective saturation is defined as follows

$S_e\left(\psi \right)=\{\begin{array}{ll}1\hfill & \text{if}\psi \ge {\psi }_e,\hfill \\ S_{e,\text{law}}\left(\psi \right)\hfill & \text{otherwise},\hfill \end{array}$ (5)

where $S_{e,\text{law}}$ is given by empirical laws, see 1 and 1.

Remark. The nonlinear behavior of the constitutive laws $S_{e,\text{law}}$ and $K_{r,\text{law}}$ (see Table 1 and Figure 1) are responsible of the fails of the convergence of the numerical methods and a particular attention have been done. In particular, we have:

• In the saturated zone, hydraulic properties remain constant and RE becomes an elliptic equation characterized by fast diffusion.

• In the unsaturated zone, hydraulic properties approach very close to zero, which halts diffusion and can cause numerical inconvenience.

• For a specific set of parameters, when $\psi \to 0^{-}$ , constitutive laws may display extremely steep gradients.

To overcome, regularization techniques can be employed as in [11], for instance, which make slight modifications to the functions to avoid some types of degeneracy to improve convergence properties. In this paper, we will see that in the framework of DG, we show that whenever some numerical parameters are well-chosen, the modification of such constitutive laws is not necessary.

Equation (1) together with Equation (2) and Equation (4) can be completed with Dirichlet and/or Neumann boundary conditions as done in this work. One can also use more realistic boundary condition in view of real life simulation, such as the seepage boundary condition (we refer to [12] for details).

3. Numerical Methods

This section focusses on the presentation of the numerical solution of RE using DG methods. The solution is sought within a trial space due to the similarity of these methods to Finite Element (FE) methods, resulting in a weak problem.

Let $d\in \left\{1,2,3\right\}$ be the space dimension, the porous medium can be represented by the computational domain $\Omega \subset {ℝ}^d$ of boundary $\partial \Omega ={\Gamma }_D\cup {\Gamma }_N$ for which the subscript $D$ and $N$ stands for, respectively, Dirichlet and Neuman. Let $T\in {ℝ}_{+}^{\text{*}}$ be the final time.

The problem is:

Find $h\left(x,t\right):\Omega \times \left(0,T\right)\to ℝ$ such that:

$\{\begin{array}{ll}\theta \left(h-z\right)-\nabla \cdot \left(\mathbb{K}\left(h-z\right)\nabla h\right)=0,\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ h=h_0,\hfill & \text{in}\Omega \times \left\{0\right\},\hfill \\ h=h_D,\hfill & \text{on}{\Gamma }_D\times \left(0,T\right),\hfill \\ -\mathbb{K}\left(h-z\right)\nabla h\cdot n=q_N,\hfill & \text{on}{\Gamma }_N\times \left(0,T\right)\hfill \end{array}$ ()

where $h\in L^2\left(\Omega \times \left(0,T\right)\right)$ represents the solution of RE. Additionally, $h_0\in L^2\left(\Omega \right)$ , $h_D\in L^2\left({\Gamma }_D;\left(0,T\right)\right)$ , and $q_N\in L^2\left({\Gamma }_N;\left(0,T\right)\right)$ correspond to the initial condition, the Dirichlet boundary condition, and the Neumann boundary condition, respectively.

The matrix-valued function $\mathbb{K}$ depends monotonically on $h$ , is symmetric positive definite, and is uniformly bounded below and above (see Equation (2), Table 1 and Figure 1). Similarly, the function $\theta$ , also depends monotonically on $h$ , is uniformly bounded below and above (see Equation (4), Table 1 and Figure 1). Both $\mathbb{K}$ and $\theta$ are continuous functions within a given porous medium but may be discontinuous at the interface of heterogeneous materials.

3.1. Settings

The time duration $\left(0,T\right)$ is subdivided into $N$ time intervals such that $0=t^0<t^1<\cdots <t^N=T$ . Let $n\in \mathbb{N}$ , $0<n<N$ , if the time interval $T^n=\left[t^n,t^{n+1}\right]$ is considered, the corresponding time step is $\text{Δ}{t}^n=t^{n+1}-t^n$ .

Let us define ${ℰ}^n$ a partition of the computational domain Ω valid for all $t\in T^n$ . For the sake of simplicity, it is assumed that Ω is a polygonal domain in two space dimensions so that ${ℰ}^n$ covers Ω exactly. The mesh ${ℰ}^n$ is composed of quadrilateral and triangular elements not necessarily conformal.

For all elements $E\in {ℰ}^n$ , $d_E$ is its diameter defined as the ratio between its surface (${\mathfrak{s}}_E$ ) and perimeter (${\mathfrak{p}}_E$ ) and $d^n:={\max }_{E\in {ℰ}^n}\left(d_E\right)$ .

The set of all open faces of all elements $E\in {ℰ}^n$ is denoted by $ℱ$ . Moreover, one can define two subsets of $ℱ$ , ${ℱ}^{\partial }$ for the boundary faces and ${ℱ}^{\text{in}}$ for the interior faces:

$ℱ:={\displaystyle \underset{F\in \partial \Omega }{\cup }F}\text{and}{ℱ}^{\text{in}}:=ℱ\backslash {ℱ}^{\partial }.$ (6)

For a given element $E\in {ℰ}^n$ , there exists a set of face ${ℱ}^E:=\left\{F\in ℱ|F\in \partial E\right\}$ which defines boundaries of $E$ . Then, for all interior faces of $E$ , i.e., $\forall F\in {ℱ}^E\cap {ℱ}^{\text{in}}$ , there exists a neighboring element $E_r$ such that $E\cap E_r=F$ . Consequently, the normal unit vector $n_F:={\left(n_x,n_y\right)}^{\text{T}}$ pointing from $E$ to $E_r$ can be defined. An example of interior face is given Figure 2(a). Moreover for all boundary faces of $E$ , i.e., $\forall F\in {ℱ}^E\cap {ℱ}^{\partial }$ , there exists $E_{\partial }$ a fictitious element such that $E\cap E_{\partial }=F$ . Consequently, the normal unit vector $n_F$ pointing always from $E$ to $E_{\partial }$ can be defined.

Example 1. Figure 2(a) gives a graphical representation for an example mesh composed of triangles and quadrilaterals. In this example, the mesh is composed of 7 elements, i.e., ${ℰ}^n=\left\{E_i,i\in 1,\cdots ,7\right\}$ . Thus, the set of faces $ℱ=\left\{F_i,i\in 1,\cdots ,19\right\}$ is defined. It can be split into two subsets, the first one ${ℱ}^{\partial }=\left\{F_i,i\in 1,\cdots ,9\right\}$ boundary faces of $ℱ$ , depicted with dashed lines on Figure 2. The second one ${ℱ}^{\text{in}}=\left\{F_i,i\in 10,\cdots ,19\right\}$ interior faces of $ℱ$ . Figure 2(b) gives graphical representation for two elements $E_5$ and $E_7$ . Faces are also depicted with their normal vectors.

Figure 2. Example of a mesh.

Let two neighbouring elements $E_l$ and $E_r$ sharing one face $F\in ℱ$ . There are two traces of a function $v$ on $E_l$ ($v_l$ ) and on $E_r$ ($v_r$ ):

$v_l\left(x\right):=\underset{\epsilon \to 0^{-}}{\text{lim}}v\left(x+\epsilon n_F\right)\text{and}{v}_r\left(x\right):=\underset{\epsilon \to 0^{+}}{\text{lim}}v\left(x+\epsilon n_F\right),\forall x\in F.$ (7)

In addition, on any boundary faces $F\in {ℱ}^{\partial }$ the trace of $v$ is only defined on the left side of the face:

$v_l\left(x\right):=\underset{\epsilon \to 0^{-}}{\text{lim}}v\left(x+\epsilon n_F\right),\forall x\in F$ (8)

Using these trace definitions, one can define the jump and the average on any face of the mesh (as displayed in 1D on Figure 3). On an interior face $F\in {ℱ}^{\text{in}}$ , the jump and the average are respectively defined as:

$\forall x\in F,〚v〛\left(x\right):=v_r\left(x\right)-v_l\left(x\right)\text{and}\left\{\left|u\right|\right\}\left(x\right):=\frac{1}{2}\left(v_r\left(x\right)+v_l\left(x\right)\right).$ (9)

Moreover, on a boundary face $F\in {ℱ}^{\partial }$ , the jump and the average are respectively defined as:

$\forall x\in F,〚v〛\left(x\right):=v_l\left(x\right)\text{and}\left\{\left|u\right|\right\}\left(x\right):=v_l\left(x\right).$ (10)

Figure 3. Definition of the mean and jump operators for two elements $E_l$ and $E_r$ in 1D.

The solution of Problem () is sought in a subspace of the well-known broken Sobolev space, taken to be:

${\mathcal{V}}^p\left({ℰ}^n\right):=\left\{v\in L^2\left(\Omega \right)|{v|}_E\in {\mathbb{P}}^p\left(E\right),\forall E\in {ℰ}^n\right\}$ (11)

where ${\mathbb{P}}^p\left(E\right)$ stands for the set of polynomial functions of degree less than or equal to $p\in \mathbb{N}$ on $E$ . It is called the DG space. For more detailed and general definitions of this set, see [13].

3.2. Semi-Discrete Weak Formulation

Keeping in mind that

$\forall u,v\in {\mathcal{V}}^P\left({ℰ}^n\right),〚uv〛=〚u〛\left\{\left|v\right|\right\}+\left\{\left|u\right|\right\}〚v〛,$ (12)

assuming that the flux of RE is continuous at the interfaces of elements:

$\forall F\in ℱ,{〚\mathbb{K}\left(h-z\right)\nabla h\cdot n_F〛|}_F=0,$ (13)

the Neumann boundary condition arises naturally in the weak formulation, multiplying Problem () by a test function $\phi \in {\mathcal{V}}^p\left({ℰ}^n\right)$ and integrating on each element of $ℰ$ , we get

$\{\begin{array}{ll}{\displaystyle \sum _{E\in ℰ}{\displaystyle {\int }_E\theta \left(h-z\right)\phi \text{d}E}}+{\displaystyle \sum _{E\in ℰ}{\displaystyle {\int }_E\left(\mathbb{K}\left(h-z\right)\nabla h\right)\cdot \nabla \phi \text{d}E}}\hfill & \hfill \\ -{\displaystyle \sum _{F\in {ℱ}^{\text{in}}}{\displaystyle {\int }_F〚\left|\left(\mathbb{K}\left(h-z\right)\nabla h\right)\cdot n_F\right|〛〚\phi 〛\text{d}F}}-{\displaystyle \sum _{F\in {ℱ}^D}{\displaystyle {\int }_F\left(\mathbb{K}\left(h-z\right)\nabla h\right)\cdot n_F\phi \text{d}F}}\hfill & \hfill \\ +{\displaystyle \sum _{F\in {ℱ}^N}{\displaystyle {\int }_F{q}_N\phi \text{d}F}}=0,\hfill & \text{on}t\in \left(0,T\right)\hfill \\ {\displaystyle \sum _{E\in ℰ}{\displaystyle {\int }_{E}h\phi \text{d}E}}={\displaystyle \sum _{E\in ℰ}{\displaystyle {\int }_E{h}_0\phi \text{d}E}},\hfill & \hfill \\ h=h_D,\hfill & \text{on}{\Gamma }_D\times \left(0,T\right).\hfill \end{array}$ (14)

To enforce the continuity of the solution and the Dirichlet boundary condition, two penalty terms are added:

$J_I\left(h,\phi \right):={\displaystyle \sum }_{F\in {ℱ}^{\text{in}}}\frac{1}{2}\left(\frac{{\sigma }_E^{in}}{d_E}+\frac{{\sigma }_{E_r}^{in}}{d_{E_r}}\right){\displaystyle {\int }_F〚h〛〚\phi 〛\text{d}F}$ (15)

$J_D\left(h,\phi \right):={\displaystyle \sum }_{F\in {ℱ}^D}\frac{{\sigma }_E^{\partial }}{d_E}{\displaystyle {\int }_F\left(h-h_D\right)\phi \text{d}F}$ (16)

where, $J_I$ represents the penalization terms that constrain the continuity of the solution on the interior of the domain, and, $J_D$ for the Dirichlet boundary conditions. ${\sigma }_E^{in}$ and ${\sigma }_E^{\partial }$ are the penalization parameters for the interior and for the Dirichlet boundary condition where, we recall that, $d_E$ is the diameter of an element $E$ .

Remark. This method is known as the IIPG method [12] [14]-[16]. The role of these parameters is essential to ensure the convergence of the method and will be studied in Section 4 for the first time, up to our knowledge, in the nonlinear case. The linear case has been dealt in [17].

Using Equation (15) and Equation (16) in Equation (14), the semi-discrete nonlinear weak formulation of Problem () is, $\forall t\in T^n$ ,

$\{\begin{array}{l}\text{Find}h\in {\mathcal{V}}^p\left({ℰ}^n\right)\text{such}\text{that}:\\ m_n\left(\theta \left(h-z\right),\phi \right)+a_n\left(h,\phi ;h\right)=l_n\left(\phi \right),\forall \phi \in {\mathcal{V}}^p\left({ℰ}^n\right),\end{array}$ ()

where $m_n$ , $a_n$ , and $l_n$ are given by:

$m_n\left(q,\phi \right)={\displaystyle \sum _{E\in {ℰ}^n}{\displaystyle {\int }_{E}q\phi \text{d}E}}$ (17)

$\begin{array}{c}{a}_n\left(h,\phi ;h\right)={\displaystyle \sum _{E\in {ℰ}^n}{\displaystyle {\int }_E\left(\mathbb{K}\left(h-z\right)\nabla h\right)\cdot \nabla \phi \text{d}E}}\\ -{\displaystyle \sum _{F\in {ℱ}^{\text{in}}}{\displaystyle {\int }_F\left\{\left|\left(\mathbb{K}\left(h-z\right)\nabla h\right)\cdot n_F\right|\right\}〚\phi 〛\text{d}F}}\\ +{\displaystyle \sum }_{F\in {ℱ}^{\text{in}}}\frac{1}{2}\left(\frac{{\sigma }_E^{in}}{d_E}+\frac{{\sigma }_{E_r}^{in}}{d_{E_r}}\right){\displaystyle {\int }_F〚h〛〚\phi 〛\text{d}F}\\ -{\displaystyle \sum _{F\in {ℱ}^D}{\displaystyle {\int }_F\left(\mathbb{K}\left(h-z\right)\nabla h\right)\cdot n_F\phi \text{d}F}}+{\displaystyle \sum }_{F\in {ℱ}^D}\frac{{\sigma }_E^{\partial }}{d_E}{\displaystyle {\int }_{F}h\phi \text{d}F}\end{array}$ (18)

$l_n\left(\phi \right)={\displaystyle \sum }_{F\in {ℱ}^D}\frac{{\sigma }_E^{\partial }}{d_E}{\displaystyle {\int }_F{h}_D\phi \text{d}F}-{\displaystyle \sum _{F\in {ℱ}^N}{\displaystyle {\int }_F{q}_N\phi \text{d}F}}.$ (19)

3.3. Time Discretization

The aim of this section is to present the time discretisation through the implicit BDF method for Problem (). In the following, we make use of notation: $\forall n\in \mathbb{N}$ , $u^n\left(x\right):=u\left(x,t_n\right)$ , for any function $u\in L^2\left(\Omega \times \left(0,T\right)\right)$ . Let us recall that the time step is defined by $\Delta t^n=t^{n+1}-t^n$ and the time interval by $T^n=\left[t^n,t^{n+1}\right]$ .

Due to their stability properties, the BDF methods are commonly used to solve stiff differential equations such as Problem (). These linear multi-step methods allow to construct time approximation up to order $q\le 6$ . The analysis of these methods can be found in [18]. The 1-step BDF method corresponds to the classical backward Euler scheme. BDF methods have been used in [19] [20] up to 6th-order. BDF methods are well-known to balance space and time errors and particularly well-designed in combination with DG methods. BDF methods can be constructed both with a constant time step [18] or a variable [21]. The case of variable time step is more pertinent for Problem () concerned. The method of order $q$ is derived from the Newton interpolation polynomial of degree $q$ , which interpolates $h^j$ at time $t^j$ for $j=n+1,\cdots ,n+1-q$ , using the method of divided difference.

The backward divided difference for a given function $y$ is defined by a recursive division process:

$\{\begin{array}{l}{\delta }^0{y}^{n+1}=\left[y^{n+1}\right]=y^{n+1},\hfill \\ {\delta }^1{y}^{n+1}=\left[y^{n+1},y^n\right]=\frac{{\delta }^0{y}^{n+1}-{\delta }^0{y}^n}{\Delta t^n}=\frac{y^{n+1}-y^n}{\Delta t^n},\hfill \\ {\delta }^2{y}^{n+1}=\left[y^{n+1},y^n,y^{n-1}\right]=\frac{{\delta }^1{y}^{n+1}-{\delta }^1{y}^n}{\Delta t^n+\Delta t^{n-1}}=\frac{\frac{y^{n+1}-y^n}{\Delta t^n}-\frac{y^n-y^{n-1}}{\Delta t^{n-1}}}{\Delta t^n+\Delta t^{n-1}},\hfill \\ ⋮\hfill \\ {\delta }^j{y}^{n+1}=\left[y^{n+1},y^n,\cdots ,y^{n+1-j}\right]=\frac{{\delta }^{j-1}{y}^{n+1}-{\delta }^{j-1}{y}^n}{{{\displaystyle \sum }}_{k=0}^{j-1}\Delta t^{n-k}}.\hfill \end{array}$ (20)

For a given ode, for instance $\frac{\text{d}u}{\text{d}t}=f\left(u,t\right)$ with initial condition, the implicit BDF method of order $q$ is given by: $\begin{array}{l}{\displaystyle \sum }_{j=1}^q\left({\displaystyle \prod }_{k=1}^{j-1}\left({\displaystyle \sum }_{l=0}^{k-1}\Delta t^{n-l}\right)\right){\delta }^j{u}^{n+1}={\displaystyle \sum }_{j=0}^q{\alpha }_{q,j}{u}^{n+1-j}=f\left(u^{n+1},t^{n+1}\right),\\ \iff {\alpha }_{q,0}{u}^{n+1}-f\left(u^{n+1},t^{n+1}\right)=-{\displaystyle \sum }_{j=0}^{q-1}{\alpha }_{q,j+1}{u}^{n-j}\end{array}$ where ${\alpha }_{q,j}$ are the linear combination coefficients obtained from the divided differences of $u$ . For instance, for the 2-order BDF method, the coefficients are:

$\{\begin{array}{l}{\alpha }_{2,0}=\frac{1}{\Delta t^n}+\frac{1}{\Delta t^n+\Delta t^{n-1}},\\ {\alpha }_{2,1}=-\frac{1}{\Delta t^n}-\frac{1}{\Delta t^n+\Delta t^{n-1}}-\frac{\Delta t^n}{\Delta t^{n-1}\left(\Delta t^n+\Delta t^{n-1}\right)},\\ {\alpha }_{2,2}=\frac{\Delta t^n}{\Delta t^{n-1}\left(\Delta t^n+\Delta t^{n-1}\right)}.\end{array}$ (21)

Remark. (Stability.) BDF methods of order 1 and 2 are $A$ -stable, and $L$ -stable [22]. BDF methods of order 3 to 6 are $A\left(\alpha \right)$ -stable where $\alpha$ decreases as the order increases [23]. BDF methods of order $q>6$ are unconditionally unstable. The use of variable time steps is recommended to enhance the stability of the method. In practical applications, variations in time step sizes are limited by an upper bound known as the swing factor to ensure stability and robustness Table 2 (see [24]). In the following, swing factors are used.

Table 2. Maximum swing $\Delta t^{n+2}/\Delta t^{n+1}$ for BDF methods with variable time steps.

Order $q$	1	2	3	4	5	6
Maximum swing $\Delta t^{n+2}/\Delta t^{n+1}$	-	2.6	1.9	1.5	1.2	1.05

Applying the BDF method to Problem (), we get

$\{\begin{array}{l}\text{Find}\text{as}\text{equence}\text{of}{\left(h^n\right)}_{0\le n\le N}\in {\mathcal{V}}^p\left({ℰ}^n\right)\text{such}\text{that}:\\ m_n\left({\frac{\partial \theta \left(\psi \right)}{\partial \psi }|}_{{\psi }^{n+1}}{\displaystyle \sum }_{j=0}^q{\alpha }_{q,j}{h}^{n+1-j},\phi \right)+a_n\left(h^{n+1},\phi ;h^{n+1}\right)=L_n\left(\phi \right),\forall \phi \in {\mathcal{V}}^p\left({ℰ}^n\right).\end{array}$ ()

where $m_n$ , $a_n$ and $l_n$ are given, respectively, by Equation (17), Equation (18), Equation (19) with $\psi =h-z$ .

The time integration method needs an initialization step to compute the solution for further time steps. The initialization uses the prescribed initial condition to start the first-time step. A direct and simple way is to write the corresponding discontinuous weak formulation:

$\text{Find}{h}^0\in {\mathcal{V}}^p\left({ℰ}^0\right)\text{such}\text{that}:m_0\left(h^0,\phi \right)=f_0\left(\phi \right),$ (22)

where $m_0$ is defined by Equation (17) and $f_0$ is the linear form defined by:

$f_0\left(\phi \right)={\displaystyle \sum _{E\in {ℰ}^0}{\displaystyle {\int }_E{h}_0\phi \text{d}E}},\forall \phi \in {\mathcal{V}}^p\left({ℰ}^0\right).$ (23)

3.4. Nonlinear Iterative Process

Problem () being nonlinear, several iterative methods can be used such as the Newton-Raphson method or the classical first-order fixed-point method Picard’s method. Due to the strong nonlinearities of the constitutive laws Equation (2) and Equation (4) (see also Remark 1), the convergence of the iterative methods may fail [25] [26]. We will see in Section 4 that in the case of IIPG methods one can enhance the convergence of the iterative methods, at least in the case of a Picard’s fixed-point method, whenever the penalization terms Equation (15) and Equation (16) are well-chosen. Therefore, in what follows, we present the Picard’s fixed-point method for Problem ().

Remark. (Choice of the Picard linearization.) Although Newton-Raphson iterations may offer quadratic convergence, we adopted a Picard fixed-point linearization for robustness in strongly nonlinear configurations such as the Vauclin infiltration case. In preliminary tests, Newton iterations often diverged without regularization of the hydraulic functions, whereas the Picard approach provided stable convergence at a moderate computational cost. Similar observations have been reported in [7] [27], highlighting that Picard iterations remain preferable when the Jacobian varies sharply near saturation thresholds.

Linearization of Problem () is done by a Picards’ iterative procedure. For $k=0,\cdots$ , the problem is:

$\{\begin{array}{l}\text{For}\text{a}\text{given}{h}^{n+1,k}\in {\mathcal{V}}^p\left({ℰ}^n\right),\text{find}{h}^{n+1,k+1}\in {\mathcal{V}}^p\left({ℰ}^n\right)\text{such}\text{that},\forall \phi \in {\mathcal{V}}^p\left({ℰ}^n\right):\\ m_n\left({\frac{\partial \theta \left(\psi \right)}{\partial \psi }|}_{{\psi }^{n+1,k}}{\alpha }_{q,0}{h}^{n+1,k+1},\phi \right)+a_n\left(h^{n+1,k+1},\phi ;h^{n+1,k}\right)\\ =l_n\left(\phi \right)-m_n\left({\frac{\partial \theta \left(\psi \right)}{\partial \psi }|}_{{\psi }^{n+1,k}}{\displaystyle \sum }_{j=0}^{q-1}{\alpha }_{q,j+1}{h}^{n-j},\phi \right).\end{array}$ ()

where $m_n$ , $a_n$ and $l_n$ are given, respectively, by Equation (17), Equation (18), Equation (19) with $\psi =h-z$ . $h^{n-j}$ stands for the solution at the rank $k$ of the iterative process (see Figure 4).

Figure 4. Scheme of the general proof.

The global algorithm of the Picard’s fixed-point iteration, for a positive $n$ , is:

1) Start with an initial guess $h^{n+1,0}$ ;

2) Compute the solution of Problem () with $h^{n+1,0}$ to get $h^{n+1,1}$ ;

3) Start again with $h^{n+1,1}$ ;

4) Compute the solution of Problem () with $h^{n+1,k}$ to get $h^{n+1,k+1}$ ;

5) Start again with $h^{n+1,k+1}$ until the stopping criteria are satisfied;

6) Set $h^{n+1}=h^{n+1,k+1}$ .

The stopping criterion is one important choice in determining accuracy for a nonlinear iterative process. For RE, the stopping criterion can be specified in terms of absolute error for pressure head or water content between two successive iterations [12]. For this study, we have used: $\frac{\Vert r_n\left(h,\phi \right)\Vert }{\Vert a_n\left(h,\phi \right)\Vert }<{\epsilon }_1$ and $\frac{\Vert {\delta }_k\Vert }{\Vert h^k\Vert }<{\epsilon }_2$ , where ${\delta }_k=h^k-h^{k-1}$ and $r_n\left(h,\phi \right)=m_n\left(h,\phi ;h\right)+a_n\left(h,\phi ;h\right)-l_n\left(\phi \right)$ . ${\epsilon }_1$ and ${\epsilon }_2$ are user-defined tolerances. These two criteria are relative and independent of the characteristic quantities of the problem.

3.5. Adaptive Time Stepping

Time adaptation is motivated by the convergence of the nonlinear solver. On one hand, transient simulations have difficulties to converge if the time step is too large but, on the other hand, shorter time steps mean more time steps and so, a longer computational time. That is the reason why time adaptation is very attractive and common for Richards’ equation. Different strategies can be used to adjust the time step [28]-[30], either heuristic and mainly based on convergence performance of the nonlinear solver or rational and based on error control. The latter ones are generally more efficient but heuristic methods remain a relevant approach due to their simplicity.

In this study, the time step is adjusted heuristically based on the number of iterations $N_{it}$ from the nonlinear solver, as discussed in [29] [31]. The size of the time step directly influences the convergence of the solver. The simulations start with a time step $\text{Δ}{t}^0$ , and subsequent time steps are calculated according to the following rule: the time step remains unchanged if convergence is achieved between $m_{it}$ and $M_{it}$ nonlinear iterations; it is increased by an amplification factor ${\lambda }_{amp}>1$ if convergence is achieved in fewer than $m_{it}$ nonlinear iterations; and it is decreased by a reduction factor ${\lambda }_{red}<1$ if convergence requires more than $M_{it}$ nonlinear iterations. If convergence fails due to solver issues (poor initial guess, bad condition number) or exceeds a prescribed maximum bound $W_{it}$ , the time step is recalculated using a reduced step size (${\lambda }_{red}<1$ ). The calculation of the next time step $\text{Δ}{t}^{n+1}$ from the previous one $\text{Δ}{t}^n$ follows this time-stepping scheme:

$\{\begin{array}{l}\text{Δ}{t}^{n+1}=\{\begin{array}{ll}{\lambda }_{amp}\text{Δ}{t}^n\hfill & \text{if}{N}_{it}\le m_{it},\hfill \\ \text{Δ}{t}^n\hfill & \text{if}{m}_{it}<N_{it}\le M_{it},\hfill \\ {\lambda }_{red}\text{Δ}{t}^n\hfill & \text{if}{M}_{it}<N_{it}\le W_{it},\hfill \end{array}\\ \text{Δ}{t}^n={\lambda }_{red}\text{Δ}{t}^n\text{if}{N}_{it}>W_{it}\text{or}\text{if}\text{the}\text{solver}\text{has}\text{failed}\left(\text{time}\text{step}\text{is}\text{started}\text{again}\right),\end{array}$ (24)

where $N_{it}$ is the number of nonlinear iterations.

Remark. By studying the full-time-dependent problem, as done in Section 4 in the case of the steady problem, the time step can be adjusted automatically and this work is in progress.

Remark. In the numerical code RIVAGE, Adaptive Mesh Refinement can be also employed. We refer to [7] [12] [32]-[34] for more details.

4. Theoretical Study and Estimation of the Optimal Penalization Parameters

In this section, we present the main result of this work, namely, the way to get a convergent iterative scheme by constructing a robust method to compute automatically the penalization parameters (see Equation (15) and Equation (16)). This is achieved by studying the theoretical properties and convergence of the solution of the discrete problem Problem () to the mathematical problem Problem (). To this end and for the sake of simplicity, we will consider a toy model similar to the stationary RE for which we study, as depicted in Section 4.

1) The existence and uniqueness of the weak solution to the nonlinear problem in Section 4.2.

2) The existence and uniqueness of the weak solution to the discrete linearized problem in Section 4.3.

3) The method to compute optimal penalization parameters to ensure the convergence of the nonlinear solver at the discrete level in Section 4.4.

4) The convergence of the discrete linearized weak problem to the continuous linearized weak problem in Section 4.5.

Proofs of this section are given in Appendix and can be easily extended to several space dimensions. However, since the computations are rather technical to get the optimal penalization parameters in the two-dimensional case, for the sake of completeness, the 2D case for the existence and uniqueness of the weak solution to the discrete linearized problem is considered in Section 4.3. We will see that the construction of the optimal penalization parameters is essentially based on the constants appearing in the discrete continuity and the discrete coercivity of the operator.

4.1. Toy Model

Let us consider the following toy problem () on the interval $\Omega =\left[a,b\right]\subset ℝ$ :

$\begin{array}{l}\text{For}\text{a}\text{given}f\in L^2\left(\Omega \right),\text{find}u\left(x\right):\Omega \to ℝ\text{such}\text{that}:\\ \{\begin{array}{ll}-{\left(A\left(x,u,u^{\prime }\right)\right)}^{\prime }=f,\hfill & \text{in}\Omega \hfill \\ u=0,\hfill & \text{on}\partial \Omega \hfill \end{array}\end{array}$ ()

with $A\left(x,s,\xi \right)=K\left(x,s\right)\xi$ where the real function $K$ intends to mimick the properties of $\mathbb{K}$ (Equation (2)). Following [35] and in view of the properties of $\mathbb{K}$ (Equation (2)), assuming that

$\{\begin{array}{lll}\exists K_0,K_1\in {ℝ}_{+}^{*},\hfill & K_0\le K\left(x,\overline{u}\right)\le K_1,\hfill & \forall x\in \Omega ,\forall \overline{u}\in ℝ\hfill \\ \exists K_{lip}\in {ℝ}_{+},\hfill & \left|K\left(x,{\overline{u}}_1\right)-K\left(x,{\overline{u}}_2\right)\right|\le K_{lip}\left|{\overline{u}}_1-{\overline{u}}_2\right|,\hfill & \forall x\in \Omega ,\forall \left({\overline{u}}_1,{\overline{u}}_2\right)\in {ℝ}^2\hfill \end{array},$ ($\mathscr{H}1$ )

we deduce that $A$ is straightforwardly a Carathéodory function, which we recall hereafter,

$\begin{array}{lll}\left(1\right)\exists \alpha >0\hfill & \text{s}.\text{t}.\hfill & \left(A\left(x,s,\xi \right)-A\left(x,s,0\right)\right)\xi \ge \alpha {\left|\xi \right|}^2,\hfill \\ \left(2\right)\exists \beta >0,\exists h\in L^2\left(\Omega \right)\hfill & \text{s}.\text{t}.\hfill & \left|A\left(x,s,\xi \right)\right|\le \beta \left(h\left(x\right)+\left|s\right|+\left|\xi \right|\right),\hfill \\ \left(3\right)\exists \gamma >0\hfill & \text{s}.\text{t}.\hfill & \left(A\left(x,s,\xi \right)-A\left(x,s,\eta \right)\right)\left(\xi -\eta \right)\ge \gamma |\xi -\eta {|}^2,\hfill \\ \left(4\right)\exists \delta >0,\exists h\in L^2\left(\Omega \right)\hfill & \text{s}.\text{t}.\hfill & \left|A\left(x,s,\xi \right)-A\left(x,t,\xi \right)\right|\le \delta \left|s-t\right|\left(h\left(x\right)+\left|\xi \right|+\left|s\right|+\left|t\right|\right).\hfill \end{array}$ ($\mathscr{H}2$ )

This problem can be cast into the weak formulation by multiplying by a test function $v\in H_0^1\left(\Omega \right)$ and integrating over Ω:

$\text{Find}u\in H_0^1\left(\Omega \right)\text{such}\text{that}:a\left(u,v\right)=l\left(v\right),\forall v\in H_0^1\left(\Omega \right)$ ($\mathcal{W}$ )

where

$a\left(u,v\right)={\displaystyle {\int }_{\Omega }K\left(x,u\right)u^{\prime }{v}^{\prime }\text{d}x},l\left(v\right)={\displaystyle {\int }_{\Omega }fv\text{d}x}.$

Problem () being nonlinear, we use the Picard’s iterations method as in Problem () to get

$\{\begin{array}{l}\text{For}\text{a}\text{given}\overline{u}\in L^{\text{2}}\left(\Omega \right),\text{find}u\in H_0^1\left(\Omega \right)\text{such}\text{that}:\\ \tilde{a}\left(u,v;\overline{u}\right)=l\left(v\right),\forall v\in H_0^1\left(\Omega \right)\end{array}$($\tilde{\mathcal{W}}$ )

with

$\tilde{a}\left(u,v;\overline{u}\right)={\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right)u^{\prime }{v}^{\prime }\text{d}x.$(25)

Given ${\overline{u}}^0$ , we solve the Problem ($\tilde{\mathcal{W}}$ ) with $\overline{u}={\overline{u}}^0$ to obtain $u^1$ . Then, we solve the Problem ($\tilde{\mathcal{W}}$ ) with $\overline{u}={\overline{u}}^1$ to obtain $u^2$ and so on. The sequence of solutions of the linearized problem is denoted by ${\left(u^n\right)}_{n\in \mathbb{N}}$ and its limit when $n$ goes to infinity is expected to be the solution to the nonlinear Problem ($\mathcal{W}$ ). In the following, we note $u^{n+1}=T\left(u^n\right)$ the fixed point.

4.2. Existence and Uniqueness of the Weak Solution to the Nonlinear Problem ($\mathcal{W}$ )

The first step is to show that Problem ($\mathcal{W}$ ) has a unique solution in $H_0^1\left(\text{Ω}\right)$ . The existence of solution of Problem ($\mathcal{W}$ ) can be achieved by using the Schauder fixed-point theorem to the operator $T$ while the uniqueness can be obtained through the technique proposed in [35].

Thus, we have

Lemma 1. (Existence of a solution to Problem ($\mathcal{W}$ ).) Under Hypothesis ($\mathscr{H}1$ ), $\exists u\in H_0^1\left(\Omega \right)$ ; $T\left(u\right)=u$ .

Then, one can obtain uniqueness through the following result

Lemma 2. (Uniqueness of the solution to Problem ($\mathcal{W}$ ).) Under Hypothesis ($\mathscr{H}1$ ), the solution $u\in H_0^1\left(\Omega \right)$ of Problem ($\mathcal{W}$ ) is unique.

These results hold for the dimension $d\le 3$ and the proofs are rather classical and left to the reader.

4.3. Existence and Uniqueness of the Weak Solution to the Discrete Linearized Problem ($\tilde{\mathcal{W}}$ )

One-dimensional case

To solve numerically Problem ($\tilde{\mathcal{W}}$ ), we use DG methods as in Section 3. Let $0=x_0<\cdots <x_N=1$ be a partition ${ℰ}_h$ of Ω (see Figure 5) and denote $I_n=\left[x_n,x_{n+1}\right]$ a sub-interval. The size of a sub-interval is defined as $\left|I_n\right|:=h=\frac{1}{N}$ , $\forall n\in \left\{0,\cdots ,N-1\right\}$ with $N$ —the number of elements in the partition. The solution is sought in the DG space ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ defined as:

${\mathcal{V}}_0^p\left({ℰ}_h\right)=\left\{v\in L^{\text{2}}\left(\Omega \right)|{v|}_{\partial \Omega }=0;{v|}_{I_n}\in {\mathbb{P}}^p\left(I_n\right),\forall I_n\in {ℰ}_h\right\}\subseteq L^{\text{2}}\left(\Omega \right)$ (26)

Figure 5. Representation of ${ℰ}_h$ in the one-dimension case.

As in Section 3, we define

$v\left(x_n^{+}\right)=\underset{\begin{array}{c}ϵ\to 0\\ ϵ>0\end{array}}{\lim }v\left(x_n+ϵ\right),v\left(x_n^{-}\right)=\underset{\begin{array}{c}ϵ\to 0\\ ϵ>0\end{array}}{\lim }v\left(x_n-ϵ\right),$ (27)

${〚v〛}_{x_n}=v\left(x_n^{-}\right)-v\left(x_n^{+}\right),{\left\{\left|v\right|\right\}}_{x_n}=\frac{1}{2}\left(v\left(x_n^{-}\right)+v\left(x_n^{+}\right)\right),\forall n\in \left\{1,\cdots ,N-1\right\},$ (28)

and

${〚v〛}_{x_0}=-v\left(x_0^{+}\right),{\left\{\left|v\right|\right\}}_{x_0}=v\left(x_0^{+}\right),{〚v〛}_{x_N}=v\left(x_N^{-}\right),{\left\{\left|v\right|\right\}}_{x_N}=v\left(x_N^{-}\right).$ (29)

The DG space ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ is associated with the norm:

${\Vert v\Vert }^2={\displaystyle \sum }_{n=0}^{N-1}{\Vert v^{\prime }\Vert }_{I_n}^2+{\displaystyle \sum }_{n=0}^N\frac{1}{h}{〚v〛}_{x_n}^2={\displaystyle \sum }_{n=0}^{N-1}{\Vert v^{\prime }\Vert }_{I_n}^2+{\left|v\right|}_J^2$ (30)

where ${\Vert \cdot \Vert }_{I_n}$ is the usual norm $L^2\left(I_n\right)$ and ${\left|v\right|}_J^2:={\displaystyle {\sum }_{n=0}^N\frac{1}{h}{〚v〛}_{x_n}^2}$ is the jump semi-norm. With this definition of the norm, jumps are controlled. One can observe that $\Vert \cdot \Vert$ is a norm on ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ . One can note that ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ is a complete Banach space, i.e., a complete normed vector space for $\Vert \cdot \Vert$ . Lastly, the concept of broken gradient is introduced to specify when only the regular part of the gradient is considered. The broken gradient ${\nabla }_h:{\mathcal{V}}_0^p\left({ℰ}_h\right)\to L^2\left(\Omega \right)$ is defined that, for all $v\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ ,

$\forall E\in {ℰ}_h{}_h,{\left({\nabla }_{h}v\right)|}_E:=\nabla \left({v|}_E\right).$ (31)

The linearized weak formulation Problem ($\tilde{\mathcal{W}}$ ) can be discretized using the IIPG formulation as in Section 3 to get

$\{\begin{array}{l}\text{For}\text{a}\text{given}\overline{u}\in {\mathcal{V}}_0^p\left({ℰ}_h\right),\text{find}{u}_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right)\text{such}\text{that}:\\ {\tilde{a}}_h\left(u_h,v_h;\overline{u}\right)=l_h\left(v_h\right),\forall v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right)\end{array}$(${\tilde{\mathcal{W}}}_h$ )

with

$\begin{array}{c}{\tilde{a}}_h\left(u_h,v_h,\overline{u}\right)={\displaystyle \sum }_{n=0}^{N-1}{\displaystyle {\int }_{I_n}K}\left(x,\overline{u}\right){u^{\prime }}_h{v^{\prime }}_h\text{d}x-{\displaystyle \sum }_{n=0}^N{\left\{\left|K\left(x,\overline{u}\right)u^{\prime }\right|\right\}}_h{}_{x_n}{〚v_h〛}_{x_n}\\ +\frac{{\sigma }_0}{h}{〚u_h〛}_{x_0}{〚v_h〛}_{x_0}+{\displaystyle \sum }_{n=1}^{N-1}\frac{{\sigma }_{n-1}+{\sigma }_n}{2h}{〚u_h〛}_{x_n}{〚v_h〛}_{x_n}\\ +\frac{{\sigma }_N}{h}{〚u_h〛}_{x_N}{〚v_h〛}_{x_N},\end{array}$(32)

$l_h\left(v_h\right)={\displaystyle {\int }_{\Omega }f{v}_h\text{d}x}.$

At the discrete level, one can write Hypothesis ($\mathscr{H}1$ ) as follows: for all $n\in \left\{0,\cdots ,N-1\right\}$ :

$\{\begin{array}{ll}\exists K_0^{\left(n\right)},K_1^{\left(n\right)}\in {ℝ}_{+}^{\text{*}},\forall x\in I_n,\forall \overline{u}\in ℝ,\hfill & K_0^{\left(n\right)}\le K\left(x,\overline{u}\right)\le K_1^{\left(n\right)};\hfill \\ \exists K_{lip}^{\left(n\right)}\in {ℝ}_{+},\forall x\in I_n,\forall \left({\overline{u}}_1,{\overline{u}}_2\right)\in {ℝ}^2,\hfill & \left|K\left(x,{\overline{u}}_1\right)-K\left(x,{\overline{u}}_2\right)\right|\le K_{lip}\left|{\overline{u}}_1-{\overline{u}}_2\right|\hfill \end{array}$ (${\mathscr{H}}_n$ )

where

$K_1:=\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\min }{K}_1^{\left(n\right)},K_0:=\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\max }{K}_0^{\left(n\right)}\text{and}{K}_{lip}=\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\max }{K}_{lip}^{\left(n\right)}.$ (33)

Existence and unicity for the solution to Problem (${\tilde{\mathcal{W}}}_h$ ) is obtained using the Lax-Milgram theorem. We have the following result.

Theorem 3. (Existence and uniqueness of the weak solution to the discrete linearized Problem (${\tilde{\mathcal{W}}}_h$ ).) Under Hypothesis (${\mathscr{H}}_n$ ) for all $n$ , for a given $\overline{u}\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ , then $\exists !u\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ such that ${\tilde{a}}_h\left(u_h,v_h;\overline{u}\right)=l_h\left(v_h\right)$ , $\forall v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ .

This existence and uniqueness result is obtained thanks to the below-following lemmas.

Lemma 4. (Discrete coercivity of ${\tilde{a}}_h$ .) Under Hypothesis (${\mathscr{H}}_n$ ) for all $n$ , for any vector of positive numbers $ϵ={\left({\epsilon }^{\left(n\right)}\right)}_{n=0,\cdots ,N-1}$ , there exists a constant $C^{*}\left(ϵ\right)>0$ such that

$\forall u_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right),{\tilde{a}}_h\left(u_h,u_h;\overline{u}\right)\ge C^{\text{*}}\left(ϵ\right){\Vert u_h\Vert }^2$

if

$\{\begin{array}{ll}{\epsilon }^{\left(n\right)}<2,\hfill & \forall n\in \left\{0,\cdots ,N-1\right\}\hfill \\ {\sigma }_n>{\sigma }_n^{\text{*}},\hfill & \forall n\in \left\{1,\cdots ,N-1\right\}\hfill \\ {\sigma }_0>{\sigma }_0^{\text{*}}\hfill & \hfill \\ {\sigma }_N>{\sigma }_N^{\text{*}}\hfill & \hfill \end{array}\text{with}\{\begin{array}{l}{\sigma }_n^{\text{*}}=\frac{{\left(K_1^{\left(n\right)}{C}_{\text{tr},p-1}^{\left(n\right)}\right)}^2}{2{\epsilon }^{\left(n\right)}{K}_0^{\left(n\right)}},\forall n\in \left\{1,\cdots ,N-1\right\}\hfill \\ {\sigma }_0^{\text{*}}=\frac{{\left(K_1^{\left(0\right)}{C}_{\text{tr},p-1}^{\left(0\right)}\right)}^2}{{\epsilon }^{\left(0\right)}{K}_0^{\left(0\right)}}\hfill \\ {\sigma }_N^{\text{*}}=\frac{{\left(K_1^{\left(N-1\right)}{C}_{\text{tr},p-1}^{\left(N-1\right)}\right)}^2}{{\epsilon }^{\left(N-1\right)}{K}_0^{\left(N-1\right)}}\hfill \end{array}$ (34)

and

$C^{\text{*}}\left(ϵ\right)=\text{min}\left\{\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{min}}\left(K_0^{\left(n\right)}\left(1-\frac{{\epsilon }^{\left(n\right)}}{2}\right)\right),{\sigma }_0-{\sigma }_0^{\text{*}},{\sigma }_N-{\sigma }_N^{\text{*}},\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{min}}\left(\frac{{\sigma }_n-{\sigma }_n^{\text{*}}}{2}\right)\right\}.$ (35)

Lemma 5. (Discrete continuity of ${\tilde{a}}_h$ .) Under Hypothesis (${\mathscr{H}}_n$ ) for all $n$ , for any vector of positive numbers $ϵ={\left({\epsilon }^{\left(n\right)}\right)}_{n=0,\cdots ,N-1}$ , there exists a constant $\tilde{C}\left(ϵ\right)>0$ such that

$\forall u_h,v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right),\left|{\tilde{a}}_h\left(u_h,v_h;\overline{u}\right)\right|\le \tilde{C}\left(ϵ\right)\Vert u_h\Vert \Vert v_h\Vert$

where

$\begin{array}{c}\tilde{C}\left(ϵ\right)=\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{max}}\left(K_1^{\left(n\right)}\right)+\sqrt{\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{max}}\left(2{\epsilon }^{\left(n\right)}{K}_1^{\left(n\right)}\right)\text{max}\left({\sigma }_0^{\text{*}},{\sigma }_N^{\text{*}},\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{max}}\left(\frac{{\sigma }_n^{\text{*}}}{2}\right)\right)}\\ +\max \left({\sigma }_0,{\sigma }_N,\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{max}}\left(\frac{{\sigma }_n}{2}\right)\right).\end{array}$(36)

Lemma 6. (Discrete continuity of $l_h$ .) There exists a constant $B>0$ such that $\forall v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ , $\Vert l_h\left(v_h\right)\Vert \le B\Vert v_h\Vert$ .

Remark. Trace constantly involved in bounds for penalization parameters are a function of the polynomial degree $p$ , the type of polynomial basis used. In the one-dimensional case, with an orthonormal basis and for $u\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ , the trace constant for $I_n$ is given by:

$C_{\text{tr},p}^{\left(n\right)}:=p+1.$ (37)

Proofs of Lemmas 4, 5, and 6 can be found in Appendix. The proof of Theorem 3 is a straightforward application of the Lax-Milgram theorem and is left to the reader.

Two-dimensional case

We propose to extend the previous results to the dimension 2. Let us consider the two-dimensional extension of Problem (${\tilde{\mathcal{W}}}_h$ )

$\{\begin{array}{l}\text{For}\text{a}\text{given}\overline{u}\in {\mathcal{V}}_0^p\left({ℰ}_h\right),\text{find}{u}_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right)\text{such}\text{that},\forall v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right):\\ {\tilde{a}}_h\left(u_h,v_h;\overline{u}\right)=l_h\left(v_h\right).\end{array}$(${\tilde{\mathcal{W}}}_h^2$ )

where

$\begin{array}{c}\tilde{a}\left(u_h,v_h;\overline{u}\right)={\displaystyle \sum _{E\in {ℰ}^n}{\displaystyle {\int }_E\left(\mathbb{K}\left(x,\overline{u}\right)\nabla u_h\right)\cdot \nabla v_h\text{d}E}}-{\displaystyle \sum _{E\in ℱ}{\displaystyle {\int }_F\left\{\left|\left(\mathbb{K}\left(x,\overline{u}\right)\nabla u_h\right)\cdot n_F\right|\right\}〚v_h〛\text{d}F}}\\ +{\displaystyle \sum _{E\in {ℱ}^{\text{in}}}\frac{1}{2}}\left(\frac{{\sigma }_E^{in}}{d_E}+\frac{{\sigma }_{E_r}^{in}}{d_{E_r}}\right){\displaystyle {\int }_F〚u_h〛〚v_h〛\text{d}F}+{\displaystyle \sum }_{E\in {ℱ}^D}\frac{{\sigma }_E^{\partial }}{d_E}{\displaystyle {\int }_F{u}_h{v}_h\text{d}F}\end{array}$

$l_h\left(v_h\right)={\displaystyle {\int }_{\Omega }f{v}_h\text{d}x}.$

The two-dimensional version of the discrete hypothesis on $\mathbb{K}$ is given by: For all $E\in ℰ$ :

$\{\exists K_0^E,K_1^E\in {ℝ}_{+}^{\text{*}},\forall x\in E,\forall \overline{u}\in ℝ,K_0^E\le {\Vert \mathbb{K}\left(x,\overline{u}\right)\Vert }_2\le K_1^E;$(${\mathscr{H}}_E^2$ )

with ${\Vert \mathbb{K}\Vert }_2={\max }_{i=1,2}\left({\mathbb{K}}_{ii}\right)$ . In addition, $K_1={\max }_{E\in ℰ}{K}_1^E$ and $K_0={\min }_{E\in ℰ}{K}_0^E$ denotes global bound of $\mathbb{K}$ .

The DG space is associated with the following norm:

${\Vert v\Vert }^2:={\displaystyle \sum }_{E\in ℰ}{\Vert v\Vert }_E^2+{\displaystyle \sum }_{E\in {ℱ}^{\text{in}}}\left(\frac{1}{d_E}+\frac{1}{d_{E_r}}\right){\Vert 〚v〛\Vert }_F^2+{\displaystyle \sum }_{E\in {ℱ}^{\partial }}\frac{1}{d_E}{\Vert 〚v〛\Vert }_F^2={\displaystyle \sum }_{E\in ℰ}{\Vert v\Vert }_E^2+{\left|v\right|}_J^2$ (38)

where ${\Vert v\Vert }_E^2$ is the usual $L^2$ norm on $E$ , ${\Vert v\Vert }_F^2$ is the $L^2$ norm on $F$ and ${\left|v\right|}_J^2$ is the jump semi-norm. This norm has the same characteristics as in the one-dimensional case. We obtain the following result.

Theorem 7. (Existence and uniqueness of the weak solution to the discrete linearized Problem (${\tilde{\mathcal{W}}}_h$ ).) If $K$ satisfies Hypothesis (${\mathscr{H}}_E^2$ ) for all $E\in ℰ$ and for a given $\overline{u}\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ , then $\exists !u\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ such that ${\tilde{a}}_h\left(u_h,v_h;\overline{u}\right)=l_h\left(v_h\right)$ , $\forall v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ .

As before, this result is a consequence of the Lax-Milgram theorem through the following lemmas:

Lemma 8. (Discrete coercivity of ${\tilde{a}}_h$ .) If $K$ satisfies Hypothesis (${\mathscr{H}}_E^2$ ) for all $E\in ℰ$ and for any vector of positive numbers $ϵ={\left({\epsilon }^E\right)}_{E\in ℰ}$ , there exists a constant $C^{*}\left(ϵ\right)>0$ such that

$\forall u_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right),{\tilde{a}}_h\left(u_h,u_h;\overline{u}\right)\ge C^{\text{*}}\left(ϵ\right){\Vert u_h\Vert }^2$

if

$\{\begin{array}{ll}{\epsilon }^E<2,\hfill & \forall E\in ℰ\hfill \\ {\sigma }_E^{in}>{\sigma }_E^{in,\ast }\text{and}{\sigma }_{E_r}^{in}>{\sigma }_{E_r}^{in,\text{*}},\hfill & \forall F\in {ℱ}^{\text{in}}\hfill \\ {\sigma }_E^{\partial }>{\sigma }_E^{\partial ,\ast },\hfill & \forall F\in {ℱ}^{\partial }\hfill \end{array}\text{with}\forall E\in ℰ\{\begin{array}{l}{\sigma }_E^{in,\ast }=\frac{D^E{\left(K_1^E{C}_{\text{tr},p-1}^E\right)}^2}{4{\epsilon }^E{K}_0^E}\\ {\sigma }_E^{\partial ,\ast }=\frac{D^E{\left(K_1^E{C}_{\text{tr},p-1}^E\right)}^2}{2{\epsilon }^E{K}_0^E}\end{array}$ (39)

and $D^E$ is the number of edges of the element $E$ . Moreover

$C^{\text{*}}\left(ϵ\right)=\text{min}\left\{\underset{E\in ℰ}{\text{min}}\left(K_0^E\left(1-\frac{{\epsilon }^E}{2}\right)\right),\underset{E\in ℰ}{\text{min}}\left(\frac{{\sigma }_E^{in}-{\sigma }_E^{in,\text{*}}}{2}\right),\underset{F\in {ℱ}^{\partial }}{\text{min}}\left({\sigma }_E^{\partial }-{\sigma }_E^{\partial ,\text{*}}\right)\right\}.$ (40)

Lemma 9. (Discrete continuity of ${\tilde{a}}_h$ .) If $K$ satisfies Hypothesis (${\mathscr{H}}_E^2$ ) for all $E\in ℰ$ and for any vector of positive numbers $ϵ={\left({\epsilon }^E\right)}_{E\in ℰ}$ , there exists a constant $\tilde{C}\left(ϵ\right)>0$ such that

$\forall u_h,v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right),\left|{\tilde{a}}_h\left(u_h,v_h;\overline{u}\right)\right|\le \tilde{C}\left(ϵ\right)\Vert u_h\Vert \Vert v_h\Vert$

where

$\begin{array}{c}\tilde{C}\left(ϵ\right)=\underset{E\in ℰ}{\text{max}}{K}_1^E+\text{max}\left\{\underset{E\in ℰ}{\text{max}}\left(\frac{{\sigma }_E^{in}}{2}\right),\underset{F\in {ℱ}^{\partial }}{\text{max}}\left({\sigma }_E^{\partial }\right)\right\}\\ +\sqrt{2\underset{E\in ℰ}{\text{max}}{\epsilon }^E{K}_1^E\text{max}\left\{\underset{E\in ℰ}{\text{max}}\left(\frac{{\sigma }_E^{in,\text{*}}}{2}\right),\underset{F\in {ℱ}^{\partial }}{\text{max}}\left({\sigma }_E^{\partial ,\text{*}}\right)\right\}}.\end{array}$(41)

Lemma 6 still holds in the two-dimensional case and is left to the reader. Proofs of Lemmas 7, 8 and 9 are similar to proofs in the one-dimensional case. The main difference is in the expression of trace constants. In two dimensions, they are linked to the element’s shape. For an orthonormal basis and for $u\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ , the trace constant of $E\in ℰ$ is given by:

$C_{\text{tr},p}^E=\{\begin{array}{ll}\sqrt{\frac{\left(p+1\right)\left(p+2\right)}{2}},\hfill & \text{if}E\text{is}\text{a}\text{triangle},\hfill \\ \frac{p+1}{2},\hfill & \text{if}E\text{is}\text{a}\text{quadrilateral}.\hfill \end{array}$ (42)

4.4. Optimal Penalization Parameters

Thanks to the previous results on the discrete linearized problem Problem (${\tilde{\mathcal{W}}}_h$ ), one can now construct a method to set automatically penalization parameters. They must be chosen to ensure the coercivity and continuity of the linearized discrete problem, i.e., $C^{\text{*}}\left(ϵ\right)>0$ and $\tilde{C}\left(ϵ\right)>0$ . Moreover, using Céa’s Lemma 10, they are set to minimize the distance between the weak and discrete solutions.

Lemma 10. (Céa’s lemma). Let $V$ be a real Hilbert space with the norm $\Vert \cdot \Vert$ . Let $a:V\times V\to ℝ$ be a bilinear form and $l:V\to ℝ$ a linear form satisfying the Lax-Milgram theorem. Let $V_h$ be a closed subspace of $V$ . Then, there exists a unique $u_h\in V_h$ such that

$\forall v_h\in V_h,a\left(u_h,v_h\right)=l\left(v_h\right)\text{and}\Vert u-u_h\Vert \le \frac{\tilde{C}}{C^{\text{*}}}\Vert u-v\Vert ,\forall v\in V_h$(43)

where $\tilde{C}$ is the continuity constant and $C^{\text{*}}$ the coercivity constant.

Firstly, as a reminder, positivity of continuity and coercivity constants enforce that for all $n\in \left\{0,\cdots ,N-1\right\}$ , ${\epsilon }^{\left(n\right)}<2$ and $\forall n\in \left\{0,\cdots ,N\right\}$ , ${\sigma }_n>{\sigma }_n^{\text{*}}$ . They are given by:

$C^{\text{*}}\left(ϵ\right)=\text{min}\left\{\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{min}}\left(K_0^{\left(n\right)}\left(1-\frac{{\epsilon }^{\left(n\right)}}{2}\right)\right),{\sigma }_0-{\sigma }_0^{\text{*}},{\sigma }_N-{\sigma }_N^{\text{*}},\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{min}}\left(\frac{{\sigma }_n-{\sigma }_n^{\text{*}}}{2}\right)\right\},$ (44)

and

$\begin{array}{c}\tilde{C}\left(ϵ\right)=\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{max}}\left(K_1^{\left(n\right)}\right)+\sqrt{\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{max}}\left(2{\epsilon }^{\left(n\right)}{K}_1^{\left(n\right)}\right)}\sqrt{\text{max}\left({\sigma }_0^{\text{*}},{\sigma }_N^{\text{*}},\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{max}}\left(\frac{{\sigma }_n^{\text{*}}}{2}\right)\right)}\\ +\text{max}\left({\sigma }_0,{\sigma }_N,\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{max}}\left(\frac{{\sigma }_n}{2}\right)\right).\end{array}$ (45)

For the sake of simplicity, let us consider that the variable $\epsilon$ is the same for every element: $\forall n\in \left\{0,\cdots ,N-1\right\}$ , ${\epsilon }^{\left(n\right)}=\epsilon <2$ , and in addition, because penalization parameters are bounded below, let us consider that they are above the lower bound of an amount $\alpha$ constant for every element:

$\forall \alpha >1,\forall n\in \left\{1,\cdots ,N-1\right\},{\sigma }_n=\frac{\alpha }{2\epsilon }{\tilde{\sigma }}_n^{\text{*}},{\sigma }_0=\frac{\alpha }{\epsilon }{\tilde{\sigma }}_0^{\text{*}},{\sigma }_N=\frac{\alpha }{\epsilon }{\tilde{\sigma }}_N^{\text{*}}\text{with}{\tilde{\sigma }}_n^{\text{*}}=\frac{{\left(K_1^{\left(n\right)}{C}_{\text{tr},p-1}^{\left(n\right)}\right)}^2}{K_0^{\left(n\right)}}.$ (46)

Using previous assumptions, it can be noticed that $C^{\text{*}}$ and $\tilde{C}$ are functions of $\epsilon$ and $\alpha$ and can be rewritten:

$C^{\text{*}}\left(\alpha ,\epsilon \right)=\text{min}\left\{K_0\left(1-\frac{\epsilon }{2}\right),\frac{\alpha -1}{\epsilon }{\tilde{\sigma }}_0,\frac{\alpha -1}{\epsilon }{\tilde{\sigma }}_N,\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{min}}\left(\frac{\alpha -1}{\epsilon }\frac{{\tilde{\sigma }}_n}{4}\right)\right\}$(47)

and

$\begin{array}{c}\tilde{C}\left(ϵ\right)=K_1+\sqrt{2\epsilon K_1}\sqrt{\frac{1}{\epsilon }\text{max}\left({\tilde{\sigma }}_0^{\text{*}},{\tilde{\sigma }}_N^{\text{*}},\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{max}}\left(\frac{{\tilde{\sigma }}_n^{\text{*}}}{4}\right)\right)}\\ +\frac{\alpha }{\epsilon }\text{max}\left({\tilde{\sigma }}_0,{\tilde{\sigma }}_N,\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{max}}\left(\frac{{\tilde{\sigma }}_n}{4}\right)\right).\end{array}$(48)

One can see that two quantities are involved in the two previous definitions:

${\tilde{\sigma }}_{\min }=\min \left\{{\tilde{\sigma }}_0,{\tilde{\sigma }}_N,\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\min }\left(\frac{{\tilde{\sigma }}_n}{4}\right)\right\}\text{and}{\tilde{\sigma }}_{\max }=\max \left\{{\tilde{\sigma }}_0,{\tilde{\sigma }}_N,\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\max }\left(\frac{{\tilde{\sigma }}_n}{4}\right)\right\}$ (49)

to have the final write:

$C^{\text{*}}\left(\alpha ,\epsilon \right)=\text{min}\left\{K_0\left(1-\frac{\epsilon }{2}\right),\frac{\alpha -1}{\epsilon }{\tilde{\sigma }}_{\min }\right\}\text{and}\tilde{C}\left(\alpha ,\epsilon \right)=K_1+\sqrt{2K_1{\tilde{\sigma }}_{\max }}+\frac{\alpha }{\epsilon }{\tilde{\sigma }}_{\max }$ (50)

These new expressions of $C^{\text{*}}$ and $\tilde{C}$ show that $C^{\text{*}}$ has two different states and $\tilde{C}$ is continuous concerning $\alpha$ and $\epsilon$ . The aim of this section can now be reformulated as find $\alpha$ and $\epsilon$ such that $\gamma \left(\alpha ,\epsilon \right)=\frac{\tilde{C}\left(\alpha ,\epsilon \right)}{C^{\text{*}}\left(\alpha ,\epsilon \right)}$ is minimal. First, $C^{\text{*}}$ and $\tilde{C}$ are studied separately, then $\gamma$ is observed. $C^{\text{*}}$ has two different states, is continuous and well defined for all $\left(\alpha ,\epsilon \right)\in \left(1,+\infty \right)\times \left(0,2\right)$ . It can be rewritten as follows:

$\forall \left(\alpha ,\epsilon \right)\in \left(1,+\infty \right)\times \left(0,2\right),$ (51)

$C^{\text{*}}\left(\alpha ,\epsilon \right)=\{\begin{array}{ll}\frac{\alpha -1}{\epsilon }{\tilde{\sigma }}_{\min },\hfill & \text{if}\alpha \le {\alpha }^{\text{*}}\left(ϵ\right)\hfill \\ K_0\left(1-\frac{\epsilon }{2}\right),\hfill & \text{otherwise}\hfill \end{array}\text{with}{\alpha }^{\text{*}}\left(ϵ\right)=\frac{K_0}{2{\tilde{\sigma }}_{\min }}\epsilon \left(2-\epsilon \right)+1.$(52)

where $\tilde{C}$ is continuous and well defined for all $\left(\alpha ,\epsilon \right)\in \left(1,+\infty \right)\times \left(0,2\right)$ . $C^{\text{*}}$ and $\tilde{C}$ are now explicitly characterized and now $\gamma \left(\alpha ,\epsilon \right)=\frac{\tilde{C}\left(\alpha ,\epsilon \right)}{C^{\text{*}}\left(\alpha ,\epsilon \right)}$ can be studied. $\left({\alpha }_{opt},{\epsilon }_{opt}\right)$ are looked for such that $\gamma$ is minimal and it is given by:

$\forall \left(\alpha ,\epsilon \right)\in \left(1,+\infty \right)\times \left(0,2\right),$ (53)

$\gamma \left(\alpha ,\epsilon \right)=\{\begin{array}{ll}\frac{\epsilon }{{\tilde{\sigma }}_{\min }\left(\alpha -1\right)}\left(K_1+\sqrt{2K_1{\tilde{\sigma }}_{\max }}\right)+\frac{\alpha }{\alpha -1}\frac{{\tilde{\sigma }}_{\max }}{{\tilde{\sigma }}_{\min }},\hfill & \text{if}\alpha \le {\alpha }^{\text{*}}\left(ϵ\right)\hfill \\ \frac{2}{K_0\left(2-\epsilon \right)}\left(K_1+\sqrt{2K_1{\tilde{\sigma }}_{\max }}+\frac{\alpha }{\epsilon }{\tilde{\sigma }}_{\max }\right),\hfill & \text{otherwise}\hfill \end{array}$(54)

where $\gamma$ is studied on its different open subdomains and the boundary between them. On ${\mathcal{D}}_1$ , for all $\left(\alpha ,\epsilon \right)\in \left({\alpha }^{\text{*}}\left(\epsilon \right),+\infty \right)\times \left(0,2\right)$ , it gives:

$\gamma \left(\alpha ,\epsilon \right)=a\frac{1}{2-\epsilon }+b\frac{\alpha }{\epsilon \left(2-\epsilon \right)}\text{with}a=2\frac{K_1+\sqrt{2K_1{\tilde{\sigma }}_{\max }}}{K_0}\text{and}b=2\frac{{\tilde{\sigma }}_{\max }}{K_0}.$(55)

Then, looking at its variations, it gives that:

${\partial }_{\epsilon }\gamma \left(\alpha ,\epsilon \right)\{\begin{array}{ll}<0,\hfill & \text{if}0<\epsilon <{\epsilon }^{\text{*}}\hfill \\ =0,\hfill & \text{if}\epsilon ={\epsilon }^{\text{*}}\hfill \\ >0,\hfill & \text{if}{\epsilon }^{\text{*}}<\epsilon <2\hfill \end{array}\text{and}{\partial }_{\alpha }\gamma \left(\alpha ,\epsilon \right)=\frac{b}{\epsilon \left(2-\epsilon \right)}>0$ (56)

with ${\epsilon }^{\text{*}}=\frac{\sqrt{b\left(2a+b\right)}-b}{a}>0$ . And finally noting that $\gamma \to +\infty$ when $\alpha \to +\infty$ and when $\epsilon \to 0$ or $\epsilon \to 2$ it gives that $\gamma$ is minimal for $\epsilon ={\epsilon }^{\text{*}}$ and $\alpha \to {\alpha }^{\text{*}}\left({\epsilon }^{\text{*}}\right)$ .

On ${\mathcal{D}}_2$ , for all $\left(\alpha ,\epsilon \right)\in \left(1,{\alpha }^{\text{*}}\left(\epsilon \right)\right)\times \left(0,2\right)$ it gives:

$\gamma \left(\alpha ,\epsilon \right)=\frac{\epsilon }{{\tilde{\sigma }}_{\min }\left(\alpha -1\right)}\left(K_1+\sqrt{2K_1{\tilde{\sigma }}_{\max }}\right)+\frac{\alpha }{\alpha -1}\frac{{\tilde{\sigma }}_{\max }}{{\tilde{\sigma }}_{\min }}.$ (57)

Then, looking at its variations, it gives that:

$\begin{array}{l}{\partial }_{\epsilon }\gamma \left(\alpha ,\epsilon \right)=\frac{K_1+\sqrt{2K_1{\tilde{\sigma }}_{\max }}}{{\tilde{\sigma }}_{\min }\left(\alpha -1\right)}>0\text{and}\\ {\partial }_{\alpha }\gamma \left(\alpha ,\epsilon \right)=-\frac{1}{{\left(\alpha -1\right)}^2}\left(\frac{\epsilon }{{\tilde{\sigma }}_{\min }}\left(K_1+\sqrt{2K_1{\tilde{\sigma }}_{\max }}\right)+\frac{{\tilde{\sigma }}_{\max }}{{\tilde{\sigma }}_{\min }}\right)<0.\end{array}$(58)

And finally, noting that $\gamma \to +\infty$ when $\alpha \to 1$ it gives that $\gamma$ is minimal for $\alpha \to {\alpha }^{\text{*}}\left(\epsilon \right)$ . On the boundary between ${\mathcal{D}}_1$ and ${\mathcal{D}}_2$ , for all $\alpha ={\alpha }^{\text{*}}\left(\epsilon \right)$ and $\epsilon \in \left(0,2\right)$ it gives:

$\begin{array}{l}\gamma \left({\alpha }^{\text{*}}\left(\epsilon \right),\epsilon \right)=a\frac{1}{2-\epsilon }+b\frac{1}{\epsilon \left(2-\epsilon \right)}+c\text{with}a=2\frac{K_1+\sqrt{2K_1{\tilde{\sigma }}_{\max }}}{K_0},\\ b=2\frac{{\tilde{\sigma }}_{\max }}{K_0}\text{and}c=\frac{{\tilde{\sigma }}_{\max }}{{\tilde{\sigma }}_{\min }}.\end{array}$(59)

Then, looking at its variations, it gives that:

${\partial }_{\epsilon }\gamma \left({\alpha }^{\text{*}}\left(\epsilon \right),\epsilon \right)\{\begin{array}{ll}<0,\hfill & \text{if}0<\epsilon <{\epsilon }_{opt}\hfill \\ =0,\hfill & \text{if}\epsilon ={\epsilon }_{opt}\hfill \\ >0,\hfill & \text{if}{\epsilon }_{opt}<\epsilon <2\hfill \end{array}\text{with}{\epsilon }_{opt}=\frac{\sqrt{b\left(2a+b\right)}-b}{a}>0.$ (60)

The expression of $\left({\alpha }_{opt},{\epsilon }_{opt}\right)$ can be summarized as follows:

$\begin{array}{l}{\epsilon }_{opt}=\frac{\sqrt{b\left(2a+b\right)}-b}{a}\text{with}a=2\frac{K_1+\sqrt{2K_1{\tilde{\sigma }}_{\max }}}{K_0}\text{and}b=2\frac{{\tilde{\sigma }}_{\max }}{K_0}\text{and}\\ {\alpha }_{opt}=\frac{K_0}{2{\tilde{\sigma }}_{\min }}{\epsilon }_{opt}\left(2-{\epsilon }_{opt}\right)+1.\end{array}$(61)

Finally, in one dimension, the auto-calibration of penalization parameters is given by:

$\begin{array}{l}\forall n\in \left\{1,\cdots ,N-1\right\},{\sigma }_n=\frac{{\alpha }_{opt}}{2{\epsilon }_{opt}}{\tilde{\sigma }}_n^{\text{*}},{\sigma }_0=\frac{{\alpha }_{opt}}{{\epsilon }_{opt}}{\tilde{\sigma }}_0^{\text{*}},{\sigma }_N=\frac{{\alpha }_{opt}}{{\epsilon }_{opt}}{\tilde{\sigma }}_N^{\text{*}}\\ \text{with}{\tilde{\sigma }}_n^{\text{*}}=\frac{{\left(K_1^{\left(n\right)}{C}_{\text{tr},p-1}^{\left(n\right)}\right)}^2}{K_0^{\left(n\right)}}.\end{array}$(62)

In two dimensions, the auto-calibration of penalization parameters is given by:

$\{\begin{array}{l}\forall F\in {ℱ}^{\text{in}},{\sigma }_E^{in}=\frac{{\alpha }_{opt}}{2{\epsilon }_{opt}}{\sigma }_E^{\text{*}},{\sigma }_{E_r}^{in}=\frac{{\alpha }_{opt}}{2{\epsilon }_{opt}}{\sigma }_{E_r}^{\text{*}}\\ \forall F\in {ℱ}^{\partial },{\sigma }_E^{\partial }=\frac{{\alpha }_{opt}}{{\epsilon }_{opt}}{\sigma }_E^{\text{*}}\end{array}\text{with}{\sigma }_E^{\text{*}}=\frac{D^E{\left(K_1^E{C}_{\text{tr},p-1}^E\right)}^2}{2{\epsilon }^E{K}_0^E}$ (63)

and $D^E$ is the number of edges of the element $E$ and $C_{\text{tr},p-1}^E$ is the trace constant defined in Equation (42).

Extension to multi-dimensional meshes. In multi-dimensional settings, the penalization term given in Equation (62) is extended by replacing the one-dimensional element length $d_E$ with a local characteristic size derived from the ratio between the element volume and its boundary surface area, i.e., $d_E=\left|E\right|/\left|\partial E\right|$ . For non-quadrilateral (triangular or polygonal) cells, this characteristic measure ensures that the penalty parameter scales consistently with the element geometry. Therefore, no additional geometric constant is required, and the same penalization formula applies naturally in two and three dimensions.

4.5. Convergence of the Discrete Linearized Weak Problem to the Continuous Linearized Weak Problem

Previously, it has been proven that the Problem (${\tilde{\mathcal{W}}}_h$ ) has a unique solution. This problem is part of a fixed-point method, and it has been proven in Section 4.2 that this fixed point has a unique solution also. To solve the nonlinear weak formulation Problem ($\mathcal{W}$ ), one step needs to be added to prove the well-posedness of the problem. It is addressed in the following; the goal is to prove that the solution of Problem (${\tilde{\mathcal{W}}}_h$ ) converges towards the solution of Problem ($\tilde{\mathcal{W}}$ ) and prove that the bilinear form ${\tilde{a}}_h$ of Problem (${\tilde{\mathcal{W}}}_h$ ) converges to Problem ($\mathcal{W}$ ).

The work in this section is based on the book of Pietro and Ern published in 2012 [13]. They proved convergence in the case of a Symmetric Interior Penalty Galerkin method and sketch the proof in the case of an Incomplete Interior Penalty method. The following study provides detailed proof of the IIPG case.

The key idea is to revisit the concept of consistency and introduce a new point of view based on asymptotic consistency. This new form of consistency and the usual stability of the discrete bilinear form are the two main ingredients for asserting convergence to the minimal regularity solutions. The discrete bilinear form $a_h$ needs to be reformulated to consider only the contribution of $K$ on the mesh elements, not the interfaces; consequently, lifting operators are introduced. They map functions defined on mesh faces to functions defined on mesh elements. In the context of DG methods, liftings act on interfaces and boundary jumps. Bassi and Rebay introduced them [36] in the context of compressible flows and analyzed by Brezzi et al. [37] in the context of the Poisson problem. Liftings have many useful applications. They can be combined with the gradient to define discrete gradients. Discrete gradients play an essential role in the design and analysis of DG methods. Indeed, they can be used to formulate the discrete problem locally on each element using numerical fluxes.

Liftings: Definition

For any point $x_n$ , and for all $\phi \in L^2\left(\left\{x_n\right\}\right)$ the lifting operator $r_n^p:L^2\left(\left\{x_n\right\}\right)\to {\mathcal{V}}_0^p\left({ℰ}_h\right)$ is defined as the solution of the following problem:

${\displaystyle {\int }_{\Omega }{r}_n^p\left(\phi \right){\tau }_h\text{d}x}={\left\{\left|{\tau }_h\right|\right\}}_{x_n}\phi \left(x_n\right),\forall {\tau }_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right).$ (64)

For any $v$ in ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ , the global lifting of its interface and boundary jumps is defined as follows:

$R_h^p\left(〚v〛\right):={\displaystyle \sum }_{n=0}^N{r}_n^p\left(〚v〛\right)\in {\mathcal{V}}_0^p\left({ℰ}_h\right).$ (65)

Discrete gradients: Definition

The discrete gradient operator $G_h^p:{\mathcal{V}}_0^p\left({ℰ}_h\right)\to L^2\left(I_n\right)$ is defined as follows: for all $v$ in ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ ,

$G_h^p\left(v\right):={\nabla }_{h}v-R_h^p\left(〚v〛\right).$ (66)

In addition, there exists a bound on the discrete gradient operator:

${\Vert G_h^p\left(v\right)\Vert }_{L^{\text{2}}\left(\Omega \right)}\le \alpha \Vert v\Vert$ (67)

where $\Vert \cdot \Vert$ is the norm associated with the IIPG formulation defined Equation (30).

Theorem 11. (Regularity of the limit and weak asymptotic consistency of discrete gradients.) Let $p\ge 0$ . Let $v_h$ be a sequence in ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ bounded by the $\Vert .\Vert$ -norm. Then, there is a function $v\in H_0^1\left(\Omega \right)$ such that as $h\to 0$ , up to a subsequence,

$v_h\to v\text{strongly}\text{in}{L}^{\text{2}}\left(\Omega \right),$ (68)

and for all $p\ge 0$ , the discrete gradients defined by Equation (66) are such that

$G_h^p\left(v_h\right)\rightharpoonup v^{\prime }\text{weakly}\text{in}{L}^{\text{2}}\left(\Omega \right).$ (69)

Proof of Theorem 11 is available in [13] (pp. 194-195).

Because of the shape of the IIPG formulation, the modified discrete gradient operator ${\widehat{G}}_h^p:{\mathcal{V}}_0^p\left({ℰ}_h\right)\to L^2\left(I_n\right)$ is defined as follows: for all $v$ in ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ ,

${\widehat{G}}_h^p\left(v\right):={\nabla }_{h}v.$ (70)

Using liftings and discrete gradients, surface contributions of the flux in Equation (32) are transformed to volume contribution. It makes working with the bilinear form ${\tilde{a}}_h$ easier. For a given $\overline{u}\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ , it can be rewritten as follows:

$\begin{array}{c}\forall u_h,v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right),\\ {\tilde{a}}_h\left(u_h,v_h;\overline{u}\right)={\displaystyle \sum }_{n=0}^{N-1}{\displaystyle {\int }_{I_n}K}\left(x,\overline{u}\right){\nabla }_h{u}_h{\nabla }_h{v}_h\text{d}x-{\displaystyle \sum }_{n=0}^N{\left\{\left|K\left(x,\overline{u}\right){\nabla }_h{u}_h\right|\right\}}_{x_n}{〚v_h〛}_{x_n}+s_h\left(u_h,v_h\right)\\ ={\displaystyle \sum }_{n=0}^{N-1}{\displaystyle {\int }_{I_n}K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right){\nabla }_h{v}_h\text{d}x-{\displaystyle \sum }_{n=0}^N{\displaystyle \sum }_{m=0}^{N-1}{\displaystyle {\int }_{I_m}K}\left(x,\overline{u}\right)r_n^p\left(〚v_h〛\right){\widehat{G}}_h^p\left(u_h\right)+s_h\left(u_h,v_h\right)\\ ={\displaystyle \sum }_{n=0}^{N-1}{\displaystyle {\int }_{I_n}K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right){\nabla }_h{v}_h\text{d}x-{\displaystyle \sum }_{n=0}^{N-1}{\displaystyle {\int }_{I_n}K}\left(x,\overline{u}\right)R_h^p\left(〚v_h〛\right){\widehat{G}}_h^p\left(u_h\right)+s_h\left(u_h,v_h\right)\\ ={\displaystyle \sum }_{n=0}^{N-1}{\displaystyle {\int }_{I_n}K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right)G_h^p\left(v_h\right)\text{d}x+s_h\left(u_h,v_h\right)\end{array}$ (71)

with

$\begin{array}{l}\forall u_h,v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right),s_h\left(u_h,v_h\right)=\frac{{\sigma }_0}{h}{〚u_h〛}_{x_0}{〚v_h〛}_{x_0}+{\displaystyle \sum }_{n=1}^{N-1}\frac{{\sigma }_{n-1}+{\sigma }_n}{2h}{〚u_h〛}_{x_n}{〚v_h〛}_{x_n}\\ +\frac{{\sigma }_N}{h}{〚u_h〛}_{x_N}{〚v_h〛}_{x_N}.\end{array}$ (72)

Consider that ${\left({\sigma }_n\right)}_{n=0,\cdots ,N}$ are chosen according to Lemma 4 that implies discrete coercivity in the $\Vert .\Vert$ -norm, and hence well-posedness of the discrete linearized problem $\left({\tilde{\mathcal{V}}}_h\right)$ .

Definition 1. (Asymptotic adjoint consistency.) The discrete bilinear form $a_h$ is asymptotically adjoint consistent with the exact bilinear form $a$ on ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ if for any subsequence $v_h$ in ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ bounded in the $\Vert .\Vert$ -norm and for any smooth function $\phi \in C_0^{\infty }\left(\Omega \right)$ , there is a subsequence ${\phi }_h$ in ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ converging to $\phi$ in the $\Vert .\Vert$ -norm and such that, up to a subsequence

$\underset{h\to 0}{\lim }{a}_h\left(v_h,{\phi }_h\right)=a\left(v,\phi \right)={\displaystyle {\int }_{\Omega }{v}^{\prime }{\phi }^{\prime }\text{d}x}$ (73)

where $v\in H_0^1\left(\Omega \right)$ is the limit of the subsequence identified in Theorem 11.

Lemma 12. (Asymptotic adjoint consistency of ${\tilde{a}}_h$ .) The discrete bilinear form ${\tilde{a}}_h$ of Problem (${\tilde{\mathcal{W}}}_h$ ) is asymptotically adjoint consistent with the exact bilinear form $\tilde{a}$ of Problem ($\tilde{\mathcal{W}}$ ) on ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ .

Finally, we deduce the following result.

Theorem 13. (Convergence to minimal regularity solutions.) Let $p\ge 1$ . Let $u_h$ be a sequence of approximate solutions generated by solving the discrete linearized problem $\left({\tilde{\mathcal{V}}}_h\right)$ with ${\tilde{a}}_h$ defined by Equation (32) and with penalty parameters ensuring coercivity. Then, as $h\to 0$

$u_h\to u\text{strongly}\text{in}{L}^{\text{2}}\left(\Omega \right)$ (74)

${\nabla }_h{u}_h\to u^{\prime }\text{strongly}\text{in}{L}^{\text{2}}\left(\Omega \right)$ (75)

${\left|u_h\right|}_J\to 0$ (76)

where $u\in H_0^1\left(\Omega \right)$ is the unique solution of the strong problem.

Proofs of Lemma 12 and Theorem 13 can be found in Appendix.

4.6. Concluding Results

In the current section, several theorems have been proven. It is proven that there exists a unique solution to Problem ($\mathcal{W}$ ) using Lemma 1 and Lemma 2. Then, it is proven that for a given $\overline{u}$ , there exists a unique solution to Problem (${\tilde{\mathcal{W}}}_h$ ) using Lemma 3. Lastly it is proven that for a given $\overline{u}$ , the solution of Problem (${\tilde{\mathcal{W}}}_h$ ) converges to the solution of Problem ($\tilde{\mathcal{W}}$ ). These results proven in a general case for a given $\overline{u}$ can be used to solve the toy problem. Figure 6 gives a graphical representation of the whole loop of resolution with different paths.

Figure 6. Scheme of the whole loop of resolution with the different linearization methods.

The nonlinear problem, Problem ($\mathcal{W}$ ) can be linearized directly at the continuous level by employing a fixed-point method. The continuous level linearization

$\begin{array}{l}T:H_0^1\left(\Omega \right)\to H_0^1\left(\Omega \right)\\ \overline{u}\mapsto T\left(\overline{u}\right)=u\end{array}$(77)

stands for: find $u$ solution of Problem ($\tilde{\mathcal{W}}$ ) for a given $\overline{u}\in H_0^1\left(\Omega \right)$ . One can define the following sequence defined by $u^0\in H_0^1\left(\Omega \right)$ an initial guess and $u^{n+1}=T\left(u^n\right)$ for $n\in \mathbb{N}$ . Lemma 1 and Lemma 2 ensure that taking ${\lim }_{n\to \infty }{u}^n$ gives the solution of Problem ($\mathcal{W}$ ).

A discretization step is needed to compute the solution of Problem ($\mathcal{W}$ ). Consequently, the projector $P_h:H_0^1\left(\Omega \right)\to {\mathcal{V}}_0^p\left({ℰ}_h\right)$ is introduced. It projects a function living in an infinite-dimensional space to a finite-dimensional space, especially it projects a function to the DG space $V_0^p\left({ℰ}_h\right)$ . Then, at a discrete level, the linearization method

$\begin{array}{l}{T}_h:V_0^p\left({ℰ}_h\right)\to V_0^p\left({ℰ}_h\right)\\ \overline{u}\mapsto T_h\left(\overline{u}\right)=u_h\end{array}$(78)

stands for: find $u_h$ discrete solution of Problem (${\tilde{\mathcal{W}}}_h$ ) for a given $\overline{u}\in V_0^p\left({ℰ}_h\right)$ . One can notice that for a given $\overline{u}\in H_0^1\left(\Omega \right)$ , it has been proven (Theorem 13) that $\left(T_h\circ P_h\right)\left(\overline{u}\right)=u_h$ converges to $u$ given by $T_C\left(\overline{u}\right)$ .

Lastly, the linearization method of Problem ($\mathcal{W}$ ) going through a discretization step is defined as

$\begin{array}{l}{T}_D:H_0^1\left(\Omega \right)\to H_0^1\left(\Omega \right)\\ \overline{u}\mapsto T_D\left(\overline{u}\right)=\underset{h\to 0}{\lim }\left(T_h\circ P_h\right)\left(\overline{u}\right)=u\end{array}$(79)

Using $T_D$ one can define a new sequence $v^0\in H_0^1\left(\Omega \right)$ an initial guess and $v^{n+1}=T_D\left(v^n\right)={\lim }_{h\to 0}\left(T_h\circ P_h\right)\left(v^n\right)$ for $n\in \mathbb{N}$ . Taking the limit when $n$ goes to infinity gives the solution of Problem ($\mathcal{W}$ ).

The previously explained method uses two limits, $h$ goes to 0 then $n$ goes to infinity. One can also consider limits in the opposite order. Using proof of Lemma 1 applied to the nonlinear discrete problem and then using Theorem 13, one can prove that the solution of the nonlinear discrete problem converges to the solution of the nonlinear continuous problem.

5. Numerical Results

Following the numerical methods and theoretical results presented in the previous sections, the RIVAGE code is validated against numerical test cases. Two analytical test cases are used to compute convergence rates and validate the code. These analytical test cases are obtained by considering the problem’s aimed solution and choosing the source term according to the solution and the hydraulic conductivity function. They are built upon the nonlinear Poisson’s equation. The first case is a nonlinear one-dimensional problem in its stationary form. The second case is a nonlinear two-dimensional problem in its stationary form. These numerical experiments are inspired by literature. In 2008, Rivière [14] and in 2021, Clément et al. [12] computed convergence rates for linear problems, also for nonlinear problems.

Stationary problems are considered since theoretical results are given on this type of problem. Moreover, they are more difficult to solve since they solve the problem at infinite time. Consequently, the nonlinear solver has to find the solution without getting time sub-steps.

Experimental test cases are solved with the RIVAGE code. These problems aim at confirming the performance of the adaptive strategy proposed in this work. Moreover, they allow to test RIVAGE of problems encountered in the hydrology field. These experiments are based on the work of Haverkamp et al. [38] and Vauclin et al. [39].

Table 3. Solver and time-integration settings used in all numerical experiments.

Parameter	Symbol/Setting	Description
Nonlinear tolerance	${\epsilon }_{\text{NL}}={10}^{-6}$	Residual tolerance for Picard iterations
Maximum Picard iterations	$N_{\text{Pic}}^{\text{max}}=40$	Upper bound before step rejection
Time-step size	adaptive	Automatically reduced if residual increases
Minimum time-step	$\text{Δ}{t}_{\text{min}}={10}^{-4}\text{s}$	Stability safeguard
Linear solver	Conjugate Gradient (CG)	Preconditioned by ILU(0)
Spatial polynomial degree	$p=1,2$	Depending on the test case
Penalty parameters	Equation (63)	Auto-calibrated during iteration

These settings (see Table 3) were kept identical for all test cases unless explicitly stated otherwise.

5.1. One-Dimensional Analytical Test Case

For this first test case, theoretical convergence rates of the IIPG methods are checked, and numerical stability is evaluated concerning penalty values and penalization methods. The following problem is considered:

$\{\begin{array}{l}-\left(K\left(u\right)u\right)=f\left(x\right)\text{in}\Omega =\left[-1,1\right]\\ u\left(-1\right)=1,\\ u\left(1\right)=-1,\end{array}$ (80)

with $K\left(u\right)=\text{tanh}\left(5u\right)+1.01$ and $f$ obtained by replacing $u$ by $u_{ex}$ in the problem. The chosen analytical solution is $u_{ex}\left(x\right)=-\text{sin}\left(\frac{\pi }{2}x\right)$ . The analytical solution is chosen not to be polynomial but to span the interval $\left[-1,1\right]$ . The hydraulic conductivity is chosen to have a nonlinear problem with a similar shape of law given in Table 1. tanh has been chosen because it is a smooth function convenient for the computation of convergence rates and looks like constitutive laws for RE. Moreover, a factor of 200 between the maximum and the minimum value of $K$ with $K_0=0.01$ . The problem is solved with the IIPG method. Three types of penalization are used. The first one ${\sigma }^E=\sigma =1$ for all $E\in ℰ$ , the second one ${\sigma }^E=\sigma =100$ for all $E\in ℰ$ and the third one ${\sigma }^E$ are auto-calibrated using the method presented in Section 4. For each type of penalization, the solution is approximated by a piecewise linear function ($p=1$ ), a piecewise quadratic function ($p=2$ ), and a piecewise cubic function ($p=3$ ). Moreover, lastly, four different mesh sizes are used $N_x=20,40,80,160$ with $N_x$ —the number of elements in the equally spaced partition of Ω.

Table 4. $L^2$ -error, convergence rates and number of iterations for the one-dimensional benchmark.

		$p=1$			$p=2$			$p=3$
$\sigma$	$N_x$	$L^2$ -error	$r$	$t\left(s\right)$	$L^2$ -error	$r$	$t\left(s\right)$	$L^2$ -error	$r$	$t\left(s\right)$
1	20	3.21 × 10−1		0.21	1.33 × 10−1		1.94	2.29 × 10−4		6.21
-	40	1.29 × 10−1	1.31	0.46	3.41 × 10−2	1.97	3.17	1.42 × 10−5	4.01	11.85
-	80	3.77 × 10−2	1.78	1.02	8.53 × 10−3	2.00	5.97	8.88 × 10−7	4.00	23.94
-	160	9.83 × 10−3	1.94	2.08	2.13 × 10−3	2.00	12.09	5.62 × 10−8	3.98	52.83
-	Fitted		1.69			1.99			4.00
100	20	8.33 × 10−3		0.21	1.36 × 10−3		1.48	2.33 × 10−6		5.89
-	40	2.10 × 10−3	1.99	0.51	3.41 × 10−4	2.00	2.97	1.44 × 10−7	4.02	11.87
-	80	5.27 × 10−4	2.00	1.03	8.53 × 10−5	2.00	5.94	9.74 × 10−9	3.89	24.03
-	160	1.31 × 10−4	2.00	2.08	2.13 × 10−5	2.00	12.16	1.37 × 10−9	2.83	53.20
-	Fitted		1.99			2.00			3.61
auto	20	3.53 × 10−2		0.24	1.69 × 10−2		1.43	3.46 × 10−6		5.88
-	40	8.88 × 10−3	1.99	0.51	4.40 × 10−3	1.94	2.86	1.61 × 10−7	4.42	11.92
-	80	2.17 × 10−3	2.03	1.05	1.13 × 10−3	1.95	5.95	9.45 × 10−9	4.10	24.07
-	160	5.32 × 10−4	2.03	2.55	2.90 × 10−4	1.97	12.13	1.33 × 10−9	2.82	53.25
-	Fitted		2.02			1.96			3.81

Figure 7. Penalization parameters for the one-dimensional test case in the case of auto-penalization.

Table 4 shows $L^2$ -error and convergence rate for each computation. It can be noticed that computed convergence rates correspond to the theoretical ones found in literature [14] and [40] (pp. 64-84). For the IIPG formulation with penalization, $p$ is odd in order $p+1$ , it is optimal, and if $p$ is even, the order is $p$ , it is suboptimal. Moreover, for a penalization speed set by the user to 1 (outside of the range specified by theoretical results), errors are about 100 times greater than in other computations. The fixed-point method converges to a less accurate solution. Computation times are also given. It can be noticed that auto-penalization is not greatly slower than user-defined penalization and can even be faster due to the quickest convergence of the iterative method.

Moreover, Figure 7 shows penalization values in the case of auto-calibration. One can observe that penalization values are not constant on Ω and vary according to the polynomial degree of approximation. On the domain, some part needs a small amount of penalization, whereas others need a higher amount.

5.2. Two-Dimensional Analytical Test Case

This second experiment focuses on the ability of the IIPG method to solve RE in two dimensions. Its convergence rates are computed, and numerical stability is evaluated concerning penalty values and penalization methods. The following problem is considered:

$\{\begin{array}{ll}-\nabla \cdot \left(K\left(u\right)\nabla u\right)=f\left(x\right)\hfill & \text{in}\Omega =\left[-1,1\right]\times \left[-1,1\right]\hfill \\ u=0\hfill & \text{on}\partial \Omega ,\hfill \end{array}$ (81)

with $K\left(u\right)=\text{tanh}\left(u\right)+1.01$ and similarly to the previous test case $f$ is obtained by replacing $u$ by $u_{ex}$ in the problem. The chosen analytical solution is $u_{ex}\left(x,y\right)=\sin \left(\frac{\pi }{2}x\right)\sin \left(\frac{\pi }{2}y\right)$ . The problem is solved similarly to the one-dimensional test case. Three types of penalization are used. The first one ${\sigma }^E=\sigma =1$ for all $E\in ℰ$ , the second one ${\sigma }^E=\sigma =100$ for all $E\in ℰ$ and the third one ${\sigma }^E$ are auto-calibrated. For each type of penalization, the solution is approximated by a piecewise linear function ($p=1$ ), a piecewise quadratic function ($p=2$ ), and a piecewise cubic function ($p=3$ ). Lastly, three different meshes are used. They are all composed of quadrilaterals of identical size, and each space direction is discretized with $N=10,20,40$ elements. It gives a mesh with $N_E=100,400,1600$ elements.

Table 5. $L^2$ -error, convergence rates and number of iterations for the two-dimensional benchmark.

		$p=1$			$p=2$			$p=3$
$\sigma$	$N_x$	$L^2$ -error	$r$	$t\left(s\right)$	$L^2$ -error	$r$	$t\left(s\right)$	$L^2$ -error	$r$	$t\left(s\right)$
1	10	6.45 × 10−2		0.29	4.83 × 10−2		1.54	8.60 × 10−4		5.07
-	20	1.51 × 10−2	1.99	1.00	1.11 × 10−2	2.11	7.21	4.69 × 10−5	4.20	27.21
-	40	3.53 × 10−3	2.10	7.57	2.65 × 10−3	2.07	59.51	2.74 × 10−6	4.09	279.59
-	Fitted		2.10			2.09			4.15
100	10	3.80 × 10−2		0.25	2.02 × 10−3		1.14	7.32 × 10−5		4.83
-	20	9.53 × 10−3	1.99	0.99	2.72 × 10−4	2.90	6.78	4.59 × 10−6	4.00	30.23
-	40	2.38 × 10−3	2.00	8.37	4.08 × 10−5	2.74	61.62	2.87 × 10−7	4.00	290.86
-	Fitted		2.00			2.82			4.00
auto	10	3.37 × 10−2		0.25	2.52 × 10−3		1.15	7.41 × 10−5		4.93
-	20	8.11 × 10−3	2.06	1.03	5.90 × 10−4	2.09	6.88	4.71 × 10−6	3.98	30.15
-	40	2.02 × 10−3	2.00	8.49	1.51 × 10−4	1.96	60.85	2.97 × 10−7	3.99	288.50
-	Fitted		2.03			2.03			3.98

Table 5 shows $L^2$ -error and convergence rate for each computation. It can be noticed that computed convergence rates correspond to the theoretical ones found in literature [14] and [40] (pp. 64-84). For the IIPG formulation with penalization, $p$ is odd in order $p+1$ , it is optimal, and if $p$ is even, the order is $p$ and suboptimal. Moreover, for a penalization speed set by the user to 1 (outside of the range specified by theoretical results), errors are about 100 times greater than in other computations. The fixed-point method converges to a less accurate solution. Computation times are also given. It can be noticed that auto-penalization is not greatly slower than user-defined penalization and can even be faster due to the quickest convergence of the iterative method as in the one-dimensional case.

Moreover, Figure 8 shows penalization values in the case of auto-calibration. One can observe that penalization values are not constant on Ω and vary according to the polynomial degree of approximation. On the domain, some part needs a small amount of penalization, whereas others need a higher amount.

Figure 8. Penalization parameters for the two-dimensional test case in the case of auto-penalization.

5.3. Application to Groundwater Flows I: Haverkamp’s Test Case

The two problems considered here, one-dimensional and two-dimensional, aim to validate the numerical resolution of RE using DG methods and auto-calibration of penalization parameters. Numerical results are compared to numerical simulations in the literature and experimental data.

The first experimental validation of solving RE with DG methods is a one-dimensional test case. The numerical results are compared with data sourced from the literature. This particular numerical test case was initially presented by Celia et al. [41]. It is based on an experiment conducted by Haverkamp et al. [38], who referred to the availability of a quasi-analytical solution provided by Philip [42]. Subsequently, it was used by others such as [43] [44], and represents a set of well-established test cases, for instance, see [30]. Despite its simplicity, this case offers insights into the fundamental physics of a wetting front resulting from infiltration.

This scenario involves the one-dimensional infiltration into a soil column measuring 40 cm in height and 8 cm in width. The hydraulic head at the top and bottom is governed by Dirichlet boundary conditions: $h_{top}=19.3\text{cm}$ and $h_{bottom}=-61.5\text{cm}$ , resulting in cumulative downward infiltration. The sides are impermeable. The initial condition is $h_0=-61.5+z\text{cm}$ . Although this case is one-dimensional, it is solved on a two-dimensional domain. Therefore, homogeneous Neumann boundary conditions are applied along the boundary in the infiltration direction. For a visual representation of this setup, refer to Figure 9.

Figure 9. Haverkamp’s test case configuration.

Hydraulic properties use Vachaud’s relations in Table 1 with $A=1.175\times {10}^6$ , $B=4.74$ , $C=1.611\times {10}^6$ , $D=3.96$ , $K_s=0.0094\text{cm}\cdot {\text{s}}^{-1}$ , ${\theta }_s=0.287$ , and ${\theta }_r=0.075$ . The simulation is done on a mesh of 160 elements along the $z$ -axis. The solution is piecewise linear ($p=1$ ), and time integration is BDF of order 2. Penalization parameters are set automatically using results from Section 4. In addition, stopping criteria are set to 10⁻⁶ for this computation. The solution to this problem is computed at $T=600\text{s}$ .

Figure 10 displays the comparison of numerical results with results from Manzini et al. [44], the pressure head distribution at $t=360\text{s}$ and the penalization parameters distribution at $t=360\text{s}$ . Numerical results are in good agreement with the literature results for this test case. The pressure head distribution shows a vertical progression of the wetting front with a steep transition from the initial $\psi$ to $\psi$ imposed at the boundary condition. Moreover, the distribution of penalization parameters shows that the penalization parameters are not constant on the whole domain and are higher on the wetting front.

This test case validates a real, evolving test case for the DG method. Moreover, it gives a good insight into the behavior of automatic penalization. Penalization parameters are auto-calibrated as long as the solution evolves. Moreover, automatic penalization impacts a full nonlinear problem because the nonlinear solver needs fewer iterations to converge to the solution.

Figure 10. Haverkamp’s test case, numerical solution for 160 elements, $p=1$ and BDF-2 method.

5.4. Application to Groundwater Flows II: Vauclin’s Test Case

Vauclin, Vachaud, and Khanji conducted a series of laboratory experiments in the 1970s, the details of which can be found in [39]. These experiments explored water table recharge and drainage in a slab of sandy soil. The work by Vauclin et al. [39] specifically focuses on simulating water flow recharge through a soil slab and provides experimental details and results. The experiment involved a 6 m by 2 m box, with only one half simulated due to symmetry. The left, top (for $x>50\text{cm}$ ), and bottom sides were impervious, with a prescribed constant flux on the top for $x\le 50\text{cm}$ of $u_g\cdot n=-14.8\text{cm}\cdot {\text{h}}^{-1}$ . The water level was maintained at a constant $h=65\text{cm}$ in the ditch on the right for $z\le 65\text{cm}$ , while the remaining boundary on the right for $z>65\text{cm}$ accounted for a seepage boundary condition. The initial state was at hydrostatic equilibrium with the water table at $z=65\text{m}$ . For further reference, please see Figure 11 for a schematic representation of the setup. The complete simulation of water table recharge by Vauclin et al. [39] has been used by numerous studies to evaluate their methods (see, for instance, [45]-[47]). The MODFLOW code validation partially relies on this experimental dataset [31].

Hydraulic properties use Vachaud’s relations in Table 1 with $A=2.99\times {10}^6$ , $B=5.0$ , $C=40000$ , $D=2.9$ , $K_s=35\text{cm}\cdot {\text{h}}^{-1}$ , ${\theta }_s=0.3$ , and ${\theta }_r=0.0$ . The simulation is carried on an evolving mesh. The mesh is adapted along the computation according to the gradient of $h$ . Mesh adaptive parameters are set to ${\beta }_c=50$ and ${\beta }_r=50$ . The solution is sought piecewise linear ($p=2$ ) and time integration is BDF of order 3. Penalization parameters are set automatically using results from Section 4. In addition, stopping criteria are set to 10⁻⁶ for this computation. The solution of this problem is computed until $T=10\text{h}$ .

Figure 11. Vauclin’s test case configuration.

In the initial mesh displayed in Figure 12, the refinement below the water entry edge aims to assist in simulating the steep wetting front. Figure 13 compares the water table’s position at $t=2,3,4,8\text{h}$ with data from Vauclin et al. [39]. The numerical results closely match the experimental profile, although there are small discrepancies in the middle of the water table, which may be due to the non-perfect isotropic and homogeneous nature of the sandy soil.

Figure 14 and Figure 15 illustrate the field distribution of hydraulic head, flux, and the positions of the water table and capillary fringe at $\theta =0.29$ . These figures also show the isolines of the hydraulic head. The numerical results are in agreement with the data from Vauclin et al. [39].

Additionally, in Figure 16, the evolution of penalization parameters during the computation is presented. At selected times, the evolution of the mesh reflects the capture of the steep front.

Finally, Figure 17 displays the evolution of time-steps and the number of elements over time. The adaptation of time steps and the number of elements is evident, with the time steps initially small due to the strong nonlinearity induced by the steep wetting front. As the front smoothens, the number of elements decreases, stabilizing at $N_{ele}=600$ after $t=3\text{h}$ .

This test case is a test case, which is a typical problem where auto-calibration of penalization parameters is essential. Since the problem is strongly nonlinear and evolving, with a basic penalization and user defined parameters, the nonlinear solver failed to capture the solution or necessitates some combination of fixed-point solver and Newton-Raphson method such as in the work of [7].

Discussions on possible limitations. A potential limitation of the proposed approach arises when hydraulic parameters (e.g., $K_s$ , $\alpha$ , $n$ ) exhibit abrupt spatial variations between adjacent elements. In such cases, the optimal scaling of the penalty parameter may deteriorate, leading to sub-optimal convergence rates. Similarly, the application of capillary-pressure regularization can alter the local nonlinearity of the constitutive laws, thereby reducing the effectiveness of the automatic penalization scheme. Future work will investigate adaptive penalty strategies that account for local parameter contrasts and regularized constitutive models.

Figure 12. Vauclin’s test case, initial mesh.

Figure 13. Vauclin’s test case, numerical water table position compared to experimental data from Vaulin et al. [39].

Figure 14. Vauclin’s test case, at $t=3\text{h}$ , spatial distribution of hydraulic head, water table position (white line), contour plot of hydraulic head (red lines) and flux (arrows).

Figure 15. Vauclin’s test case, at $t=8\text{h}$ , spatial distribution of hydraulic head, water table position (white line), contour plot of hydraulic head (red lines) and flux (arrows).

Figure 16. Vauclin’s test case, spatial distribution of penalization parameters and mesh at selected times.

Figure 17. Vauclin’s test case, evolution along time of time-steps (left) and number of elements (right).

Acknowledgements

This work has been supported by the ADEN-MED project (Adaptability to Extreme events and Natural risks-application to the Mediterranean and Djibouti), funded by the Région Sud Provence-Alpes-Côte d’Azur under the AAP MEDCLIMAT “Natural risks and food sovereignty”, and by France 2030 through the Priority Research Program and Equipment (PEPR) “Maths-Vives-Mathematics in Interactions”, targeted project HYDRAUMATH (ANR-23-EXMA-007), operated by ANR.

Appendix

Proofs on Theoretical Results

Proof of Lemma 4. For a given $\overline{u}\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ and choosing $v_h=u_h$ in (32) yields

$\begin{array}{l}\forall u_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right),{\tilde{a}}_h\left(u_h,u_h\right)={\displaystyle \sum }_{n=0}^{N-1}{\displaystyle {\int }_{I_n}K}\left(x,\overline{u}\right){\left({u^{\prime }}_h\right)}^2\text{d}x-{\displaystyle \sum }_{n=0}^N{\left\{\left|K\left(x,\overline{u}\right){u^{\prime }}_h\right|\right\}}_{x_n}{〚u_h〛}_{x_n}\\ +\frac{{\sigma }_0}{h}{〚u_h〛}_{x_0}^2+{\displaystyle \sum }_{n=1}^{N-1}\frac{{\sigma }_{n-1}+{\sigma }_n}{2h}{〚u_h〛}_{x_n}^2+\frac{{\sigma }_N}{h}{〚u_h〛}_{x_N}^2\end{array}$ (82)

An upper bound to the term ${\displaystyle {\sum }_{n=0}^N{\left\{\left|K\left(x,\overline{u}\right){u^{\prime }}_h\right|\right\}}_{x_n}{〚u_h〛}_{x_n}}$ needs to be established to prove the coercivity of ${\tilde{a}}_h$ . Using Hypothesis (${\mathscr{H}}_n$ ) and definition of average:

$\begin{array}{l}\forall n\in \left\{1,\cdots ,N-1\right\}\\ \left|{\left\{\left|K\left(x,\overline{u}\right){u^{\prime }}_h\right|\right\}}_{x_n}\right|\le \frac{1}{2}\left(\left|K\left(x_n^{-},\overline{u}\left(x_n^{-}\right)\right){u^{\prime }}_h\left(x_n^{-}\right)\right|+\left|K\left(x_n^{+},\overline{u}\left(x_n^{+}\right)\right){u^{\prime }}_h\left(x_n^{+}\right)\right|\right)\\ \le \frac{K_1^{\left(n-1\right)}}{2}\left|{u^{\prime }}_h\left(x_n^{-}\right)\right|+\frac{K_1^{\left(n\right)}}{2}\left|{u^{\prime }}_h\left(x_n^{+}\right)\right|\end{array}$(83)

Recalling the trace inequality [48] in the case of an orthonormal polynomial basis: for an interval $I_n$ ,

$\forall u\in {\mathbb{P}}^p\left(I_n\right),\left|u\left(x_n^{+}\right)\right|\le C_{\text{tr},p}\frac{{\Vert u\Vert }_{L^2\left(I_n\right)}}{\sqrt{h}},\left|u\left(x_{n+1}^{-}\right)\right|\le C_{\text{tr},p}\frac{{\Vert u\Vert }_{L^2\left(I_n\right)}}{\sqrt{h}}$ (84)

we get, $\forall n\in \left\{1,\cdots ,N-1\right\}$ :

$\begin{array}{c}\left|{\left\{\left|K\left(x,\overline{u}\right){u^{\prime }}_h\right|\right\}}_{x_n}\right|\left|{〚u_h〛}_{x_n}\right|\le \left(\frac{K_1^{\left(n-1\right)}}{2}\frac{C_{\text{tr},p-1}^{\left(n-1\right)}}{\sqrt{h}}{\Vert {u^{\prime }}_h\Vert }_{I_{n-1}}+\frac{K_1^{\left(n\right)}}{2}\frac{C_{\text{tr},p-1}^{\left(n\right)}}{\sqrt{h}}{\Vert {u^{\prime }}_h\Vert }_{I_n}\right)\left|{〚u_h〛}_{x_n}\right|\\ \le \sqrt{{\epsilon }^{\left(n-1\right)}}\sqrt{K_0^{\left(n-1\right)}}{\Vert {u^{\prime }}_h\Vert }_{I_{n-1}}\frac{K_1^{\left(n-1\right)}}{2\sqrt{{\epsilon }^{\left(n-1\right)}}\sqrt{K_0^{\left(n-1\right)}}}\frac{C_{\text{tr},p-1}^{\left(n-1\right)}}{\sqrt{h}}\left|{〚u_h〛}_{x_n}\right|\\ +\sqrt{{\epsilon }^{\left(n\right)}}\sqrt{K_0^{\left(n\right)}}{\Vert {u^{\prime }}_h\Vert }_{I_n}\frac{K_1^{\left(n\right)}}{2\sqrt{{\epsilon }^{\left(n\right)}}\sqrt{K_0^{\left(n\right)}}}\frac{C_{\text{tr},p-1}^{\left(n\right)}}{\sqrt{h}}\left|{〚u_h〛}_{x_n}\right|\end{array}$ (85)

At the boundary nodes $x_0$ and $x_N$ , we have

$\left|{\left\{\left|K\left(x,\overline{u}\right){u^{\prime }}_h\right|\right\}}_{x_0}\right|\left|{〚u_h〛}_{x_0}\right|\le \sqrt{{\epsilon }^{\left(0\right)}}\sqrt{K_0^{\left(0\right)}}{\Vert {u^{\prime }}_h\Vert }_{I_0}\frac{K_1^{\left(0\right)}}{\sqrt{{\epsilon }^{\left(0\right)}}\sqrt{K_0^{\left(0\right)}}}\frac{C_{\text{tr},p-1}^{\left(0\right)}}{\sqrt{h}}\left|{〚u_h〛}_{x_0}\right|$(86)

$\left|{\left\{\left|K\left(x,\overline{u}\right){u^{\prime }}_h\right|\right\}}_{x_N}\right|\left|{〚u_h〛}_{x_N}\right|\le \sqrt{{\epsilon }^{\left(N-1\right)}}\sqrt{K_0^{\left(N-1\right)}}{\Vert {u^{\prime }}_h\Vert }_{I_{N-1}}\frac{K_1^{\left(N-1\right)}}{\sqrt{{\epsilon }^{\left(N-1\right)}}\sqrt{K_0^{\left(N-1\right)}}}\frac{C_{\text{tr},p-1}^{\left(N-1\right)}}{\sqrt{h}}\left|{〚u_h〛}_{x_N}\right|$ (87)

Gathering the bounds on the boundary and the interior nodes, we get

$\begin{array}{l}{\displaystyle \sum }_{n=0}^{N-1}\left|{\left\{\left|K\left(x,\overline{u}\right){u^{\prime }}_h\right|\right\}}_{x_n}\right|\left|{〚u_h〛}_{x_n}\right|\\ \le \sqrt{{\epsilon }^{\left(0\right)}}\sqrt{K_0^{\left(0\right)}}{\Vert {u^{\prime }}_h\Vert }_{I_0}\frac{K_1^{\left(0\right)}}{\sqrt{{\epsilon }^{\left(0\right)}}\sqrt{K_0^{\left(0\right)}}}\frac{C_{\text{tr},p-1}^{\left(0\right)}}{\sqrt{h}}\left|{〚u_h〛}_{x_0}\right|\\ +{\displaystyle \sum }_{n=1}^{N-1}(\sqrt{{\epsilon }^{\left(n-1\right)}}\sqrt{K_0^{\left(n-1\right)}}{\Vert {u^{\prime }}_h\Vert }_{I_{n-1}}\frac{K_1^{\left(n-1\right)}}{2\sqrt{{\epsilon }^{\left(n-1\right)}}\sqrt{K_0^{\left(n-1\right)}}}\frac{C_{\text{tr},p-1}^{\left(n-1\right)}}{\sqrt{h}}\left|{〚u_h〛}_{x_n}\right|\\ +\sqrt{{\epsilon }^{\left(n\right)}}\sqrt{K_0^{\left(n\right)}}{\Vert {u^{\prime }}_h\Vert }_{I_n}\frac{K_1^{\left(n\right)}}{2\sqrt{{\epsilon }^{\left(n\right)}}\sqrt{K_0^{\left(n\right)}}}\frac{C_{\text{tr},p-1}^{\left(n\right)}}{\sqrt{h}}\left|{〚u_h〛}_{x_n}\right|)\\ +\sqrt{{\epsilon }^{\left(N-1\right)}}\sqrt{K_0^{\left(N-1\right)}}{\Vert {u^{\prime }}_h\Vert }_{I_{N-1}}\frac{K_1^{\left(N-1\right)}}{\sqrt{{\epsilon }^{\left(N-A\right)}}\sqrt{K_0^{\left(N-1\right)}}}\frac{C_{\text{tr},p-1}^{\left(N-1\right)}}{\sqrt{h}}\left|{〚u_h〛}_{x_N}\right|\end{array}$ (88)

Then, using Cauchy-Schwarz’s and Young’s inequality, we have:

$\begin{array}{l}{\displaystyle \sum }_{n=0}^{N-1}{\left\{\left|K\left(x,\overline{u}\right){u^{\prime }}_h\right|\right\}}_{x_n}{〚u_h〛}_{x_n}\\ \le {\displaystyle \sum }_{n=0}^{N-1}\frac{{\epsilon }^{\left(n\right)}{K}_0^{\left(n\right)}}{2}{\Vert {u^{\prime }}_h\Vert }_{I_n}^2+\frac{{\left(K_1^{\left(0\right)}{C}_{\text{tr},p-1}^{\left(0\right)}\right)}^2}{{\epsilon }^{\left(0\right)}{K}_0^{\left(0\right)}}\frac{{〚u_h〛}_{x_0}^2}{h}+\frac{{\left(K_1^{\left(N-1\right)}{C}_{\text{tr},p-1}^{\left(N-1\right)}\right)}^2}{{\epsilon }^{\left(N-1\right)}{K}_0^{\left(N-1\right)}}\frac{{〚u_h〛}_{x_N}^2}{h}\\ +{\displaystyle \sum }_{n=1}^{N-1}\left(\frac{{\left(K_1^{\left(n-1\right)}{C}_{\text{tr},p-1}^{\left(n-1\right)}\right)}^2}{2{\epsilon }^{\left(n-1\right)}{K}_0^{\left(n-1\right)}}\frac{{〚u_h〛}_{x_n}^2}{2h}+\frac{{\left(K_1^{\left(n\right)}{C}_{\text{tr},p-1}^{\left(n\right)}\right)}^2}{2{\epsilon }^{\left(n\right)}{K}_0^{\left(n\right)}}\frac{{〚u_h〛}_{x_n}^2}{2h}\right)\end{array}$(89)

From the above inequality, we deduce a lower bound of ${\tilde{a}}_h\left(u_h,u_h;\overline{u}\right),\forall u_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$

(90)

where

$\{\begin{array}{l}{\sigma }_n^{\text{*}}=\frac{{\left(K_1^{\left(n\right)}{C}_{\text{tr},p-1}^{\left(n\right)}\right)}^2}{2{\epsilon }^{\left(n\right)}{K}_0^{\left(n\right)}},\forall n\in \left\{1,\cdots ,N-1\right\}\hfill \\ {\sigma }_0^{\text{*}}=\frac{{\left(K_1^{\left(0\right)}{C}_{\text{tr},p-1}^{\left(0\right)}\right)}^2}{{\epsilon }^{\left(0\right)}{K}_0^{\left(0\right)}}\hfill \\ {\sigma }_N^{\text{*}}=\frac{{\left(K_1^{\left(N-1\right)}{C}_{\text{tr},p-1}^{\left(N-1\right)}\right)}^2}{{\epsilon }^{\left(N-1\right)}{K}_0^{\left(N-1\right)}}\hfill \end{array}$ (91)

Finally, thanks to the inequality (90), ${\tilde{a}}_h$ (32) is coercive if

$\{\begin{array}{ll}{\epsilon }^{\left(n\right)}<2,\hfill & \forall n\in \left\{0,\cdots ,N-1\right\}\hfill \\ {\sigma }_n>{\sigma }_n^{\text{*}},\hfill & \forall n\in \left\{1,\cdots ,N-1\right\}\hfill \\ {\sigma }_0>{\sigma }_0^{\text{*}}\hfill & \hfill \\ {\sigma }_N>{\sigma }_N^{\text{*}}\hfill & \hfill \end{array}$ (92)

which ends the proof. □

Proof of Lemma 5. For a given $\overline{u}\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ , an upper bound for $\left|{\tilde{a}}_h\left(u_h,v_h;\overline{u}\right)\right|$ , $\forall u_h,v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ needs to be established in order to prove continuity of ${\tilde{a}}_h$ . Firstly, start bounding above the volume contribution using Hypothesis (${\mathscr{H}}_n$ ):

$\begin{array}{c}\left|{\displaystyle \sum }_{n=0}^{N-1}{\displaystyle {\int }_{I_n}K\left(x,\overline{u}\right){u^{\prime }}_h{v^{\prime }}_h\text{d}x}\right|\le {\displaystyle \sum }_{n=0}^{N-1}{K}_1^{\left(n\right)}\left|{\displaystyle {\int }_{I_n}{u^{\prime }}_h{v^{\prime }}_h\text{d}x}\right|\le {\displaystyle \sum }_{n=0}^{N-1}\sqrt{K_1^{\left(n\right)}}{\Vert {u^{\prime }}_h\Vert }_{I_n}\sqrt{K_1^{\left(n\right)}}{\Vert {v^{\prime }}_h\Vert }_{I_n}\\ \le {\left({\displaystyle \sum }_{n=0}^{N-1}{K}_1^{\left(n\right)}{\Vert {u^{\prime }}_h\Vert }_{I_n}^2\right)}^{\frac{1}{2}}{\left({\displaystyle \sum }_{n=0}^{N-1}{K}_1^{\left(n\right)}{\Vert {v^{\prime }}_h\Vert }_{I_n}^2\right)}^{\frac{1}{2}}.\end{array}$(93)

Then, penalization terms are bounded above

$\begin{array}{l}\left|\frac{{\sigma }_0}{h}{〚u_h〛}_{x_0}{〚v_h〛}_{x_0}+{\displaystyle \sum }_{n=1}^{N-1}\frac{{\sigma }_{n-1}+{\sigma }_n}{2h}{〚u_h〛}_{x_n}{〚v_h〛}_{x_n}+\frac{{\sigma }_N}{h}{〚u_h〛}_{x_N}{〚v_h〛}_{x_N}\right|\\ \le {\left(\frac{{\sigma }_0}{h}{〚u_h〛}_{x_0}^2+{\displaystyle \sum }_{n=1}^{N-1}\frac{{\sigma }_{n-1}+{\sigma }_n}{2h}{〚u_h〛}_{x_n}^2+\frac{{\sigma }_N}{h}{〚u_h〛}_{x_N}^2\right)}^{\frac{1}{2}}\end{array}$ (94)

$\begin{array}{l}{\left(\frac{{\sigma }_0}{h}{〚v_h〛}_{x_0}^2+{\displaystyle \sum }_{n=1}^{N-1}\frac{{\sigma }_{n-1}+{\sigma }_n}{2h}{〚v_h〛}_{x_n}^2+\frac{{\sigma }_N}{h}{〚v_h〛}_{x_N}^2\right)}^{\frac{1}{2}}\\ \le \text{max}\left({\sigma }_0,{\sigma }_N,\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{max}}\left(\frac{{\sigma }_n}{2}\right)\right)\Vert u_h\Vert \Vert v_h\Vert \end{array}$ (95)

and one can write

$\begin{array}{l}{\displaystyle \sum }_{n=0}^N\left|{\left\{\left|K\left(x,\overline{u}\right){u^{\prime }}_h\right|\right\}}_{x_n}{〚v_h〛}_{x_n}\right|\\ \le {\left(2{\displaystyle \sum }_{n=0}^{N-1}{\epsilon }^{\left(n\right)}{K}_0^{\left(n\right)}{\Vert {u^{\prime }}_h\Vert }_{I_n}^2\right)}^{\frac{1}{2}}(\frac{{\left(K_1^{\left(0\right)}{C}_{\text{tr},p-1}^{\left(0\right)}\right)}^2}{{\epsilon }^{\left(0\right)}{K}_0^{\left(0\right)}}\frac{{〚v_h〛}_{x_0}^2}{h}+\frac{{\left(K_1^{\left(N-1\right)}{C}_{\text{tr},p-1}^{\left(N-1\right)}\right)}^2}{{\epsilon }^{\left(N-1\right)}{K}_0^{\left(N-1\right)}}\frac{{〚v_h〛}_{x_N}^2}{h}\\ +{{\displaystyle \sum }_{n=1}^{N-1}\left(\frac{{\left(K_1^{\left(n-1\right)}{C}_{\text{tr},p-1}^{\left(n-1\right)}\right)}^2}{2{\epsilon }^{\left(n-1\right)}{K}_0^{\left(n-1\right)}}+\frac{{\left(K_1^{\left(n\right)}{C}_{\text{tr},p-1}^{\left(n\right)}\right)}^2}{2{\epsilon }^{\left(n\right)}{K}_0^{\left(n\right)}}\right)\frac{v_h{}_{x_n}^2}{2h})}^{\frac{1}{2}}.\end{array}$(96)

From those inequalities, we obtain an upper bound $\forall u_h,v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ , as follows

$\begin{array}{c}\left|{\tilde{a}}_h\left(u_h,v_h;\overline{u}\right)\right|\le {\left({\displaystyle \sum }_{n=0}^{N-1}{K}_1^{\left(n\right)}{\Vert {u^{\prime }}_h\Vert }_{I_n}^2\right)}^{\frac{1}{2}}{\left({\displaystyle \sum }_{n=0}^{N-1}{K}_1^{\left(n\right)}{\Vert {v^{\prime }}_h\Vert }_{I_n}^2\right)}^{\frac{1}{2}}\\ +{\left({\displaystyle \sum }_{n=0}^{N-1}2{\epsilon }^{\left(n\right)}{K}_1^{\left(n\right)}{\Vert {u^{\prime }}_h\Vert }_{I_n}^2\right)}^{\frac{1}{2}}{\left(\frac{{\sigma }_0^{\text{*}}}{h}{〚v_h〛}_{x_0}^2+{\displaystyle \sum }_{n=1}^{N-1}\frac{{\sigma }_{n-1}^{\text{*}}+{\sigma }_n^{\text{*}}}{2h}{〚v_h〛}_{x_n}^2+\frac{{\sigma }_N^{\text{*}}}{h}{〚v_h〛}_{x_N}^2\right)}^{\frac{1}{2}}\\ +\max \left({\sigma }_0,{\sigma }_N,\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\max }\left(\frac{{\sigma }_n}{2}\right)\right)\Vert u_h\Vert \Vert v_h\Vert \\ \le \underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\max }\left(K_1^{\left(n\right)}\right)\Vert u_h\Vert \Vert v_h\Vert \\ +\sqrt{\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{max}}\left(2{\epsilon }^{\left(n\right)}{K}_1^{\left(n\right)}\right)\text{max}\left({\sigma }_0^{\text{*}},{\sigma }_N^{\text{*}},\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{max}}\left(\frac{{\sigma }_n^{\text{*}}}{2}\right)\right)}\Vert u_h\Vert \Vert v_h\Vert \\ +\max \left({\sigma }_0,{\sigma }_N,\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\max }\left(\frac{{\sigma }_n}{2}\right)\right)\Vert u_h\Vert \Vert v_h\Vert \\ \le \tilde{C}\Vert u_h\Vert \Vert v_h\Vert \end{array}$ (97)

where

$\begin{array}{c}\tilde{C}\left(ϵ\right)=\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{max}}\left(K_1^{\left(n\right)}\right)+\sqrt{\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{max}}\left(2{\epsilon }^{\left(n\right)}{K}_1^{\left(n\right)}\right)\text{max}\left({\sigma }_0^{\text{*}},{\sigma }_N^{\text{*}},\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\text{max}}\left(\frac{{\sigma }_n^{\text{*}}}{2}\right)\right)}\\ +\max \left({\sigma }_0,{\sigma }_N,\underset{\begin{array}{c}n=1,\cdots ,N-1\end{array}}{\max }\left(\frac{{\sigma }_n}{2}\right)\right).\end{array}$(98)

□

Proof of Lemma 6. An upper bound for $\left|l\left(v_h\right)\right|$ is established using Poincaré inequality and Cauchy Schwarz: $\forall v_h\in {\mathcal{V}}_0^p{(}_h),$

$\begin{array}{c}\left|{\displaystyle \sum }_{n=0}^{N-1}{\displaystyle {\int }_{I_n}fv\text{d}x}\right|\le {\displaystyle \sum }_{n=0}^{N-1}{\Vert f\Vert }_{I_n}{\Vert v_h\Vert }_{I_n}\le {\displaystyle \sum }_{n=0}^{N-1}{\Vert f\Vert }_{I_n}{\beta }_n{\Vert {v^{\prime }}_h\Vert }_{I_n}\\ \le {\left({\displaystyle \sum }_{n=0}^{N-1}{\Vert f\Vert }_{I_n}^2\right)}^{\frac{1}{2}}{\left({\displaystyle \sum }_{n=0}^{N-1}{\beta }_n^2{\Vert {v^{\prime }}_h\Vert }_{I_n}^2\right)}^{\frac{1}{2}}\le B\Vert v_h\Vert \end{array}$ (99)

with $B=\underset{\begin{array}{c}n=0,\cdots ,N-1\end{array}}{\text{max}}\left({\beta }_n\right){\left({\displaystyle {\sum }_{n=0}^{N-1}{\Vert f\Vert }_{I_n}^2}\right)}^{\frac{1}{2}}$ . □

Proof of Lemma 12. For a given $\overline{u}\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ , let $v_h$ be a sequence in ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ bounded in the $\Vert .\Vert$ -norm and let $\phi \in C_0^{\infty }\left(\Omega \right)$ . For all $h\in {ℝ}_{+}^{\text{*}}$ , set ${\phi }_h={\pi }_h\phi$ where ${\pi }_h$ denotes the $L^2$ -orthogonal projection onto ${\mathcal{V}}_0^p\left({ℰ}_h\right)$ . Since $p\ge 1$ , infer $\Vert \phi -{\pi }_h\phi \Vert \underset{h\to 0}{\to }0$ . Owing to Equation (67) and since $G_h^p\left(\phi \right)={\phi }^{\prime }$ because $\phi \in C_0^{\infty }\left(\Omega \right)$ , obtain for all $p\ge 0$

$G_h^p\left({\pi }_h\phi \right)\to {\phi }^{\prime }\text{strongly}\text{in}{L}^{\text{2}}\left(\Omega \right)$ (100)

One can observe that

${\tilde{a}}_h\left(v_h,{\pi }_h\phi ;\overline{u}\right)={\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(v_h\right)G_h^p\left({\pi }_h\phi \right)\text{d}x+s_h\left(v_h,{\pi }_h\phi \right):={\mathfrak{T}}_1+{\mathfrak{T}}_2$(101)

Clearly as $h\to 0$ , ${\mathfrak{T}}_1\to {\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right)v^{\prime }{\phi }^{\prime }\text{d}x$ owing to the weak convergence of ${\widehat{G}}_h^p\left(v_h\right)$ to $v^{\prime }$ and to the strong convergence of $G_h^p\left({\pi }_h\phi \right)$ to ${\phi }^{\prime }$ . Furthermore, using Cauchy-Schwarz inequality yields:

$\begin{array}{c}\left|{\mathfrak{T}}_2\right|=\left|s_h\left(v_h,{\pi }_h\phi \right)\right|\\ =\left|\frac{{\sigma }_0}{h}{〚v_h〛}_{x_0}{〚{\pi }_h\phi 〛}_{x_0}+{\displaystyle \sum }_{n=1}^{N-1}\frac{{\sigma }_{n-1}+{\sigma }_n}{2h}{〚v_h〛}_{x_n}{〚{\pi }_h\phi 〛}_{x_n}+\frac{{\sigma }_N}{h}{〚v_h〛}_{x_N}{〚{\pi }_h\phi 〛}_{x_N}\right|\\ \le {\left({\sigma }_0^2\frac{{〚v_h〛}_{x_0}^2}{h}+{\displaystyle \sum }_{n=1}^{N-1}\frac{{\left({\sigma }_{n-1}+{\sigma }_n\right)}^2}{4}\frac{{〚v_h〛}_{x_n}^2}{h}+{\sigma }_N^2\frac{{〚v_h〛}_{x_N}^2}{h}\right)}^{\frac{1}{2}}{\left({\displaystyle \sum }_{n=0}^N\frac{{〚{\pi }_h\phi 〛}_{x_n}^2}{h}\right)}^{\frac{1}{2}}\\ \le \mathcal{C}{\left|v_h\right|}_J{\left|{\pi }_h\phi \right|}_J\end{array}$ (102)

where

$\mathcal{C}=\text{max}\left\{{\sigma }_0^2,\frac{{\left({\sigma }_{n-1}+{\sigma }_n\right)}^2}{4},{\sigma }_N^2\right\}.$ (103)

Since ${\left|v_h\right|}_J$ is bounded by assumption and since ${\left|{\pi }_h\phi \right|}_J={\left|\phi -{\pi }_h\phi \right|}_J\underset{h\to 0}{\to }0$ , infer ${\mathfrak{T}}_2\underset{h\to 0}{\to }0$ . □

Proof of the Theorem 13. For a given $\overline{u}\in {\mathcal{V}}_0^p\left({ℰ}_h\right)$ , owing to the discrete coercivity of ${\tilde{a}}_h$ , the sequence $u_h$ is bounded in the $\Vert .\Vert$ -norm. Theorem 11 implies that there is $v\in H_0^1\left(\Omega \right)$ such that up to a subsequence, $u_h\to v$ in $L^{\text{2}}\left(\Omega \right)$ and for all $p\ge 0$ , $G_h^p\left(u_h\right)\rightharpoonup v^{\prime }$ weakly in $L^{\text{2}}\left(\Omega \right)$ as $h\to 0$ . Let $\phi \in C_0^{\infty }\left(\Omega \right)$ . Owing Lemma 12, ${\tilde{a}}_h\left(u_h,{\pi }_h\phi ;\overline{u}\right)\to \tilde{a}\left(v,\phi \right)$ as $h\to 0$ . Since $u_h$ solves the discrete linearized problem $\left({\tilde{\mathcal{W}}}_h\right)$ , infer as $h\to 0$

$\begin{array}{c}{\tilde{a}}_h\left(u_h-{\pi }_h\phi ,u_h-{\pi }_h\phi ;\overline{u}\right)={\tilde{a}}_h\left(u_h,u_h-{\pi }_h\phi ;\overline{u}\right)-{\tilde{a}}_h\left({\pi }_h\phi ,u_h-{\pi }_h\phi ;\overline{u}\right)\\ \to \tilde{a}\left(v,v-\phi \right)-\tilde{a}\left(\phi ,v-\phi \right)\\ \to {\displaystyle {\int }_{\Omega }\left(v-\phi \right)f\text{d}x}-{\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\phi }^{\prime }{\left(v-\phi \right)}^{\prime }\text{d}x\end{array}$(104)

Hence, using ${\tilde{a}}_h\left(v_h,v_h;\overline{u}\right)\ge C^{\text{*}}{\Vert v_h\Vert }^2$ from Lemma 4

$\begin{array}{l}{C}^{\text{*}}\Vert u_h-{\pi }_h\phi \Vert \le {\tilde{a}}_h\left(u_h-{\pi }_h\phi ,u_h-{\pi }_h\phi ;\overline{u}\right)\\ \iff \underset{h\to 0}{\lim \sup }{C}^{*}\Vert u_h-{\pi }_h\phi \Vert \le \underset{h\to 0}{\lim \sup }{\tilde{a}}_h\left(u_h-{\pi }_h\phi ,u_h-{\pi }_h\phi ;\overline{u}\right)\\ \le \left|{\displaystyle {\int }_{\Omega }\left(v-\phi \right)f\text{d}x}-{\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\phi }^{\prime }{\left(v-\phi \right)}^{\prime }\text{d}x\right|\\ \le {\Vert f\Vert }_{L^{\text{2}}\left(\Omega \right)}{\Vert v-\phi \Vert }_{L^{\text{2}}\left(\Omega \right)}+K_1{\Vert {\phi }^{\prime }\Vert }_{L^{\text{2}}\left(\Omega \right)}{\Vert {\left(v-\phi \right)}^{\prime }\Vert }_{L^{\text{2}}\left(\Omega \right)}\\ \le C_{f,\phi }{\left({\Vert v-\phi \Vert }_{L^{\text{2}}\left(\Omega \right)}^2+{\Vert {\left(v-\phi \right)}^{\prime }\Vert }_{L^{\text{2}}\left(\Omega \right)}^2\right)}^{\frac{1}{2}}\\ \le C_{f,\phi }{\Vert v-\phi \Vert }_{H^1\left(\Omega \right)}\end{array}$(105)

with $C_{f,\phi }={\left({\Vert f\Vert }_{L^{\text{2}}\left(\Omega \right)}^2+K_1{\Vert {\phi }^{\prime }\Vert }_{L^{\text{2}}\left(\Omega \right)}^2\right)}^{\frac{1}{2}}$ . As a consequence,

$\underset{h\to 0}{\lim \sup }\Vert u_h-{\pi }_h\phi \Vert \le \frac{1}{C^{*}}{C}_{f,\phi }{\Vert v-\phi \Vert }_{H^1\left(\Omega \right)}.$ (106)

One can observe that the choice for ${\widehat{G}}_h^p$ satisfy the stability property

$\forall v_h\in {\mathcal{V}}_0^p\left({ℰ}_h\right),{\Vert {\widehat{G}}_h^p\left(v_h\right)\Vert }_{L^{\text{2}}\left(\Omega \right)}\le \widehat{C}\Vert v_h\Vert$(107)

for $\widehat{C}$ independent of $h$ . As a result,

$\underset{h\to 0}{\lim \sup }{\Vert {\widehat{G}}_h^p\left(u_h\right)-{\widehat{G}}_h^p\left({\pi }_h\phi \right)\Vert }_{L^2\left(\Omega \right)}\le \widehat{C}\frac{1}{C^{*}}{C}_{f,\phi }{\Vert v-\phi \Vert }_{H^1\left(\Omega \right)}$ (108)

because

$\begin{array}{l}{\Vert {\widehat{G}}_h^p\left(u_h\right)-{\widehat{G}}_h^p\left({\pi }_h\phi \right)\Vert }_{L^{\text{2}}\left(\Omega \right)}\le \widehat{C}\Vert u_h-{\pi }_h\phi \Vert \\ \iff \underset{h\to 0}{\lim \sup }{\Vert {\widehat{G}}_h^p\left(u_h\right)-{\widehat{G}}_h^p\left({\pi }_h\phi \right)\Vert }_{L^2\left(\Omega \right)}\le \widehat{C}\underset{h\to 0}{\lim \sup }\Vert u_h-{\pi }_h\phi \Vert \\ \le \widehat{C}\frac{1}{C^{*}}{C}_{f,\phi }{\Vert v-\phi \Vert }_{H^1\left(\Omega \right)}\end{array}$(109)

And since ${\widehat{G}}_h^p\left({\pi }_h\phi \right)$ strongly converges to ${\phi }^{\prime }$ in $L^{\text{2}}\left(\Omega \right)$ , this yields

$\underset{h\to 0}{\lim \sup }{\Vert {\widehat{G}}_h^p\left(u_h\right)-{\phi }^{\prime }\Vert }_{L^2\left(\Omega \right)}\le \widehat{C}\frac{1}{C^{*}}{C}_{f,\phi }{\Vert v-\phi \Vert }_{H^1\left(\Omega \right)}.$ (110)

Since $\phi$ is arbitrary in $C_0^{\infty }\left(\Omega \right)$ , and since this space is dense in $H_0^1\left(\Omega \right)$ , the term on the right hand side can be made as small as desired taking $\phi =v$ , infer

${\widehat{G}}_h^p\left(u_h\right)\underset{h\to 0}{\to }{v}^{\prime }\text{strongly}\text{in}{L}^{\text{2}}\left(\Omega \right)$(111)

As a result, taking $\phi$ arbitrary in $C_0^{\infty }\left(\Omega \right)$ yields

${\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right)v^{\prime }{\phi }^{\prime }\text{d}x\underset{h\to 0}{\leftarrow }{\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){u^{\prime }}_h{\pi }_h{\phi }^{\prime }\text{d}x={\tilde{a}}_h\left(u_h,{\pi }_h\phi \right)={\displaystyle {\int }_{\Omega }f{\pi }_h\phi }\underset{h\to 0}{\to }{\displaystyle {\int }_{\Omega }f\phi \text{d}x}$ (112)

using Lemma 12., i.e., $v$ solves the Poisson problem by density of $C_0^{\infty }\left(\Omega \right)$ in $H^1\left(\Omega \right)$ . Since the solution $u$ to the Poisson problem is unique, the whole sequence $u_h$ strongly converges to $u$ in $L^{\text{2}}\left(\Omega \right)$ and, for all $p\ge 0$ , the sequence ${\left(G_h^p\left(u_h\right)\right)}_{h\in {ℝ}_{+}^{\text{*}}}$ weakly converges to $u^{\prime }$ in $L^{\text{2}}\left(\Omega \right)$ .

$\begin{array}{c}{\tilde{a}}_h\left(u_h,u_h;\overline{u}\right)={\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right)G_h^p\left(u_h\right)\text{d}x+s_h\left(u_h,u_h\right)\\ \ge {\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right)G_h^p\left(u_h\right)\text{d}x\end{array}$(113)

Thus

$\underset{h\to 0}{\lim \inf }{\tilde{a}}_h\left(u_h,u_h;\overline{u}\right)\ge \underset{h\to 0}{\lim \inf }{\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right)G_h^p\left(u_h\right)\text{d}x\ge {\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right)u^{\prime }{u}^{\prime }\text{d}x$ (114)

Furthermore

${\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right)G_h^p\left(u_h\right)\text{d}x\le {\tilde{a}}_h\left(u_h,u_h;\overline{u}\right)={\displaystyle {\int }_{\Omega }f{u}_h\text{d}x}$(115)

yielding with Equation (112)

$\begin{array}{c}\underset{h\to 0}{\lim \sup }{\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right)G_h^p\left(u_h\right)\text{d}x\le \underset{h\to 0}{\lim \sup }{\tilde{a}}_h\left(u_h,u_h;\overline{u}\right)\\ =\underset{h\to 0}{\lim \sup }{\displaystyle {\int }_{\Omega }f{u}_h\text{d}x}\le {\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right)u^{\prime }{u}^{\prime }\text{d}x\end{array}$ (116)

Thus, ${\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right)G_h^p\left(u_h\right)\text{d}x\underset{h\to 0}{\to }{\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right)u^{\prime }{u}^{\prime }\text{d}x$ strongly. Moreover, ${\tilde{a}}_h\left(u_h,u_h;\overline{u}\right)\underset{h\to 0}{\to }{\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right)u^{\prime }{u}^{\prime }\text{d}x$ strongly. Owing that

$\begin{array}{l}{\tilde{a}}_h\left(u_h,u_h;\overline{u}\right)={\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right)G_h^p\left(u_h\right)\text{d}x+s_h\left(u_h,u_h\right)\\ \ge {\displaystyle {\int }_{\Omega }K}\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right)G_h^p\left(u_h\right)\text{d}x+\underset{\begin{array}{c}n=0,\cdots ,N\end{array}}{\text{min}}\left({\sigma }_n\right){\left|u_h\right|}_J^2\\ \iff \underset{\begin{array}{c}n=0,\cdots ,N\end{array}}{\text{min}}\left({\sigma }_n\right){\left|u_h\right|}_J^2\le a_h\left(u_h,u_h\right)-{\displaystyle {\int }_{\Omega }K\left(x,\overline{u}\right){\widehat{G}}_h^p\left(u_h\right)G_h^p\left(u_h\right)\text{d}x}\end{array}$(117)

and since $\underset{\begin{array}{c}n=0,\cdots ,N\end{array}}{\text{min}}\left({\sigma }_n\right)>0$ and the right-hand side tends to zero, ${\left|u_h\right|}_J\to 0$ .

${\Vert {u^{\prime }}_h-u^{\prime }\Vert }_{L^{\text{2}}\left(\Omega \right)}={\Vert {\widehat{G}}_h^p\left(u_h\right)-u^{\prime }\Vert }_{L^{\text{2}}\left(\Omega \right)}\to 0,$(118)

□

Conflicts of Interest

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

References

[1]	Darcy, H. (1856) Les Fontaines Publiques de La Ville de Dijon. Librairie des Corps Impériaux des Ponts et des Chaussées et des Mines.
[2]	Buckingham, E. (1907) Studies on the Movement of Soil Moisture, Volume 38. Bulletin Edition, Washington, Government Publishing Office.
[3]	Richardson, L.F. (1922) Weather Prediction by Numerical Process. Cambridge University Press.
[4]	Richards, L.A. (1931) Capillary Conduction of Liquids through Porous Mediums. Physics, 1, 318-333.[CrossRef]
[5]	Rubin, J. (1968) Theoretical Analysis of Two-Dimensional, Transient Flow of Water in Unsaturated and Partly Unsaturated Soils. Soil Science Society of America Journal, 32, 607-615.[CrossRef]
[6]	Cooley, R.L. (1971) A Finite Difference Method for Unsteady Flow in Variably Saturated Porous Media: Application to a Single Pumping Well. Water Resources Research, 7, 1607-1625.[CrossRef]
[7]	Clément, J.-B. (2021) Numerical Simulation of Flows in Unsaturated Porous Media by an Adaptive Discontinuous Galerkin Method: Application to Sandy Beaches. Ph.D. Thesis, Université de Toulon.
[8]	Irmay, S. (1954) On the Hydraulic Conductivity of Unsaturated Soils. Eos, Transactions American Geophysical Union, 35, 463-467.
[9]	Vachaud, G. and Thony, J. (1971) Hysteresis during Infiltration and Redistribution in a Soil Column at Different Initial Water Contents. Water Resources Research, 7, 111-127.[CrossRef]
[10]	van Genuchten, M.T. (1980) A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Science Society of America Journal, 44, 892-898.[CrossRef]
[11]	Dolejší, V., Kuraz, M. and Solin, P. (2019) Adaptive Higher-Order Space-Time Discontinuous Galerkin Method for the Computer Simulation of Variably-Saturated Porous Media Flows. Applied Mathematical Modelling, 72, 276-305.[CrossRef]
[12]	Clément, J., Golay, F., Ersoy, M. and Sous, D. (2021) An Adaptive Strategy for Discontinuous Galerkin Simulations of Richards’ Equation: Application to Multi-Materials Dam Wetting. Advances in Water Resources, 151, Article ID: 103897.[CrossRef]
[13]	Pietro, D.A.D. and Ern, A. (2012) Mathematical Aspects of Discontinuous Galerkin Methods. Springer.
[14]	Rivière, B. (2008) Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations. Society for Industrial and Applied Mathematics.[CrossRef]
[15]	Shin, H.-G. (2025) A Posteriori Error Estimates and Adaptivity for Solving the Richards Equation. Ph.D. Thesis, Charles University. https://dspace.cuni.cz/bitstream/handle/20.500.11956/198009/140128914.pdf?sequence=1
[16]	Congreve, S., Dolejší, V. and Sakić, S. (2024) Error Analysis for Local Discontinuous Galerkin Semidiscretization of Richards’ Equation. IMA Journal of Numerical Analysis, 45, 580-630.[CrossRef]
[17]	Epshteyn, Y. and Rivière, B. (2007) Estimation of Penalty Parameters for Symmetric Interior Penalty Galerkin Methods. Journal of Computational and Applied Mathematics, 206, 843-872.[CrossRef]
[18]	Süli, E. and Mayers, D.F. (2003) An Introduction to Numerical Analysis. Cambridge University Press.[CrossRef]
[19]	Li, H., Farthing, M.W., Dawson, C.N. and Miller, C.T. (2007) Local Discontinuous Galerkin Approximations to Richards’ Equation. Advances in Water Resources, 30, 555-575.[CrossRef]
[20]	Farthing, M.W., Kees, C.E. and Miller, C.T. (2003) Mixed Finite Element Methods and Higher Order Temporal Approximations for Variably Saturated Groundwater Flow. Advances in Water Resources, 26, 373-394.[CrossRef]
[21]	Hay, A., Etienne, S., Pelletier, D. and Garon, A. (2015) Hp-Adaptive Time Integration Based on the BDF for Viscous Flows. Journal of Computational Physics, 291, 151-176.[CrossRef]
[22]	Dahlquist, G.G. (1963) A Special Stability Problem for Linear Multistep Methods. BIT Numerical Mathematics, 3, 27-43.[CrossRef]
[23]	Hairer, E. and Wanner, G. (1996) Solving Ordinary Differential Equations II, Volume 14 of Springer Series in Computational Mathematics. Springer.
[24]	Skelboe, S. (1977) The Control of Order and Steplength for Backward Differentiation Methods. BIT, 17, 91-107.[CrossRef]
[25]	Lehmann, F. and Ackerer, P. (1998) Comparison of Iterative Methods for Improved Solutions of the Fluid Flow Equation in Partially Saturated Porous Media. Transport in Porous Media, 31, 275-292.[CrossRef]
[26]	List, F. and Radu, F.A. (2016) A Study on Iterative Methods for Solving Richards’ Equation. Computational Geosciences, 20, 341-353.[CrossRef]
[27]	Costa-Solé, A., Ruiz-Gironés, E. and Sarrate, J. (2021) High-Order Hybridizable Discontinuous Galerkin Formulation with Fully Implicit Temporal Schemes for the Simulation of Two-Phase Flow through Porous Media. International Journal for Numerical Methods in Engineering, 122, 3583-3612.[CrossRef]
[28]	Farthing, M.W. and Ogden, F.L. (2017) Numerical Solution of Richards’ Equation: A Review of Advances and Challenges. Soil Science Society of America Journal, 81, 1257-1269.[CrossRef]
[29]	Bergamaschi, L. and Putti, M. (1999) Mixed Finite Elements and Newton-Type Linearizations for the Solution of Richards’ Equation. International Journal for Numerical Methods in Engineering, 45, 1025-1046.[CrossRef]
[30]	Miller, C.T., Abhishek, C. and Farthing, M.W. (2006) A Spatially and Temporally Adaptive Solution of Richards’ Equation. Advances in Water Resources, 29, 525-545.[CrossRef]
[31]	Brad Thoms, R., Johnson, R.L. and Healy, R.W. (2006) User’s Guide to the Variably Saturated Flow (VSF) Process to MODFLOW.[CrossRef]
[32]	Ersoy, M., Golay, F. and Yushchenko, L. (2013) Adaptive Multiscale Scheme Based on Numerical Density of Entropy Production for Conservation Laws. Open Mathematics, 11, 1392-1415.[CrossRef]
[33]	Golay, F., Ersoy, M., Yushchenko, L. and Sous, D. (2015) Block-Based Adaptive Mesh Refinement Scheme Using Numerical Density of Entropy Production for Three-Dimensional Two-Fluid Flows. International Journal of Computational Fluid Dynamics, 29, 67-81.[CrossRef]
[34]	Altazin, T., Ersoy, M., Golay, F., Sous, D. and Yushchenko, L. (2016) Numerical Investigation of BB-AMR Scheme Using Entropy Production as Refinement Criterion. International Journal of Computational Fluid Dynamics, 30, 256-271.[CrossRef]
[35]	Boccardo, L., Thierry, G. and Murat, F. (1992) Unicité de la solution de certaines équations elliptiques non linéaires. Comptes rendus de l’Académie des Sciences Paris, 315, 1159-1164.
[36]	Bassi, F. and Rebay, S. (1997) A High-Order Accurate Discontinuous Finite Element Method for the Numerical Solution of the Compressible Navier-Stokes Equations. Journal of Computational Physics, 131, 267-279.[CrossRef]
[37]	Brezzi, F., Manzini, G., Marini, D., Pietra, P. and Russo, A. (2000) Discontinuous Galerkin Approximations for Elliptic Problems. Numerical Methods for Partial Differential Equations, 16, 365-378.[CrossRef]
[38]	Haverkamp, R., Vauclin, M., Touma, J., Wierenga, P.J. and Vachaud, G. (1977) A Comparison of Numerical Simulation Models for One-Dimensional Infiltration. Soil Science Society of America Journal, 41, 285-294.[CrossRef]
[39]	Vauclin, M., Khanji, D. and Vachaud, G. (1979) Experimental and Numerical Study of a Transient, Two-Dimensional Unsaturated-Saturated Water Table Recharge Problem. Water Resources Research, 15, 1089-1101.[CrossRef]
[40]	Dolejší, V. and Feistauer, M. (2015) Discontinuous Galerkin Method. Springer International Publishing.
[41]	Celia, M.A., Bouloutas, E.T. and Zarba, R.L. (1990) A General Mass-Conservative Numerical Solution for the Unsaturated Flow Equation. Water Resources Research, 26, 1483-1496.[CrossRef]
[42]	Philip, J.R. (2006) The Theory of Infiltration. 1. The Infiltration Equation and Its Solution. Soil Science, 171, S34-S46.[CrossRef]
[43]	Sochala, P. (2008) Méthodes Numériques Pour Les Écoulements Souterrains et Couplage Avec Le Ruissellement. Ph.D. Thesis, Paris Est.
[44]	Manzini, G. and Ferraris, S. (2004) Mass-Conservative Finite Volume Methods on 2-D Unstructured Grids for the Richards’ Equation. Advances in Water Resources, 27, 1199-1215.[CrossRef]
[45]	Dogan, A. and Motz, L.H. (2005) Saturated-Unsaturated 3D Groundwater Model. II: Verification and Application. Journal of Hydrologic Engineering, 10, 505-515.[CrossRef]
[46]	Twarakavi, N.K.C., Šimůnek, J. and Seo, S. (2008) Evaluating Interactions between Groundwater and Vadose Zone Using the Hydrus-Based Flow Package for Modflow. Vadose Zone Journal, 7, 757-768.[CrossRef]
[47]	Zha, Y., Shi, L., Ye, M. and Yang, J. (2013) A Generalized Ross Method for Two-and Three-Dimensional Variably Saturated Flow. Advances in Water Resources, 54, 67-77.[CrossRef]
[48]	Warburton, T. and Hesthaven, J.S. (2003) On the Constants in Hp-Finite Element Trace Inverse Inequalities. Computer Methods in Applied Mechanics and Engineering, 192, 2765-2773.[CrossRef]