Elyes Ahmed
email: elyes.ahmed@sintef.no

Saber Amdouni
email: saber.amdouni@enit.utm.tn

Efficient Adaptive Solvers for Doubly Degenerate Elliptic Equations using A Posteriori Error Estimates

Keywords: degenerate elliptic equations, fast and slow diffusion, a posteriori error analysis, local efficiency, adaptive linearization, switching criteria

In this work, we derive a posteriori error estimates for a class of doubly nonlinear and degenerate elliptic equations, including the Stefan problem and fast and slow diffusion in porous media. Our approach employs equilibrated flux reconstructions, providing guaranteed and fully computable upper bounds on an energy-type norm and local efficiency. These bounds remain independent of the strength of nonlinearity and degeneracy rates. These estimators drive an adaptive solver, dynamically switching between nonlinear solvers to achieve optimal iterations. The adaptive algorithm accounts for discretization, regularization, quadrature, and linearization error components. When Newton's method encounters challenges in achieving convergence, the adaptive algorithm transitions to the L-scheme solver. This solver optimally precomputes the stabilization (or tuning) parameter L > 0 during an offline phase, mirroring the behavior of the Jacobian. The adaptive algorithm is exemplified through four prototypical examples, showcasing its effective error control and notable computational savings.

Introduction

In this article, we present an adaptive solver designed to approximate doubly degenerate elliptic equations of the Stefan and porous media forms, characterized by the expression β(u) -∆ϕ(u) = f + ∇•ξ. Our paper focuses specifically on addressing scenarios where both ϕ and β lack uniform Lipschitz continuity. This is particularly relevant in cases such as the fast diffusion model [START_REF] Droniou | The gradient discretisation method[END_REF], which is non-Lipschitz at infinity, or when these functions are not Lipschitz-continuous at specific points, as observed in slow diffusion. We also allow for the possibility of plateaus in these functions. These degenerate elliptic equations hold significant relevance in various fields, particularly in porous media applications such as density-dependent groundwater flow modeling [START_REF] Chen | Fully discrete finite element analysis of multiphase flow in groundwater hydrology[END_REF] and immiscible two-phase flow problems [START_REF] Galusinski | On a degenerate parabolic system for compressible, immiscible, twophase flows in porous media[END_REF][START_REF] Girault | A finite element method for degenerate two-phase flow in porous media. part i: Well-posedness[END_REF]. Additionally, they find crucial applications in biology [START_REF] Epstein | Degenerate Diffusion Operators Arising in Population Biology[END_REF][START_REF] Gurtin | On the diffusion of biological populations[END_REF] chemistry [START_REF] Truskey | Transport phenomena in biological systems[END_REF], and material sciences [START_REF] Bonetti | Local existence for frémond's model of damage in elastic materials[END_REF]. It is of utmost importance to highlight that the solution operator of this model not only showcases nonlinearity and degeneracy but also may displays a striking transition between two distinct degeneration rates, which is contingent on the specific choices of ϕ and β. Consequently, the diffusion process becomes singularly perturbed, potentially resulting in collapsing effects. To effectively address this challenge and ensure numerical stability, common remedies include exploring a regularized version of the problem [START_REF] Droniou | Uniform-in-time convergence of numerical methods for non-linear degenerate parabolic equations[END_REF] or reformulating the system with a different form of nonlinearity [START_REF] Brenner | Improving newton's method performance by parametrization: The case of the richards equation[END_REF]. These approaches are employed, possibly in combination, to address degenerate cases, especially those featuring double degeneracy with potentially distinct rates, and effectively manage the resulting singular effects.

The purpose of numerical methods for Partial Differential Equations (PDEs) is to approximate the solution, as precisely as possible, to the underlying continuous solution. Accurately approximating complex PDEs can be challenging, especially when errors are intertwined, arising from various sources such as discretization, regularization, and linearization errors. To effectively address these challenges, the use of a posteriori error estimates becomes crucial. A key property of a reliable a posteriori estimate [START_REF] Ern | Eléments finis: théorie, applications, mise en oeuvre[END_REF][START_REF] Ern | A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF][START_REF] Ern | Adaptive inexact newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF][START_REF] Ern | Guaranteed, locally space-time efficient, and polynomialdegree robust a posteriori error estimates for high-order discretizations of parabolic problems[END_REF][START_REF] Ern | Guaranteed and robust discontinuous galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF] is its ability to quantify the error in the numerical approximation accurately, making it suitable for guiding adaptive algorithms. Adaptive algorithms enable efficient approximation of localized features in the solution even with intricate features. Considerable progress has been made in recent years regarding a posteriori analysis for linear elliptic equations [START_REF] Ern | Discrete p-robust H(div)-liftings and a posteriori estimates for elliptic problems with H -1 source terms[END_REF][START_REF] Smears | Simple and robust equilibrated flux a posteriori estimates for singularly perturbed reaction-diffusion problems[END_REF][START_REF] Cheddadi | Guaranteed and robust a posteriori error estimates for singularly perturbed reaction-diffusion problems[END_REF]. However, rigorous a posteriori error estimates for nonlinear elliptic problems, particularly in degenerate cases, remain a less-explored area [START_REF] Mitra | Guaranteed, locally efficient, and robust a posteriori estimates for nonlinear elliptic problems in iteration-dependent norms. An orthogonal decomposition result based on iterative linearization[END_REF][START_REF] Harnist | Robust energy a posteriori estimates for nonlinear elliptic problems[END_REF]. Explorations into degenerate parabolic equations have also been undertaken [START_REF] Cancès | Convergence and a posteriori error analysis for energystable finite element approximations of degenerate parabolic equations[END_REF][START_REF] Dolejší | A posteriori error estimate and mesh adaptation for the numerical solution of the richards equation[END_REF], typically under the assumption that nonlinearities are globally Lipschitz-continuous.

Ongoing efforts are dedicated to designing robust and efficient linearization schemes tailored for degenerate elliptic and parabolic problems [START_REF] Slodicka | A robust and efficient linearization scheme for doubly nonlinear and degenerate parabolic problems arising in flow in porous media[END_REF][START_REF] Bergamaschi | Mixed finite elements and newton-type linearizations for the solution of richard's equation[END_REF][START_REF] Hao | A gradient descent method for solving a system of nonlinear equations[END_REF][START_REF] Pop | Mixed finite elements for the richards' equation: linearization procedure[END_REF]. While Newton's method [START_REF] Droniou | Non-conforming Finite Elements on Polytopal Meshes[END_REF] is commonly the default choice due to its potential for quadratic convergence, it encounters challenges when the initial solution is not sufficiently close. In order to control degeneracy, Newton's method often involves the introduction of regularized versions of nonlinear functions. Another strategy is to explore globally convergent solvers such as the L-scheme, Picard, or damped Newton's methods; however, this comes at the trade-off of slower convergence. For instance, the L-scheme has demonstrated unconditional convergence for many model problems, but it's slow and only linearly convergent [START_REF] Mitra | A modified l-scheme to solve nonlinear diffusion problems[END_REF]. The quest for a balance between robustness and speed has spurred the development of hybrid solvers, as explored in [START_REF] Bergamaschi | Mixed finite elements and newton-type linearizations for the solution of richard's equation[END_REF][START_REF] Illiano | Iterative schemes for surfactant transport in porous media[END_REF]. In this context, multiple iterations of a slower scheme are employed to generate an initial approximation for Newton's method. This strategy was further refined in [START_REF] Stokke | An adaptive solution strategy for richards' equation[END_REF], where a switching criteria between the L-scheme and Newton's method for the Richards equation was guided by an a posteriori error estimator.

Addressing doubly degenerate elliptic problems presents an additional challenge for nonlinear solvers, demanding robustness in handling scenarios involving singular and double degeneracy [START_REF] Mitra | Guaranteed, locally efficient, and robust a posteriori estimates for nonlinear elliptic problems in iteration-dependent norms. An orthogonal decomposition result based on iterative linearization[END_REF]. Specifically, ϕ and β do not necessarily have to be Lipschitz-continuous or single-valued. This entails developing or combining solvers capable of adaptive switching and termination, all while selecting the appropriate regularization parameters and/or formulation [START_REF] Arioli | Stopping criteria for iterations in finite element methods[END_REF][START_REF] Di Pietro | Adaptive regularization, linearization, and discretization and a posteriori error control for the two-phase Stefan problem[END_REF][START_REF] Brenner | Improving newton's method performance by parametrization: The case of the richards equation[END_REF]. Figuring out when to stop or switch solvers depends on a key connection to the error in how the problem is discretized. A posteriori error analysis stands as a powerful tool for addressing these tasks [START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF][START_REF] Ahmed | A posteriori error estimates and stopping criteria for space-time domain decomposition for two-phase flow between different rock types[END_REF][START_REF] Ahmed | Adaptive asynchronous time-stepping, stopping criteria, and a posteriori error estimates for fixed-stress iterative schemes for coupled poromechanics problems[END_REF]. It not only guides the entire nonlinear and regularization process but also certifies errors [START_REF] Chamoin | A pedagogical review on a posteriori error estimation in finite element computations[END_REF]. A comprehensive review of controlling both linear and nonlinear solvers through adaptive stopping criteria can be found in [START_REF] Haberl | Convergence and quasi-optimal cost of adaptive algorithms for nonlinear operators including iterative linearization and algebraic solver[END_REF][START_REF] Févotte | Adaptive regularization, discretization, and linearization for nonsmooth problems based on primal-dual gap estimators[END_REF][START_REF] Becker | Cost-optimal adaptive iterative linearized fem for semilinear elliptic pdes[END_REF] and the references therein. When evaluating the effectiveness of a posteriori error estimators, key attributes to consider include reliability (upper bound for error) and efficiency (lower bound for error). Ideally, the constants in the upper bounds should be explicit, independent of PDE data, and finite-dimensional approximation parameters. The requirement of having a constant-free upper bound, coupled with a lower bound featuring constants ideally independent of model parameters, is particularly appealing for adaptively balancing error sources.

Building upon prior contributions, we extend our focus to adaptive solvers tailored for addressing nonlinear doubly degenerate elliptic equations within conforming methods. The proposed solver dynamically alternates between two linearization strategies: the linearly and globally L-scheme linearization and the quadratic and locally convergent Newton's method. Notably, in instances where Newton's method struggles with non-convergence, the L-scheme is enlisted to initiate the linearization process until an adequately precise starting point for Newton's method is attained. This decision is contingent on achieving a balance between spatial and nonlinear error components. Specifically, the L-scheme persists until the linearization error estimator drops below the spatial error estimator or becomes sufficiently small in comparison. Once this criterion is met, the solver transitions to the Newton method. Furthermore, the choice of the regularization parameter, final stopping criteria, and adaptive grid refinement/coarsening is informed by a posteriori indicators derived from local flux reconstructions and saturation post-processing.

This paper is structured as follows: In Section 2, we introduce essential notation, outline the finite element spaces employed, and describe the nonlinear degenerate elliptic model. Section 3 is focused on presenting the discrete linearized scheme. Therein, we discuss various linearization methods and introduce our adaptive solver, built through the formulation of a posteriori error estimators and the design of balancing criteria. Section 4 is dedicated to the construction of the a posteriori error estimators, leveraging two error measure based on the dual norm of the residual. We first introduce the residual and energy norms used to evaluate the distance between the approximate solution and the exact solution and demonstrate equivalence between the error and dual norm of the residual. Subsequently, we define the local a posteriori estimators, which rely on the calculated equilibrated flux and post-processed saturation. We also furnish proof of the guaranteed upper bound, local efficiency, and robustness (with respect to the strength of the nonlinearities) of our a posteriori estimate. Finally, in Section 5, we validate our adaptive solver through numerical experiments involving four illustrative examples.

Model problem

Sobolev spaces and notations

Let Ω ⊂ R d , with d ∈ {2, 3}, represent an open, bounded, and connected domain. We define H 1 (Ω) as the standard Sobolev space of L 2 (Ω) functions with weak gradients in L 2 (Ω)

d , and V := H 1 0 (Ω) as its subspace consisting of functions with zero trace. The space H(div) denotes the set of L 2 (Ω)

d functions with weak divergences in L 2 (Ω). The symbols ∇ and ∇• (or div) respectively represent the weak gradient and divergence operators in Ω. For any subdomain ω of Ω, we use (•, •) ω to denote the L 2 -inner product, || • || ω to indicate the corresponding norm (with the index omitted when ω := Ω), and |ω| to represent the Lebesgue measure on ω. Furthermore, we define H -1 (ω) as the topological dual space of H 1 0 (ω), and we express the duality bracket between these spaces as v, w -1,1,ω , which reduces to the L 2 -inner product if v, w ∈ L 2 (ω). When ω = Ω, the duality bracket is denoted simply as . , . -1,1 .

Partition of Ω and finite element spaces

Let T h be a family of closed simplices, such that Ω = K∈T h K. The partition is assumed to be conforming, meaning that if K and L are in T h , with K = L, then K ∩ L is either an empty set, a common face, edge, or vertex of both K and L. We denote the set of all faces of T h by F h , which is composed of F int h ∪ F ext h , where F int h and F ext h represent the set of internal and external (boundary) faces, respectively. The faces of an element K are collected in the set F K . The diameter of an element K is denoted by h K , and the diameter of a face σ ∈ F h is denoted by h σ . We set h = max K∈T h h K . The set of vertices of the mesh T h is denoted by V h , with V int h representing the set of interior vertices and V ext h for the set of boundary vertices. For each a ∈ V h , we use T a h to denote the patch of the vertex a, i.e., all the elements K ∈ T h that share the vertex a. We denote by ω a the corresponding open subset of Ω with diameter h ωa . Additionally, we define F a ⊂ F h as the set of all faces in T a h , with F int a and F ext a representing its interior and boundary faces, respectively. we define F ∂Ω a = F a F ext h . For each element K ∈ T h , we collect in V K the set of vertices belonging to K from V h . Furthermore, for an element K ∈ T h , and for any edge σ ∈ ∂K, we define n σ as the unit normal vector on σ pointing outside K.

Let P p (K) represent the space of polynomials on K ⊂ T h of total degree less than or equal to p. The set P p (T h ) denotes the corresponding space of piecewise p-degree polynomials on T h , meaning

P p (T h ) = {p h ∈ L 2 (Ω); p h | K ∈ P p (K), ∀K ∈ T h }. Additionally, we consider the piecewise Raviart-Thomas-Nédélec space RTN p (T h ) ⊂ L 2 (Ω), defined by RTNp(T h ) = {v h ∈ H(div); v h | K ∈ RTN p (K), ∀K ∈ T h }, where RTN p (K) = [P p (K)] d + P p (K)x. Let Π h : L 2 (Ω) → P p (T h ) denote the L 2 -orthogonal projection operator onto P p (T h ) or (P p (T h ; R d )). We will also need Π RTN h : L 2 (Ω) → RTN p (T h ) denote the L 2 -orthogonal projection operator onto RTN p (T h ).
For a vertex a ∈ V h , we introduce the hat operator ψ a ∈ P 1 (T h )∩H 1 (Ω) such that ψ a (a) = 1 and is zero otherwise (at other vertices), thus, supp(ψ a ) = T a h . For each a ∈ V h , RTN p (T a h ) and P p (T a h ) denote the restriction of RTN p (T h ) and P p (T h ) to the patch T a h , respectively. We also introduce V h := P p (T h ) ∩ H 1 0 (Ω) and let Q h := P p (T h ) and W h := RTN p (T h ).

Model problem and variational formulation

We consider in this paper the following nonlinear and doubly degenerate elliptic equation as the basis of our discussion [START_REF] Carrillo | Entropy solutions for nonlinear degenerate problems[END_REF][START_REF] Droniou | High-order mass-lumped schemes for nonlinear degenerate elliptic equations[END_REF]: find u such that

β(u) -∆ϕ(u) = f + ∇•ξ, in Ω, (2.1a) u(•, 0) = u 0 , in Ω, (2.1b) ϕ(u) = 0, on ∂Ω. (2.1c)
As a preliminary step towards establishing the weak formulation of (2.1), we address the necessary assumptions that must hold: (A1) ϕ : R → R is continuous and non-decreasing, further ϕ(0) = 0, and there exist M 0 , M 1 > 0 such that |ϕ(s)| ≥ M 0 |s| -M 1 , for all s ∈ R .

(A2) β : R → R is continuous non-decreasing, and β(0) = 0, and there exist

K 0 , K 1 > 0 such that |β(s)| ≤ K 0 |s| + K 1 , for all s ∈ R. (A3) ϕ + β : R → R is strictly increasing. (A4) The source terms are such that f ∈ L 2 (Ω) and ξ ∈ (L 2 (Ω)) d .
The superlinearity of ϕ(u) in (A1) ensures that u belongs to L 2 (Ω) (at least). Thus, the sub-linearity of β assumed by (A2) implies that β(u) ∈ L 2 (Ω) as well. That, assumption (A4) implies that the resulting source term f + ∇•ξ should be understood in the sense of distribution [START_REF] Cohen | Convergence rates of AFEM with h -1 data[END_REF][START_REF] Ern | Discrete p-robust H(div)-liftings and a posteriori estimates for elliptic problems with H -1 source terms[END_REF]. Note that further generalizations are possible, bringing more technicalities. Assume (A1)-(A4) hold true, we introduce the weak solution to problem (2.1):

find u ∈ L 2 (Ω) with ϕ(u) ∈ V such that (β(u), v) + (∇ϕ(u), ∇v) = (f, v) -(ξ, ∇v) , ∀v ∈ V. (2.2)
The theoretical aspects of (2.2) can be found in the seminal paper [START_REF] Carrillo | Entropy solutions for nonlinear degenerate problems[END_REF]. The existence and uniqueness of the weak solution to (2.2) are established in [START_REF] Droniou | High-order mass-lumped schemes for nonlinear degenerate elliptic equations[END_REF], where detailed regularity properties are also presented (see also [START_REF] Droniou | The gradient discretization method for slow and fast diffusion porous media equations[END_REF]). Still following [START_REF] Droniou | The gradient discretisation method[END_REF], one can opt for a reformulation of 2.2 using a new variable, say w, to serve as the primary unknown in the new system, in which Lipschitz-continuous counterparts of ϕ and β are obtained. Under Assumption (A1)-(A4), we observe that β + ϕ : R → R is a bijective function. This allows us to define µ(s

) = ζ (β + ζ) -1 (s) and ρ(t) = t -µ(t) = β (β + ζ) -1 (t)
. These functions are non-decreasing and have 1-Lipschitz continuity (see [START_REF] Droniou | High-order mass-lumped schemes for nonlinear degenerate elliptic equations[END_REF]). Problem 2.2 can be equivalently expressed as

finding w ∈ L 2 (Ω) with ρ(w) ∈ V such that (µ(w), v) + (∇ρ(w), ∇v) = (f, v) -(ξ, ∇v) , ∀v ∈ V. (2.3)
This reformulation brings advantages for a priori error analysis in the context of conforming methods. This is mainly because it provides more flexibility in using the chain rule and choosing suitable test functions compared to 2.2. Specifically, the Lipshitz continuity of µ and ρ makes formulation 2.2 more suitable for conforming methods. However, it's worth noting that computing the inverse transform w -→ u can be computationally demanding, and w doesn't have a direct physical interpretation. The equivalence between the formulations comes from the uniqueness of solutions, making both applicable for our a posteriori error analysis. For a detailed exploration of the strengths and weaknesses in these two formulations, along with the numerical approaches applied to (2.1), we refer the reader to [START_REF] Droniou | High-order mass-lumped schemes for nonlinear degenerate elliptic equations[END_REF]. This work introduces the Gradient Discretisation Method (GDM), a generic numerical analysis framework enabling a unified convergence analysis of different discretisations [START_REF] Droniou | The gradient discretisation method[END_REF].

Discrete setting and adaptive algorithm

In this section, we provide the discretized, linearized, and regularized version of the variational formulation (2.2). We also detail the methodology and algorithm utilized to solve this formulation.

Discretization

To illustrate our framework, we present a conforming Galerkin method applied to (2.2). We aim to find u h ∈ V h satisfying the following discrete formulation:

(β(u h ), v) + (∇ϕ(u h ), ∇v) = (f, v) -(ξ, ∇v) , ∀v ∈ V h . (3.1)
Yet, due to the necessity of mass-lumping to ensure a monotone scheme [START_REF] Eymard | Gradient schemes for the stefan problem[END_REF][START_REF] Nochetto | An adaptive finite element method for two-phase stefan problems in two space dimensions. i. stability and error estimates[END_REF], there is limited interest in employing high-order methods for 3.1. However, one can still rely on the Gradient Discretisation Method (GDM) framework corresponding to MLFE to ensure higher-order convergence, as detailed in [START_REF] Droniou | Non-conforming Finite Elements on Polytopal Meshes[END_REF]. Furthermore, to derive any form of a posteriori error estimate from these conforming schemes, we are constrained to assume ϕ(u h ) ∈ V h , a condition that cannot be guaranteed by the above formulation due to the lack of global Lipschitz continuity in ϕ. To address this, we may use the discrete counterpart of (2.3), which provide the necessary V -conformity for ρ(w h ) and, by equivalence, for ϕ(u h ). We shall later comment on the implications of each formulation on the a posteriori error analysis. Nevertheless, we maintain the generality of our presentation and analysis. Our framework covers conforming schemes where the discrete solution

u h ∈ L 2 (Ω) is such that ϕ(u h ) ∈ V .
Examples of methods fitting into our analysis include conforming Finite Element methods (FE), including the Mass-Lumped version (MLFE) and Vertex Approximate Gradient (VAG) methods [START_REF] Droniou | The gradient discretisation method[END_REF].

Linearization

The primary challenge in addressing the non-linear and doubly degenerate problem (3.1) lies in devising a robust linearization scheme that converges, particularly in cases where ϕ (•) and/or β (•) approach(es) zero. Equally significant is the scenario where ϕ (•) and/or β (•) become(es) unbounded. This is of important interest due to its potential impact on the convergence properties of the solution method. Let us assume a linearization process is adopted to deal with nonlinearities within (3.1): At the iteration k ≥ 1, find u k h ∈ V h such that

β k-1 (u k h ), v + ∇ϕ k-1 (u k h ), ∇v = (f, v) -(ξ, ∇v) , ∀v ∈ V h , (3.2)
where the functions ϕ k-1 and β k-1 are respectively linear approximations of ϕ(u h ) and β(u h ). If Newton-Raphson's method is used, these approximations are of the form:

ζ k-1 (u k h ) := ζ (u k-1 h ) u k h -u k-1 h + ζ(u k-1 h ), k ≥ 1. (3.3a)
Although, Newton-Raphson's method is quadratically convergent, it is only locally convergent and involves the computation of derivatives; when ζ(•) degenerates, the method might fail to converge (algebraic residuals are either stagnating, diverging or oscillating) and a remedy based on relaxation techniques is needed. We might also need smoothing functions of ϕ(•) and β(•).

An alternative to Newton-Raphson's method is to consider a quasi-Newton version that replaces the exact derivative ζ with an approximation. For instance, we adopt the so-called L-scheme, a quasi-Newton method which lies between fixed-point and Newton solvers [START_REF] Pop | Mixed finite elements for the richards' equation: linearization procedure[END_REF][START_REF] Slodicka | A robust and efficient linearization scheme for doubly nonlinear and degenerate parabolic problems arising in flow in porous media[END_REF][START_REF] Yong | A numerical approach to porous medium equations[END_REF], in which the approximations in (3.2) are simply given by

ζ k-1 (u k h ) := L u k h -u k-1 h + ζ(u k-1 h ), k ≥ 1, (3.3b)
with L > 0 is a constant mimicking the derivative of ζ. The L-scheme is linearly convergent assuming Lipschitz continuity of ζ(•). The only requirement for the convergence of the L-scheme is that the parameter L should satisfy L ≥ L ζ := sup s |ζ (s)|. It has been shown that the convergence can also be obtained for smaller values of L, e.g. for L = L ζ /2 convergence is guaranteed (see [START_REF] Mitra | A modified l-scheme to solve nonlinear diffusion problems[END_REF]) while the contraction proof can not be ensured in this case. Also, one can choose L = max s |ζ (s)|) and that's the classical Picard iterations (see [START_REF] Illiano | Iterative schemes for surfactant transport in porous media[END_REF] for details). Therefore, the L-scheme is robust when ζ (•) ≈ 0, however, when ζ (•) becomes unbounded (ζ is not Lipschitz), a regularization procedure is also needed. We refer to [START_REF] Yang | Nature-inspired optimization algorithms[END_REF] for general sub-classes of quasi-Newton methods and particularly to [START_REF] List | A study on iterative methods for solving richards' equation[END_REF][START_REF] Bergamaschi | Mixed finite elements and newton-type linearizations for the solution of richard's equation[END_REF][START_REF] Hao | A gradient descent method for solving a system of nonlinear equations[END_REF] for combinations of Newton and L-scheme or Picard iterations.

Regularization

The regularization approach involves introducing a smoothing function, denoted as ζ , for a parameter > 0.

This smoothing function ζ (•) : R → R is monotone and Lipschitz continuous,

0 ≤ |ζ (s)| ≤ L ζ < ∞, ∀s ∈ R. (3.4)
Let u h be the regularized solution, defined by the recursive sequence u ,k h ∈ V h , ∀k ≥ 1, such that

β k-1 (u ,k h ), v + ∇ϕ k-1 (u ,k h ), ∇v = (f, v) -(ξ, ∇v) , ∀v ∈ V h . (3.5)
Although this regularization procedure could prove effective for an L-scheme, strict monotonicity might be a necessary condition for the application of Newton's method. Furthermore, if the L-scheme is employed for the alternative formulation 2.3, regularization may not be necessary at all. It is evident that developing an efficient solver that maintains a balance between robustness and efficiency while addressing all nonlinear aspects of the problem necessitates combining nonlinear solvers and/or formulations. This can be achieved through the so-called adaptive solver/variable switching [START_REF] Stokke | An adaptive solution strategy for richards' equation[END_REF] or parameterization [START_REF] Brenner | Improving newton's method performance by parametrization: The case of the richards equation[END_REF], enabling us to leverage the strengths of each method or formulation. In our context, we will focus exclusively on adaptive solver switching to optimize the performance of the standard nonlinear solver tailored for our specific problem.

Adaptive algorithm

The purpose of this section is to reduce as much as possible the computational effort to solve (3.5). The estimation and localization of distinct error components lead to the development of an adaptive nonlinear solver, enabling optimal balancing of different error contributions. In doing so, one in particular can supervise Newton iterations with an adaptive switch back-and-forth to a slower but safer algorithm, if a convergence failure is detected.

Balancing criteria

Let η k, disc , η k, reg , η k, qd , η k, osc and η k, lin , be respectively fully computable a posteriori estimators of discretization error, regularization error, quadrature error, data-oscillation error, and linearization error at the k-th iteration, such that

error(u, u ,k h ) ≤ η k, disc + η k, reg + η k, qd + η k, osc + η k, lin , (3.6)
where error(u, u ,k h ) represents an energy norm, and its definition will be provided in a posteriori error analysis in Section 4.1. To minimize efficiently the right-hand side of (3.6), we adopt suitable strategy for each estimator so that each of will looks forward the smallest possible cost to himself, and such that the ideal estimate is obtained. The strategy we adopt for η k, lin is to use an adaptive switching from the L-scheme to Newton-Raphson's method, i.e., we switch from the L-scheme to Newton iteration when the following criteria is satisfied

η k, lin ≤ Γ sw η k, disc + η k, reg + η k, qd + η k, osc , (3.7)
where Γ sw < 1 is a user-given weight to be chosen adaptively. The use of (3.7) implies that the L-scheme is pursued until the L-scheme error approaches enough from the total error, which hopefully is sufficiently good (but perhaps not!) to (re)-launch the Newton solver. This later is then continued until the following stopping criteria

η k, lin ≤ Γ lin η k, disc + η k, reg + η k, qd + η k, osc , (3.8)
is thrown. The parameter Γ lin is to be given in (0, Γ sw ). If the L-scheme fails to deliver a good initial guess for Newton, the L-scheme will be pursued either with smaller parameters Γ sw or until the stopping criteria (3.8) is satisfied. The strategy we adopt for minimizing η k, reg is based on reducing adaptively the regularization parameter so that

η k, reg ≤ Γ reg η k, disc + η k, qd + η k, osc , (3.9)
with a user-given parameter Γ reg ∈ (0, 1). The strategy we adopt for η k, disc is to adapt the space mesh; refinement/coarsening so that the local discretization error estimators are distributed equally: for all

K 1 , K 2 ∈ T h , η k, disc,K1 ≈ η k, disc,K2 . (3.10)
To address restrictions in computing resources, we establish a fixed refinement threshold h min > 0 for the mesh size. This involves ensuring that min

K∈T h h K ≥ h min .
We also assume that the adaptive mesh is shape regular, i.e., that there exists a constant κ T such that, for all K ∈ T h , max

K∈T h h K ρ K ≤ κ T ,
where h K is the diameter of K, and ρ K is the radius of the largest inscribed ball of K.

Algorithm

Recall that u ,k h is an approximation of the solution u obtained after k linearization iterations using a regularization parameter . Our algorithm, incorporating the aforementioned balancing criteria, is as follows:

Algorithm 3.1: The adaptive algorithm (1) Choose an initial approximation u ,-1 h ∈ V h of u h , switching parameters (Γ sw , Γ lin , Γ reg ), and a regularization parameter > 0. Set k = -1.
(2) Do {Spatial adaptivity}

(I) Do {Regularization} (A) Set k = -1. (B) Do {Newton-Raphson} (i) Increase k = k + 1. (ii) Compute u ,k h using (3.3a) and (3.5). (iii) Compute the estimators (η k, disc , η k, lin , η k, reg , η k, qd , η k, osc ). (iv) if η k, lin increases, stagnates, or oscillates. /*Switching solver*/ (a) Reset k = -1. (b) Do {Quasi-Newton} (1) Increase k = k + 1.
(2) Compute u ,k h using (3.3b) and (3.5). (3) Compute the estimators

(η k, disc , η k, lin , η k, reg , η k, qd , η k, osc ). While (3.7) or (3.8) is not satisfied. End While (3.8) is not satisfied. /*Stopping solver*/ (C) if (3.9) is not satisfied. (i) Replace by /2. End While (3.9) is not satisfied. (II) Refine or coarsen the cells K ∈ T h until (3.10) is satisfied or h K < h min .
While (3.10) is not satisfied.

Remark 3.1 (Offline L calculation). Given that the classical stabilization parameters, such as L ϕ or

Lϕ 2 , may not always yield optimal results, an enhancement to the adaptive algorithm 3.1 can be achieved through a 'brute optimization' process (see [START_REF] Storvik | On the optimization of the fixed-stress splitting for biot's equations[END_REF]). This process aims to determine the optimal parameters (L opt β , L opt ϕ ) by solving the optimization problem:

(L opt β , L opt ϕ ) := arg min a≤L ζ ≤b,ζ∈{ϕ,β} η k, lin (L β , L ϕ ; T H , ) (3.11)
This optimization is performed over a coarse grid T H and for a sufficiently small regularization parameter . The advantages lie in the discretization-error-independent nature of η k, lin and the almost regularizationindependent performance of the L-scheme. Remark 3.2 (Localized balancing criteria). We can express localized versions of the criteria on an elementby-element basis using local estimators, ensuring their validity for all K ∈ T h . The switching criteria (3.7) can also be tailored to balance linearization and spatial errors:

γ sw • η k, • ≤ η k, lin ≤ Γ sw • η k, • , (3.12)
where Γ sw and γ sw are set close to one. Balancing linearization and spatial errors, especially with a small η k, reg through (3.9), prevents specific errors from dominating. This delivers a smoother starting point for Newton, avoiding oscillations and promoting stability.

Remark 3.3 (Adaptive switching).

The switching parameter Γ sw is to be chosen smaller if quasi-Newton's method fails to provide a good initialization for Newton's method.

Remark 3.4 (Descent methods). The adaptive algorithm is valid if we replace the L-scheme with damped

Newton's methods or any other descent solver with explicit formula of the step size [START_REF] Hao | A gradient descent method for solving a system of nonlinear equations[END_REF]. The damping parameter can be computed in the offline phase. For optimal descent methods, we must constitute a line search rule to infer the ideal step size at each iteration.

Energy-type a posteriori error estimates

In the present work, we adopt energy-type a posteriori error estimates [START_REF] Ahmed | A posteriori error estimates and stopping criteria for space-time domain decomposition for two-phase flow between different rock types[END_REF][START_REF] Smears | Simple and robust equilibrated flux a posteriori estimates for singularly perturbed reaction-diffusion problems[END_REF][START_REF] Ern | Guaranteed, locally space-time efficient, and polynomialdegree robust a posteriori error estimates for high-order discretizations of parabolic problems[END_REF]. Classically, the derivation of energy a posteriori error estimates is based on the definition the residual (and its dual norm) and the mathematical link between the energy norm of the error and the dual norm of the residual.

Residuals and energy norms

Following [START_REF] Ern | Guaranteed, locally space-time efficient, and polynomialdegree robust a posteriori error estimates for high-order discretizations of parabolic problems[END_REF], we introduce the residual operator R u (u h ) ∈ H -1 (Ω), defined by

R u (u h ), v := (β(u) -β(u h ), v) + ∇ ϕ(u) -ϕ(u h ) , ∇v , ∀v ∈ V. (4.1)
Next, we let the parameter δ be equal to 0 or 1, depending on the norm in which the error is measured. We then equip the V space with the energy norm

ψ := ∇ψ 2 + δ ψ 2 1 2 , (4.2)
where clearly the natural choice being ψ := ∇ψ with δ = 0. The dual norm of the residual is then given through

|||R u (u h )||| := sup ψ∈H 1 0 (Ω) ψ =1 R u (u h ), ψ , (4.3)
which expresses by how much the trial function u h fails to satisfy (2.2). We use simply |||R(u h )||| for δ = 0. Note that, in the case of ϕ(u) := u and β(u

) := u + g, with given g ∈ L 2 (Ω), that is, when (2.
2) is linear and strongly elliptic, the norm • with δ = 1 is a natural energy norm on H 1 0 (Ω). Remark also that the chosen H 1 0 (Ω)-norm delivers different continuity estimates for the equation (2.2). The sharpness of the continuity estimate is in fact of particular interest for designing goal-oriented error estimation. For energyerror estimation, the choice of the H 1 0 (Ω)-norm affects the dual norm of the residual and leads to different error measures. Evidently, the dual norm of the residual |||R u (u h )||| vanishes for u h = u (see Corollary (4.4) for details), and this is when δ = 0 or 1, but its actual size is unknown and its value is incomputable. Thus, the link between these dual norms and the underlying error energy norms is necessary. To show this link, we introduce the energy norm

||u -u h || 2 # = ||β(u) -β(u h )|| 2 H -1 (Ω) + ||∇ (ϕ(u) -ϕ(u h )) || 2 . (4.4)
We also introduce the weaker distance

J u (u h ) := 2 β(u) -β(u h ), ϕ(u) -ϕ(u h ) . (4.5)
We have this result.

Proposition 4.1 (Residual-Energy norm). Let u be the weak formulation of (2.2) with V equipped with the energy-norm (4.2) and let u h ∈ L 2 (Ω) be an arbitrary function such that ϕ(u h ) ∈ V . Then, if δ = 0, there holds

J u (u h ) + ||u -u h || 2 # = |||R u (u h )||| 2 . (4.6) For δ = 1, define the metric ||u -u h || 2 Σ, := ||u -u h || 2 # + ||ϕ(u) -ϕ(u h )|| 2 H -1 (Ω) -||ϕ(u) -ϕ(u h )|| 2 . (4.7)
Then, there holds,

||u -u h || 2 Σ, = |||R u (u h )||| 2 . (4.8)
Proof. To obtain the above relations, we use duality techniques as in [START_REF] Ern | Guaranteed, locally space-time efficient, and polynomialdegree robust a posteriori error estimates for high-order discretizations of parabolic problems[END_REF][START_REF] Brun | Modeling the process of speciation using a multiscale framework including a posteriori error estimates[END_REF]. Utilizing Assumption (A2), which ensures β(u) and β(u h ) are in L 2 (Ω), we proceed by letting ψ ∈ V be the solution to the auxiliary problem:

(∇ψ, ∇v) = (β(u) -β(u h ), v) , ∀v ∈ V. (4.9)
Clearly, we have

||∇ψ|| = ||β(u) -β(u h )|| H -1 (Ω) . (4.10)
We replace (4.9) in (4.1) to get

R u (u h ), v = (∇ (ψ + ϕ(u) -ϕ(u h )) , ∇v) , ∀v ∈ V. (4.11) so that from (4.3), |||R u (u h )||| = ||∇ (ψ + ϕ(u) -ϕ(u h )) ||. (4.12)
By expanding the square, we obtain

|||R u (u h )||| 2 = ||∇ψ|| 2 + ||∇ (ϕ(u) -ϕ(u h )) || 2 + 2(∇ψ, ∇ (ϕ(u) -ϕ(u h ))). (4.13)
We replace the first term using (4.10) as well as we make use of (4.9) to get (4.6).

For the second identity (4.8), we let ψ ∈ V be the solution of the auxiliary problem

(∇ψ, ∇v) + (ψ, v) = (β(u) -β(u h ) -ϕ(u) + ϕ(u h ), v) , ∀v ∈ V. (4.14)
We have first

ψ = ||β(u) -β(u h ) -ϕ(u) + ϕ(u h )|| H -1 (Ω) . (4.15)
Then it holds

R u (u h ), v = (ψ + ϕ(u) -ϕ(u h ), v) + (∇ (ψ + ϕ(u) -ϕ(u h )) , ∇v) , ∀v ∈ V.
Thus, also

|||R u (u h )||| 2 = ψ + ϕ(u) -ϕ(u h ) 2 .
Expanding the terms, utilizing the definitions of ψ and d Σ , we finally obtain

|||R u (u h )||| 2 = ψ 2 + ϕ(u) -ϕ(u h ) 2 + 2 ψ, ϕ(u) -ϕ(u h ) + 2 ∇ψ, ∇ (ϕ(u) -ϕ(u h )) = ψ 2 + ϕ(u) -ϕ(u h ) 2 + 2 β(u) -β(u h ), ϕ(u) -ϕ(u h ) -2 ϕ(u) -ϕ(u h ) 2 = ||u -u h || 2 Σ, .
This completes the proof.

Remark 4.2 (About J u (u h )). Note that J u (u h ) is positive due to the monotonicity of β and ϕ. As Arbogast [START_REF] Arbogast | A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media[END_REF] points out, the nonlinear form J u (u h ) tells us something about the error of the scheme [START_REF] Woodward | Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media[END_REF][START_REF] Brun | Modeling the process of speciation using a multiscale framework including a posteriori error estimates[END_REF]. Particularly, if β = I and ϕ is Lipschitz-continuous with Lipschitz constant L ϕ , then,

1 L ϕ ||ϕ(u) -ϕ(u h )|| 2 ≤ J u (u h ). (4.16)
Remark 4.3 (On norm(u, u h )). Note that in the linear case, as indicated by [START_REF] Smears | Simple and robust equilibrated flux a posteriori estimates for singularly perturbed reaction-diffusion problems[END_REF], we have:

||u -u h || 2 + ||∇ (u -u h ) || 2 = |||R u (u h )||| 2 , ∀u, u h ∈ V. (4.17)
For the general nonlinear case, there holds

||u -u h || 2 Σ, ≤ ||u -u h || 2 # ≤ ||β(u) -β(u h )|| 2 + ||∇ (ϕ(u) -ϕ(u h )) || 2 , ∀u, u h ∈ V. (4.18)
Corollary 4.4 (On the weak solution). There exists one solution to problem (2.2). Furthermore, there holds

δ ||u|| Σ, + (1 -δ ) J 0 (u) + ||u|| 2 # 1 2 = sup v∈H 1 0 (Ω) v =1 {(f, v) -(ξ, ∇v)}. (4.19)
Proof. The existence of a solution is a consequence of [24, Theorem 2.9]. To show the uniqueness of the solution under Assumption (A1)-(A4), we will see that if R(u h ) = 0, we then have u h = u. Suppose that R(u h ) = 0, then from (4.6) we have immediately that ϕ(u) -ϕ(u h ) = 0 and β(u) -β(u h ) = 0. The use of Assumption (A3) implies that u h = u and the uniqueness of the solution. To infer the uniqueness of the solution from the weaker identity (4.8), we consider u, u h as two solutions of problem (2.2). We then subtract their respective equations and let

v = ϕ(u) -ϕ(u h ), to get J u (u h )/2 + ||∇(ϕ(u) -ϕ(u h ))|| 2 = 0.
Due to the positivity of J u (u h ), it follows that ϕ(u)-ϕ(u h ) = 0. The weaker error measure (4.8) also results in β(u) -β(u h ) = 0. Again, the application of Assumption (A3) implies that u h = u. Finally, obtaining the identity 4.19 is straightforward, given that ϕ(0) = β(u) = 0 (and zero is an admissible element of V ). We achieve this by substituting (4.3) into (4.6) and (4.8), respectively.

Equilibrated flux σ ,k h and post-processed saturation Φ ,k h

The constructions introduced in this section follows closely [START_REF] Ern | Adaptive inexact newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF], in generalization of [START_REF] Smears | Simple and robust equilibrated flux a posteriori estimates for singularly perturbed reaction-diffusion problems[END_REF]. Let us introduce for each a ∈ T a h , the local spatial mixed finite element spaces posed on the patch domain ω a :

Q a h := P p (T a h ) if δ > 0 or a ∈ V ext h , {q ∈ P p (T a h ), (q, 1) ωa = 0} if δ = 0 and a ∈ V int h , (4.20a)
W a h := {v h ∈ RTN p (T a h ) ∩ H(div, ω a ), v h • n K = 0 on ∂ω a \ ∂Ω} if a ∈ V ext h , {v h ∈ RTN p (T a h ) ∩ H(div, ω a ), v h • n K = 0 on ∂ω a } if a ∈ V int h .
(4.20b) Definition 4.5 (Notations). We denote by I ,k h the linearized approximate flux ∇ ϕ k-1 (u ,k h ) within (3.5). The constructions given next are based on independent local mixed finite element approximations of residual problems from (3.5) over patches of elements around mesh vertices (See Section (3.2) for notations). As used in the above discrete spaces, the reconstructions will depend on the choice of the error measure. To simplify the notation, we introduce the residual function r

,k h ∈ P p (T h ), r ,k h | ∂Ω = 0, represented by (r ,k h , v) := (f, v) -( β (u ,k h ), v) -((ξ + ∇ ϕ (u ,k h )), ∇v), ∀v ∈ V h . (4.21)
Then, for a given parameter δ = 1 or 0, the equilibrated flux σ ,k h and postprocessed saturation Φ ,k h are obtained using the following reconstructions (See Appendix A. for the local minimization problems):

Reconstruction 1 ((σ ,k a , Φ ,k a )). At each iteration k ≥ 1, and for each patch ω a , a ∈ V h , we construct (σ ,k a , r ,k a ) := (σ ,k a , r ,k a ) disc + (σ ,k a , r ,k a ) lin ∈ W a h × Q a h such that: (σ ,k a,disc + σ ,k a,lin , v) ωa -(∇•v, r ,k a ) ωa = (λ a , v) ωa , ∀v ∈ W a h , (4.22a)
(∇•(σ ,k a,disc + σ ,k a,lin ), q) ωa + δ (r ,k a , q) ωa = (g a , q) ωa , ∀q ∈ Q a h , (4.22b)
with

g a := Π h ψ a f -β k-1 (u ,k h ) | ωa -∇ψ a • (Π h ξ + I ,k h )| ωa , and λ a := -(ψ a ξ + ψ a I ,k h )| ωa .
Then, we define the postprocessed saturation by

Φ ,k a = Φ ,k a,disc + Φ ,k a,lin := δ r ,k a + Π h (ψ a β k-1 (u ,k h )) and we extend each (σ ,k a , Φ ,k a ) by zero in Ω \ ω a , (σ ,k h , Φ ,k h ) ∈ RTN p (T h ) × P p (T h ) is given by (σ ,k h , Φ ,k h ) := a∈V h (σ ,k a , Φ ,k a ). (4.23) Reconstruction 2 ((σ ,k a , Φ ,k a ) disc ). At each iteration k ≥ 1, and for each patch ω a , a ∈ V h , construct (σ ,k a , r ,k a ) disc ∈ W a h × Q a h such that: (σ ,k a,disc , v) ωa -(∇•v, r ,k a,disc ) ωa = (λ ,k a , v) ωa , ∀v ∈ W a h , (4.24a)
(∇•σ ,k a,disc , q) ωa + δ (r ,k a,disc , q) ωa = (g ,k a -ψ a r ,k h , q) ωa , ∀q ∈ Q a h , (4.24b)
with

g ,k a := Π h ψ a f -β (u ,k h ) | ωa -∇ψ a • Π h ξ + ∇ ϕ (u ,k h ) | ωa , λ ,k a := -ψ a ξ + ψ a ∇ ϕ (u ,k h ) | ωa .
Then, we define the postprocessed saturation by Φ ,k a,disc := δ r ,k a,disc + Π h (ψ a β (u ,k h )) and we extend each

(σ ,k a , Φ ,k a ) disc by zero in Ω \ ω a , (σ ,k h , Φ ,k h ) disc ∈ RTN p (T h ) × P p (T h ) is given by (σ ,k h , Φ ,k h ) disc := a∈V h (σ ,k a , Φ ,k a ) disc . (4.25)
The key properties of the above reconstructions with δ ∈ {0, 1} are summarized in the following result.

Lemma 4.6 (Properties of (σ ,k h , Φ ,k h )). Let (σ ,k h , Φ ,k h ) be given by Reconstruction 1 or 2. Then it hold

1. (σ ,k h , Φ ,k h ) ∈ H(div, Ω) × L 2 (Ω), 2. ∇•σ ,k h + Φ ,k h = Π h f , with 2.a. ∇•σ ,k h,disc + Φ ,k h,disc = Π h f -r ,k h , 2.b. ∇•σ ,k h,lin + Φ ,k h,lin = r ,k h , 3. (∇•σ ,k a,disc , 1) ωa = (∇•σ ,k a,lin , 1) ωa = 0, ∀a ∈ V int h , 4. (Φ ,k a,disc , 1) ωa = (ψ a β (u ,k h ), 1) ωa and (Φ ,k a,lin , 1) ωa = (ψ a ( β k-1 (u ,k h ) -β (u ,k h )), 1) ωa , ∀a ∈ V int h .
Proof. Its straightforward to see that for any a ∈ V h the zero extension of σ ,k a is H(div,

, ∇•σ ,k h,disc + Φ ,k h,disc = a∈V h ∇•σ ,k a,disc + Φ ,k a,disc = a∈V h ∇•σ ,k a,disc + δ r ,k a,disc + Π h (ψ a β (u ,k h )) = a∈V h {Π h (ψ a f - ψ a β (u ,k h ) + ψ a β (u ,k h )) -ψ a r ,k h -∇ψ a • Π h ξ + ∇ ϕ (u ,k h ) } = Π h f -r ,k h . To show 3, we have for a ∈ V int h , for (∇•σ ,k a, , 1) ωa = (σ ,k a,
• n ωa , 1) ∂ωa = 0 because of the zero Neumann boundary condition in W a h . The property 4 is trivial when δ = 0, and for δ = 1, we can remark that the source term g ,k a in (4.24b) satisfies

(g ,k a -ψ a r ,k h , 1) ωa = 0, ∀a ∈ V int h . (4.26)
As a result, we can use (4.24b) and property 3 to see that (r ,k a,disc , 1) ωa = 0 as well, and thus (Φ ,k a,disc , 1) ωa = (ψ a β (u ,k h ), 1) ωa . We repeat the same steps to get that (Φ ,k a,lin , 1

) ωa = ψ a β k-1 (u ,k h )-β (u ,k h ) , 1 ωa .
Remark 4.7 (Cheaper flux and saturation reconstruction). In accordance with the approach outlined in [START_REF] Cancès | Guaranteed and robust a posteriori bounds for laplace eigenvalues and eigenvectors: Conforming approximations[END_REF], it is possible to perform a reconstruction with one degree lower in the space RTN p-1 (T h ) × P p-1 (T h ).

While this method can be satisfactory in practical applications, it is important to note that the theoretical polynomial degree robustness (as expressed in Eq. (4.49)) may be compromised.

Remark 4.8 (On the (ρ -µ) formulation). The above reconstruction and the analysis linking the dual norm of the residual with the underlying error energy norms in Proposition 4.1 can be adapted for the (ρ -µ) formulation. This adjustment simply requires replacing the nonlinear functions with their counterparts from (2.3).

Guaranteed reliability

We recall here two basic inequalities necessary in order to obtain our a posteriori error estimates. Let K be an element of T h , we recall the Poincaré-Friedrichs inequality,

||v -Π h v|| K ≤ C P,K ||∇v|| K , for any v ∈ H 1 (K). (4.27)
Since simplices are convex, there holds C P,K = h K π . The generalized Friedrichs inequality states that

||v|| K ≤ h Ω π ||∇v|| K , for any v ∈ H 1 (K). (4.28)
Theorem 4.9 (A basic upper bound). Let u be the weak solution of (2.2) and let u ,k h ∈ H 1 0 (Ω) be the solution to (3.5) at the iteration k ≥ 1. Let σ ,k h and Φ ,k h be, respectively, equilibrated flux and saturation reconstructions given in Section 4.2. Then, with κ = (1 -δ ) hΩ π + δ , for δ ∈ {0, 1}, we have

δ ||u -u ,k h || 2 Σ, + (1 -δ ) J u (u ,k h ) + ||u -u ,k h || 2 # ≤ K∈T h h K π ||f -Π h f || K + κ ||Φ ,k h -β(u ,k h )|| K + ||σ ,k h + ξ + ∇ϕ(u ,k h )|| K 2 . ( 4

.29)

Proof. Recalling the flux equilibration identity 2, we get

R u (u ,k h ), v = (f, v) -(ξ, ∇v) -β(u ,k h ), v -∇ϕ(u ,k h ), ∇v , = (f -Π h f, v) + Φ ,k h -β(u ,k h ), v -σ ,k h + ξ + ∇ϕ(u ,k h ), ∇v , (4.30)
where we inserted (σ ,k h , ∇v)+(∇•σ ,k h , v) = 0. The next step is to estimate each term in (4.30) elementwise. For each element K ∈ T h , we use the identity (g -Π h g, v) K = (g -Π h g, v -Π h v) K , for all g ∈ L 2 (Ω), to obtain the following bound using Poincaré-Friedrichs inequality (4.27),

f -Π h f, v K ≤ ||f -Π h f || K h K π ∇v K . (4.31)
The upper bound of the second term will depend on the considered error measure. If δ = 0, we can apply the Cauchy-Schwarz, and the Poincaré-Friedrichs inequality (4.28) to infer

Φ ,k h -β(u ,k h ), v K ≤ ||Φ ,k h -β(u ,k h )|| K h Ω π ||∇v|| K . (4.32)
For δ = 1, we can simply use the Cauchy-Schwarz inequality,

Φ ,k h -β(u ,k h ), v K ≤ ||Φ ,k h -β(u ,k h )|| K ||v|| K . (4.33)
For the last term of (4.30), one can apply the Cauchy-Schwarz inequality to get

σ ,k h + ξ + ∇ϕ(u ,k h ), ∇v K ≤ ||σ ,k h + ξ + ∇ϕ(u ,k h )|| K ||∇v|| K . (4.34)
Therefore, we deduce from (4.30) and the above inequalities that

R u (u ,k h ), v ≤ K∈T h h K π ||f -Π h f || K + κ ||Φ ,k h -β(u ,k h )|| K + ||σ ,k h + ξ + ∇ϕ(u ,k h )|| K ||v|| ,K . (4.35)
By applying the Cauchy-Schwarz inequality, the definitions of the dual norm of the residual (4.3) and (4.6) and (4.8) together with the fact that K∈T ||v|| ,K :

= K∈T ||∇v|| 2 K + δ ||v|| 2 K = v 2
, to get (4.29). We based our developments of Algorithm 3.1 on an energy-type a posteriori error estimate distinguishing the different error components at each nonlinear iteration k ≥ 1. To prove this estimate, let us introduce, for an element K ∈ T h , the local discretization, quadrature, regularization, linearization, and data oscillation estimators as follows:

η ,k disc,K := ||σ ,k h,disc + Π RTN h ξ + ∇ ϕ (u ,k h )|| K + δ ||Φ ,k h,disc -β (u ,k h )|| K (4.36a) η ,k qd,K := ||∇ ϕ(u ,k h ) -∇ϕ(u ,k h )|| K + h K π || β(u ,k h ) -β(u ,k h )|| K , (4.36b)
η ,k reg,K := ||∇ ϕ(u ,k h ) -∇ ϕ (u ,k h )|| K + κ || β(u ,k h ) -β (u ,k h )|| K , (4.36c)
η ,k lin,K := ||σ ,k h,lin || K + κ ||Φ ,k h,lin || K , (4.36d)
η ,k osc,K := h K π ||f -Π h f || K + ||ξ -Π RTN h ξ|| K . (4.36e)
The global versions of the above estimators are given by

η ,k • := K∈T h η ,k •,K 2 1 2
.

(4.37)

Then we have:

δ ||u -u ,k h || Σ, + (1 -δ ) J u (u ,k h ) + ||u -u ,k h || 2 # 1 2 ≤ η ,k disc + η ,k lin + η ,k reg + η ,k osc . (4.38)
Proof. Similar to Theorem 4.9, we use the equilibration property (4.30) to get

R u (u ,k h ), v = (f -Π h f, v) + ξ -Π RTN h ξ, ∇v + β(u ,k h ) -β(u ,k h ), v + Φ ,k h -β(u ,k h ), v + ∇ ϕ(u ,k h ) -∇ϕ(u ,k h ), ∇v -σ ,k h + Π RTN h ξ + ∇ ϕ(u ,k h ), ∇v . (4.39)
where we added and subtracting ϕ(u ,k h ) and β(u ,k h ). To estimate each term, we proceed as in Theorem 4.9 to first get

R u (u ,k h ), v ≤ K∈T h h K π ||f -Π h f || K + ||ξ -Π RTN h ξ|| K + h K π || β(u ,k h ) -β(u ,k h )|| K + ||∇ ϕ(u ,k h ) -∇ϕ(u ,k h )|| K +κ ||Φ ,k h -β(u ,k h )|| K + ||σ ,k h + Π RTN h ξ + ∇ ϕ(u ,k h )|| K ||v|| ,K . (4.40)
By adding and subtracting ϕ (u ,k h ) and β (u ,k h ) to the last two terms respectively, applying the triangular inequality, and using the fact that Φ ,k h = Φ ,k h,disc + Φ ,k h,lin and σ ,k h = σ ,k h,disc + σ ,k h,lin , we obtain

R u (u ,k h ), v ≤ K∈T h η ,k disc,K + η ,k lin,K + η ,k reg,K + η ,k osc,K ||v|| ,K ,
after assembling the terms using (4.36)-(4.37). We conclude the estimate using Cauchy-Schwarz inequality after recalling the identity (4.6) for δ = 0 and (4.8) for δ = 0, while we use the H 1 -energy norm (4.7) so that

K∈T ||∇v|| 2 K + δ ||v|| 2 K = v 2 .

Local efficiency and robustness

This section proves the local and therefrom the global efficiency of our a posteriori error estimate. To accomplish this, we must introduce the Sobolev spaces defined on a patchwise basis We also recall the Poincaré-Friedrichs inequality for functions in H 1 + (ω a ):

H 1 + (ω a ) := {H 1 (ω a ); v = 0 on ∂ω a ∩ ∂Ω} if a ∈ V ext h , {H 1 (ω a ); (v, 1) ωa = 0} if a ∈ V int h , (4.41
||v|| ωa h ωa ||∇v|| ωa , for all a ∈ V. (4.44)
Additionally, we need to introduce the patchwise counterpart of the oscillation estimator (4.36e):

η ,k osc,ωa := K∈T a h h K π ||Π h (ψ a f ) -ψ a f || K + ||ξ -Π h ξ|| K + ||ψ a ξ -Π RTN h (ψ a ξ)|| K , (4.45)
and let η ,k

osc 2 := a∈V h η ,k osc,ωa 2 . (4.46)
Referring to (4.36e), it is evident that for K ∈ T h , we have η ,k osc,K ≤ η ,k osc,ωa due to the partition of unity. Consequently, we have η ,k osc ≤ η ,k osc . We obtain the following stability result by following the procedure in [START_REF] Ern | Discrete p-robust H(div)-liftings and a posteriori estimates for elliptic problems with H -1 source terms[END_REF][START_REF] Mitra | A modified l-scheme to solve nonlinear diffusion problems[END_REF].

Proposition 4.11 (Stability of discretization-flux equilibration)

. Let f ∈ L 2 (Ω) and ξ ∈ L 2 (Ω). Consider u as the weak solution of (2.2). Let u ,k h ∈ L 2 (Ω) be the solution to (3.5) at iteration k ≥ 1 such that ϕ(u ,k h ) ∈ H 1 0 (Ω), and (σ ,k a,disc , Φ ,k a,disc ) be the equilibrated flux and saturation reconstructions given by Reconstruction 2 in accordance with Definition A.2. Let a ∈ V h be fixed, and define τ ,k a,disc

:= Π RTN h (ψ a ξ) + ψ a ∇ ϕ (u ,k h ) ∈ RTN p (T a h ).
Then, for δ ∈ {0, 1}, there holds

||σ ,k a,disc + τ ,k a,disc || 2 ωa + δ r ε,k a,disc 2
ωa 1 2 min v h ∈W a h , ∇•v h =g ,k a -Π h (ψar ,k h ) ||v h + τ ,k a,disc || ωa max |∇v|=1 v∈H 1 + (ωa)\{0} g ,k a -ψ a r ,k h , v ωa -τ ,k a,disc , ∇v ωa , (4.47)
where the constant in depends only on the space dimension d, the shape-regularity constant of T h , and the polynomial degree p. Theorem 4.12 (Local efficiency). Let the assumptions of Proposition 4.11 be satisfied. For a fixed a ∈ V h , the following patchwise efficiency bound holds

||σ ,k a,disc + τ ,k a,disc || 2 ωa + δ r ε,k a,disc 2
ωa 1 2 max v∈H 1 0 (ωa)\{0} R u (u ,k h ), ψ a v ||ψ a v|| + η ,k lin,ωa + η ,k reg,ωa + η ,k qd,ωa + η ,k osc,ωa . (4.48)
Assuming that the adaptive criteria given by (3.8)-(3.10) are satisfied. For all K ∈ T h , let T K = ∪ a∈V K T a h , and ω K be the corresponding open subdomain. Then, we obtain the following local efficiency:

η ,k disc,K + η ,k lin,K + η ,k reg,K max w∈H 1 0 (ω K )\{0} R u (u ,k h ), ψ a w ||ψ a w|| + η ,k qd,ω K + η ,k osc,ω K , (4.49)
where the constant in depends only on the balancing parameters Γ lin and Γ reg , the space dimension d, the shape-regularity constant of T h , and the polynomial degree p.

Lemma 4.13 (Global efficiency).

Let the assumptions of Proposition 4.11 be satisfied. Then, for δ ∈ {0, 1}, the global efficiency can be summarized as

η ,k disc + η ,k lin + η ,k reg δ ||u -u ,k h || Σ, + (1 -δ ) J u (u ,k h ) + ||u -u ,k h || 2 # 1 2 + η ,k qd + η ,k osc . (4.50) Proof. For each K ∈ T h , we recall that η ,k disc,K = ||σ ,k h,disc + Π RTN h ξ + ∇ ϕ (u ,k h )|| K + δ ||r ,k h,disc || K .
Noting that K has (d + 1) vertices collected in the set V K for the simplicial mesh we consider, so that using the partition of unity we have

K∈T h η ,k disc,K 2 K∈T h ||σ ,k h,disc + Π RTN h ξ + ∇ ϕ (u ,k h )|| 2 K + δ ||r ,k h,disc || 2 K , ≤ (d + 1) K∈T h a∈V K ||σ ,k a,disc + τ ,k a,disc || 2 K + δ r ε,k a,disc 2
K , = (d + 1) a∈V h ||σ ,k a,disc + τ ,k a,disc || 2 ωa + δ r ε,k a,disc 2 ωa
.

We now rely on (4.48) to infer

K∈T h η ,k disc,K 2 1 2 a∈V h max v∈H 1 0 (ωa)\{0} R u (u ,k h ), ψ a v ||ψ a v|| + η ,k lin,ωa + η ,k reg,ωa + η ,k qd,ωa + η ,k osc,ωa .
Therefrom, we use the stopping criteria (3.8) and (3.9) with sufficiently small Γ lin and Γ reg to bound the error estimators with the discretization estimator:

η ,k disc + η ,k lin + η ,k reg a∈V h sup v∈H 1 0 (ωa)\{0} R u (u ,k h ), ψ a v ||ψ a v|| + η ,k qd + η ,k osc δ ||u -u ,k h || Σ, + (1 -δ ) J u (u ,k h ) + ||u -u ,k h || 2 # 1 2 + η ,k qd + η ,k osc , (4.51)
Therein, we use the the localization of dual norms with ψ a -Galerkin orthogonality as in [START_REF] Ciarlet | Localization of global norms and robust a posteriori error control for transmission problems with sign-changing coefficients[END_REF], to verify

a∈V h sup v∈H 1 0 (ωa)\{0} R u (u ,k h ), ψ a v ||ψ a v|| R u (u ,k h ) .
The proof is complete by employing (4.6) and (4.8).

Numerical results

In this section, we illustrate our theoretical results through a series of numerical experiments using the finite element method. The discretization is conducted using Lagrange finite elements. Our focus in this paper is on the linearization and regularization processes, with no detailed reporting on quadrature rules, associated errors, or oscillation errors. It's important to note that, in the literature, standard finite elements for (3.2) may require specific quadrature rules or mass-lumped versions of high-order methods to ensure convergence properties.

Given that we exclusively deal with analytical solutions, we compute the exact error components for discretization and linearization as follows:

error k lin = ||u h -u k h || # and error k disc = ||u -u h || # . Specifically, we assume ||ϕ(u) -ϕ(u h )|| H -1 (Ω) = ||ϕ(u) -ϕ(u h )||,
which, in fact, results in a slightly overestimation of the exact error (see Remark 4.3). All presented results are derived using Reconstruction 1. Furthermore, for approximations with orders p ≥ 2, cost-effective error estimations can still be achieved by utilizing a reconstruction of one degree lower.

Test case 1: An illustrative example

Offline optimization of the L-scheme

We set Ω = [0, 1] 2 and let β(u) = u and ϕ(u) = log(u + 1), ξ = 0, and f complying to the solution u(x) = x(1 -x)y(1 -y) on Ω. In particular, due to the Lipschitz continuity of ϕ, there is no need for regularization or variable transformation. The discretization is achieved using first-order Lagrange conforming finite elements, for which we conducted tests with both first-order and one-degree-lower reconstruction (cheaper), both yielding comparable results. In this initial test case, we utilize developed error estimates to compare the L-scheme and Newton's method. Prior to conducting the comparison, we capitalize on the advantage of splitting error sources to optimize the L-scheme parameter.

To achieve this optimization, we implement a "brute force" iterative method that utilizes the linearization estimator to find the optimal parameter L during an offline phase, as detailed in Remark 3.1. This process involves exploring various values of L to determine the one that leads to the minimum iteration number, effectively reducing the overall computational cost. The offline computation of L is carried out on coarse grids, which further contributes to its negligible impact on total computational costs, thanks to the error decomposition technique employed. Overall, this approach allows us to fine-tune the L-scheme, leading to improved efficiency and convergence in the subsequent calculations.

Figure 1 (left) illustrates the L-scheme iteration as a function of the parameter L, comparing both standard and adaptive (with Γ lin ∈ {0.02, 0.2}) stopping criteria. The corresponding error estimators over the L-scheme iterations are depicted in Figure 1 (right) for Γ lin = 0.02 and various values of L (with final iterations utilizing the standard stopping criteria). The optimal value for L is determined to be approximately 0.975, which is then used in Algorithm 3.1 to perform computations on the desired fine grids. Note that this value is close but not equal to the Lipschitz constant, which is 1 for this case. Employing the adaptive stopping criteria leads to significant savings in the number of iterations, and it also helps mitigate the impact of non-optimal values of L. The effectiveness of our balancing strategy is evident from the results shown in Table 1, where the exact and estimated errors exhibit an effectivity index close to 1.

To further compare Newton's method with the L-scheme, we utilize Algorithm 3.1 and plot the various error components for both methods in Figure 2. It is important to note that the discretization estimators are equal in both methods, while the linearization estimator η k lin differs. Therefore, the results presented in Table 1 are equally applicable when using Newton's method. This observation further supports the error splitting techniques employed. In comparing the exact and estimated linearization errors with the two methods, we iterate each of the nonlinear solvers until η k lin ≤ 10 -10 , yielding an almost converged solution. Figure 3 showcases the exact and estimated error components after applying our adaptive stopping criteria, demonstrating a strong agreement between the two. For both methods, the linearization estimator becomes negligible, contributing minimally to the total error.

In conclusion, the above results reaffirm that the optimization of the L-parameter in combination with our adaptive stopping criteria makes the L-scheme an effective alternative to Newton's method.

Mesh adaptivity for an L-shaped domain

In this scenario, we examine an L-shaped domain Ω = (-1, 1) × (-1, 1) \ [0, 1] × [-1, 0] and reverse the choices for β and ϕ from the previous case. Specifically, we set ϕ(u) = u and β(u) = log(u + 1), with ξ = 0.

The exact solution is expressed in polar coordinates:

u(r, θ) = r 2 3 sin 2θ 3 , ∀θ ∈ (0, 3π 2
). We observe that the exact solution features a singularity at the origin, with u ∈ H 5 3 -ε (Ω) for arbitrarily small ε > 0. The corresponding source term is f (u) = log(1 + u). The singularity at the origin penalizes the convergence rate for uniformly refined grids [START_REF] Dolejší | hp$-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF].

For our comparisons, we set the parameter Γ lin = 0.1 and assess the effectiveness of adaptive stopping criteria against standard ones. Subsequently, we examine the performance of adaptive mesh refinement in contrast to standard uniform mesh refinement. Here, regularization is find unnecessary due to the Lipschitz continuity of β and the non-singularity of the derivative. Our adaptive strategy aims to achieve equi-distribution of discretization errors, as outlined in (3.10). In Table 2, we present the total iteration numbers for the L-scheme and Newton methods with and without adaptive stopping criteria. Although this case is not particularly demanding in terms of nonlinear iterations, our stopping criteria lead to a 50% reduction in iterations. Figure 4 illustrates the errors and estimators for uniform and adaptive mesh refinements. A strong agreement is observed between the estimated and exact errors, with adaptive grids showcasing improved performance. With fewer elements, we achieve a superior convergence rate and smaller discretization errors compared to all runs with uniform refinement. This emphasizes the effectiveness of our localized estimators in controlling and reducing singular effects, consequently enhancing convergence rates for such scenarios. In Figure 5, we illustrate the distribution of the exact and estimated errors following the application of our adaptive stopping criteria. An excellent agreement is observed between the estimated and exact error components. In the non-adaptive case, localized and significant discretization errors around the corner are evident. However, these errors are effectively redistributed evenly over the domain by refining and coarsening the grids, in accordance with our discretization criteria.

As a result, the discretization error becomes equidistributed, and the maximum error improves by a factor of 2.5 after grid adaptivity. The linearization errors and estimators for both methods exhibit a sufficiently small contribution from the linearization into the overall error, as required by our adaptive stopping criteria. Clearly, the distribution of the error differs between the two methods, with the Newton method displaying a smoother distribution compared to the L-scheme. The differences in error distribution between the two methods can be attributed to their inherent characteristics. The Newton method, with its second-order convergence and full Jacobian, facilitates a more gradual and smoother adjustment to the correct solution. In contrast, the L-scheme's stabilization amendment may introduce oscillations and localized features, resulting in a less smooth error distribution. adaptively with the solver. We also set β(u) = u, leading to a zero source term in (2.2). In Figure 6 (left), we depict the functions ϕ and its regularized counterpart:

ϕ (u) =        u if u < -, (2 u -1) 1 + 2 if -≤ u ≤ 1 + , u -1 if u ≥ 1 + , (5.4)

Overall algorithm assessment

We used Algorithm 3.1 for p = 1 while setting Γ sw = 0.9, Γ lin = 0.01 and Γ reg = 0.1. The offline computation of the optimal L-scheme parameter resulted in L opt ϕ = 0.6225, as illustrated in Figure 6 (right). We also set a fixed regularization parameter = 0.05. Figure 7 (left) depicts the evolution of the spatial, nonlinear, and regularization error estimators as a function of the number of iterations for a fixed mesh and fixed regularization parameter . As anticipated, the linearization estimator consistently diminishes, whereas the other components remain stagnant, underscoring the distinct error component separation feature. By employing a switching criterion with Γ sw = 0.9, we leverage improved initialization for Newton after completing 9 iterations with the L-scheme. The L-scheme iterations in our adaptive solver are stopped when the linearization estimator falls below the discretization estimator. Although the L-scheme may appear slow initially, it serves as a beneficial initialization for Newton, which exhibits rapid convergence towards the end. Coupled with the adaptive criterion, our solver necessitated a total of 13 iterations with Γ lin = 0.01 and 11 iterations with Γ lin = 0.1. It is crucial to emphasize that a standard Newton solver may experience oscillations or fail to converge when not provided with a sufficiently close starting point, especially for refined grids. Even with a close initialization, our adaptive Newton method significantly reduces the number of unnecessary iterations compared to a conventional criterion based on a fixed threshold of 10 -6 . Given the slower convergence of the L-scheme compared to the faster convergence of Newton, especially towards the end, the selection of Γ sw plays a crucial role. Opting for a larger or smaller Γ sw can introduce numerous additional iterations, including oscillating ones for Newton. In Figure 7 (right), analogous findings are depicted regarding the criterion for the selection of the regularization parameter . As decreases, the regularization error estimator decreases as expected, while the spatial estimator plateaus. The regularization criterion with Γ reg = 0.1 results in stopping the iterations after the fourth step, and with Γ reg = 0.01, it stops after six steps.

Mesh adaptivity

The subsequent set of experiments is designed to evaluate the spatial balancing criterion (mesh adaptivity) by demonstrating its influence on the estimated error and effectivity index in relation to the total number of unknowns, and comparing it to fixed and uniform meshes. We commenced the experiment with a fixed mesh, utilizing consistent parameters for the adaptive, switching, and regularization criteria, as done previously. In Figure 8, we present the distribution of discretization (top left), regularization (top right), and linearization (bottom left) estimators using the adaptive algorithm. The spatial balancing criterion has resulted in a more evenly distributed error, with the maximum error concentrated at the end-points of the interface γ. This distribution is attributed to the still-discontinuous behavior of the boundary conditions initially computed from the exact and non-regularized solution.

Moreover, it is clear that both regularization and linearization errors are sufficiently small to have a negligible influence on the total error. The regularization error predominantly traces the interface γ, where the discontinuity is prominent. To evaluate the overall convergence of the adaptive solver, we present in Figure 8 (bottom right) the error and estimates as a function of the total number of unknowns in both the fully adaptive case and with uniform meshes (with p = 2). Even in the uniform case, Algorithm 3.1 is employed to adaptively control the linearization and regularization processes. Clearly, the adaptive approach outperforms the uniform one, delivering enhanced convergence rates with reduced computational costs. Note that the convergence rates in the uniform case was observed to consistently maintain at or slightly exceed O(h 0.8 ) for p = 1, only slightly declining to around O(h 1.5 ) for p = 2. The convergence rates were clearly enhanced through mesh adaptivity, although the specific rate improvement has not been calculated.

Test case 3: A porous media model

We consider now porous media models by setting β(u) = u and ϕ(u) = |u| m-1 u, where m > 0. This particular selection of β and ϕ leads to a degenerate elliptic equation with intriguing properties depending on the value of m. We set f to comply with the solution u(x) = sin(πx) sin(πy) on Ω = [0, 1] 2 . To investigate the behavior of Algorithm 3.1 in handling the different choices of m, we conduct extensive testing on the model problem. This testing involves numerical simulations with various values of m, spanning both the slow and fast diffusion regimes. By employing a posteriori error estimation, adaptive schemes, and regularization techniques provided by Algorithm 3.1, we aim to obtain accurate and efficient approximations for the diverse cases of interest.

Models for slow regime

When m > 1, the model represents the slow diffusion case, which finds application in gas diffusion within porous media, often referred to as the slow diffusion model. In this scenario, the equation exhibits degeneracy at unknown points where u = 0. This means that the diffusion process becomes slow or even stagnant in regions where the solution approaches zero, leading to spatially localized behavior with potential concentration gradients near these points. Table 3 offers a comprehensive comparison of the L-scheme, Newton, and L-scheme-Newton solvers across various parameter values of m and different switching criteria. The obtained effectivity indices underscore the high quality and robustness of our estimators. Particularly, the estimators demonstrate their efficacy even in cases of reaction dominance (larger m). The consistent trend observed in both our numerically estimated errors and their resulting effectivity indices is an increase as m also increases. The near-singular regions near the boundary do not influence the estimators, highlighting their reliability. We observe that, in almost all cases, the L-scheme-Newton approach outperforms both the individual L-scheme and Newton methods in terms of the total number of iterations. Particularly for fine meshes and m ≥ 3, Newton fails to converge, and the L-scheme alone requires a significant number of iterations. It is noteworthy that the optimal L parameter is found to be close to L opt ϕ = L ϕ /2 = m/2. As an example, we plot in Figure 9 the total iterations number in function of Γ sw for the case of m = 4 and h = 1/16. The results reveal that choosing Γ sw larger or smaller than 0.25 leads to an increase in the number of nonlinear iterations. However, when opting for a small Γ sw , the solver tends to exhibit enhanced stability, effectively mitigating oscillations in the Newton method (left figure). In conclusion, the search for the optimal switching parameter Γ sw alongside the appropriate choice of L significantly enhances the efficiency and stability of the L-scheme-Newton solver.

Models with fast regime

When 0 < m < 1, the model aligns with the fast diffusion case, a prevalent scenario in diverse fields such as plasma physics, the kinetic theory of gases, or fluid transportation in porous media [START_REF] Droniou | The gradient discretization method for slow and fast diffusion porous media equations[END_REF]. In this regime, the equation becomes singular as |u| m-1 becomes unbounded whenever u approaches zero. This behavior indicates that the diffusion process undergoes rapid spreading, causing solutions to diffuse swiftly across the domain and potentially leading to sharp transitions and boundary layer effects. Here, as ϕ is not Lipschitz for arguments approaching 0, we contemplate regularization through ε > 0 and the introduction of the functions ϕ ε : R → R. Specifically, ϕ ε (u) = ε m-1 u if u lies in the interval (0, ε) and ϕ ε (u) = ϕ(u) elsewhere. The Lipschitz constant is L ϕε = ε m-1 , and even ϕ ε is Lipschitz and L b ε = m(m -1)|ε| m-2 . We can employ a transformation by defining w = |u| m-1 u and addressing the problem for w using linear diffusion and a nonlinear reaction term represented by |w| 1/m-1 w. This poses no challenges since 1/m > 1. For this formulation, we apply the L-scheme without regularization. Figure 10 presents a comparative analysis of effectivity indices for both formulations, considering varying values of m and mesh sizes. In the standard formulation (left figure), the estimator tends to overestimate the exact error, especially for small m values, causing the effectivity index to surpass the optimal value of 1. However, as m approaches 1 (linear case), the overestimation diminishes, and the effectivity index progressively converges towards 1. Overall, the effectivity index for the standard formulation remains within a reasonable range close to the optimal value of 1. In contrast, the transformed problem exhibits excellent results (right figure), maintaining an effectivity index consistently close to 1, irrespective of the m value. This highlights the enhanced properties of the transformed problem. Importantly, the transformed formulation demonstrates efficiency across various m values, requiring only a modest number of iterations-typically averaging 2 to 3 iterations. On the other hand, the original formulation demands a higher iteration count, ranging from 5 to 6 iterations when optimizing the stability parameter L β . Without optimizing the Lscheme, the iteration count could significantly increase, as depicted in Figure 11. Figure 12 showcases the spatial distribution of discretization and linearization errors, along with estimators, for the original formulations with m = 0.5 after convergence. The results reveal a remarkable agreement between estimated and exact errors. This robust agreement underscores the effectiveness of employing these estimators to control nonlinearity and certify the accuracy of numerical simulations, especially in scenarios involving degeneracy.

Conclusions

Our work presents an adaptive solver for a class of doubly nonlinear and degenerate elliptic equations, leveraging the strengths of the L-scheme and Newton methods. To this end, we developed a posteriori error estimators providing a guaranteed global upper bound on the error. We also proved the efficiency of the estimates and showed their robustness with respect to mesh size and nonlinearity strengths, and how to evaluate separately the various error components, including regularization and linearization errors.

In our numerical experiments focused on the Stefan problem and scenarios involving fast and slow diffusion in porous media, the adaptive L-scheme-Newton solver consistently outperforms the L-scheme and demonstrates greater stability compared to standard Newton methods, particularly in challenging convergence scenarios. The L-scheme's role as a robust initialization for Newton is crucial in mitigating oscillations and unnecessary iterations. The adaptive solver stands out by providing users with accurate error estimations for different error components at each nonlinear iteration.

Appendix A. Patchwise Euler-Lagrange Equations for Flux Equilibration

The a posteriori error estimates presented in this paper rely on a locally computable equilibrated flux, denoted as σ ,k h , and the post-processed saturation, denoted as Φ ,k h . For any a ∈ V h , the local equilibrated flux σ ,k a and the local post-processed saturation Φ ,k a are defined over vertex patches. These reconstructions correspond to the Euler-Lagrange equations resulting from the constrained minimization within vertex patches.

Definition A.1 (Reconstruction of (σ ε,k a , Φ ε,k a )). For each a ∈ V h , let (σ ε,k a , Φ ε,k a ) ∈ W a h × Q a
h be the minimizers of the following minimization problem:

(σ ε,k a , Φ ε,k a ) := arg min (σ h ,Φ h )∈W a h ×Q a h ∇•σ h +Φ h =b H(σ h ,Φ h ) σ h + λ a 2 ωa + δ Φ h -Π h (ψ a β k-1 (u ,k h )) 2 ωa , (A.1)
where

b := g a + Π h (ψ a β k-1 (u ,k h )) = Π h (ψ a f ) -∇ψ a • (Π h ξ + I ,k h ).
To solve the minimization problem A.1, we employ the method of Lagrange multipliers associated with the Lagrangian function:

L(σ h , Φ h , r) = H(σ h , Φ h ) -(∇ • σ h , r) ωa -(Φ h , r) ωa + (b, r) ωa ∀(σ, Φ, r) ∈ W a h × Q a h × Q a h , (A.2)
where r is the Lagrange multiplier used to impose the needed constraint on (σ h , Φ h ). Denote by (σ ,k a , Φ ,k a , r ,k a ) the saddle point of the Lagrangian function L. Then, for all (v, r, q) ∈ W

a h × Q a h × Q a h , we have (σ ,k a , v) ωa -(∇ • v, r ,k a ) ωa = (-ψ a (ξ + I ,k h ), v) ωa , (Φ ,k a , r) ωa = Π h (ψ a β k-1 (u ,k h )), r ωa + δ (r ,k a , r) ωa , (∇ • σ ,k a , q) ωa + (Φ ,k a , q) ωa = (Π h (ψ a f ), q) ωa -∇ψ a • (Π h ξ + I ,k h ), q ωa .
The Euler-Lagrange equations can then be simplified to two coupled equations as given by (4.22), followed by a post-processing step:

Φ ,k a = Π h (ψ a β k-1 (u ,k h )) + δ r ,k
a , for δ ∈ {0, 1}. The global reconstructions are constructed afterward using the partition of unity property of the hat functions through (4.23).

Definition

A.2 (Reconstruction of (σ ε,k a,disc , Φ ε,k a,disc )). For each a ∈ V h , let (σ ε,k a,disc , Φ ε,k a,disc ) ∈ W a h × Q a
h be the minimizers of the following minimization problem:

(σ ε,k a,disc , Φ ε,k a,disc ) := arg min (σ h ,Φ h )∈W a h ×Q a h ∇•σ h +Φ h =b H disc (σ h ,Φ h ) σ h + λ k, a 2 ωa + δ Φ h -Π h (ψ a β (u ,k h )) 2 ωa , (A.3) where b := g ,k a + Π h (ψ a β (u ,k h )) -ψ a r ,k h = Π h (ψ a f ) -∇ψ a • (Π h ξ + ∇ ϕ (u ,k h )) -Π h (ψ a r ,k h ).
Repeating the same steps as for (A.1), one can obtain the Euler-Lagrange equations (4.24). Therefore, the global constructions are obtained through (4.25).

Case with δ = 0: For this case, we recall that η Case with δ = 1: We leverage the results obtained for δ = 0 along with techniques from [START_REF] Smears | Simple and robust equilibrated flux a posteriori estimates for singularly perturbed reaction-diffusion problems[END_REF] to verify the estimate (4.47) for this case. The idea relies on the fact that (σ ε,k a,disc , Φ ε,k a,disc ) are minimizers of the problem A.3:

,k disc,K = ||σ ,k h,disc + Π RTN h ξ + ∇ ϕ (u ,k h )|| K ,
(σ ε,k a,disc , Φ ε,k a,disc ) = arg min (v h ,Φ h )∈W a h ×Q a h ∇•v h +Φ h =g ,k a +Π h (ψa β (u ,k h ))-Π h (ψar ,k h ) τ ,k a,disc + v h 2 ωa + Φ h -Π h ψ a β (u ,k h ) 2 ωa H disc (v h ,Φ h ) . (B.
3

)
The functional H disc (v h , Φ h ) : W a h × Q a h → R is convex with respect to both variables. Therefore, by fixing Φ * * h := Π h ψ a β (u ,k h )), which is an admissible element of Q a h , we obtain:

min v h ,Φ h H disc (v h , Φ h ) ≤ min v h H disc (v h , Φ * * h ), (B.3) = min v h ∈W a h ∇•v h =g ,k a -Π h (ψar ,k h ) τ ,k a,disc + v h 2 ωa . (B.4)
Here, the constraint is derived from ∇•v h = g ,k a + Π h (ψ a β (u ,k h )) -Π h (ψ a r ,k h ) -Φ * * h = g ,k a -Π h (ψ a r ,k h ) (see Definition A.3). Let us define

v * * h = arg min v h ∈W a h ∇•v h =g ,k a -Π h (ψar ,k h ) ||τ ,k a,disc + v h || ωa , ( = σ ,k a,disc ) (B.5)
Here again, we observe that (τ ,k a,disc , g ,k a -Π h (ψ a r ,k h )) ∈ RTN p (T a h ) × P p (T a h ), and the Neumann compatibility condition (g ,k a -Π h (ψ a r ,k h ), 1) ωa for all interior vertices. Notably, we can verify that the constraint Since (σ ε,k a,disc , Φ ε,k a,disc ) are minimizers of problem A.3, we conclude first that above estimate holds true as well for (σ ε,k a,disc , Φ ε,k a,disc ). The proof is complete.

∇•v * * h + Φ * * h = g ,k a + Π h (ψ a β (u ,k h )) -Π h (ψ a r ,k h ) of problem (B.

Appendix C. Proof of Theorem 4.12

We choose a ∈ V K and consider v ∈ H 1 + (ω a ) then we express the right-hand side of (4.47) as g ,k a -ψ a r ,k h , v with η ,k disc, K for all K ∈ T a h . We conclude the estimate by employing the balancing criteria (3.10) so that η ,k disc, K ≈ η ,k disc,K for all K ∈ T a h . The proof is complete.

|T 2 | + |T 3 | + |T 4 | η ,k osc,ωa    K∈T a h h K π ||Π h (ψ a f ) -ψ a f || K + ||ξ -Π h ξ|| K + ||ψ a ξ -Π h (ψ a ξ)|| K    ||∇v|| ωa (C.

  Ω)-conforming due to the Neumann boundary condition specified in W a h . The properties 2.a.-2.b. use the fact that a∈V K ψ a | K = 1 and a∈V K ∇ψ a | K = 0. To get 2.a. (2.b. uses the same arguments), we have from (4.24b)

Theorem 4 .

 4 10 (A posteriori estimate distinguishing the different error components). Let u be the weak solution of (2.2) and let u ,k h ∈ H 1 0 (Ω) be the solution to (3.5) at the iteration k ≥ 1. Let σ ,k h and Φ ,k h be, respectively, equilibrated flux and saturation reconstructions given in Section 4.2. Then, with κ = (1 -δ ) hΩ π + δ , for δ ∈ {0, 1}, we have

Figure 1 :

 1 Figure 1:Offline computation of optimal L (left) alongside estimated error components (right). The error estimates are derived for h = 1/16 using the L-scheme, where the discretization estimator aligns with either the Newton or L-scheme.

  . Newton lin. est. L-scheme lin. est.

Figure 2 :

 2 Figure 2: Comparing Newton and optimized L-scheme for various parameters Γ lin . Both exact and estimated discretization errors coincides with either the Newton or L-scheme.

Figure 3 :

 3 Figure 3: Comparison of discretization errors (top left) and estimators (top right), along with linearization errors (bottom left) and estimators (bottom right) using the L-scheme. The discretization errors remain identical when employing the Newton method.

Figure 5 :

 5 Figure 5: Spatial distribution of errors (left) and corresponding estimators (right). The top four figures illustrate the discretization error component, showcasing the use of a uniform mesh (top) and an adaptive mesh (bottom).The remaining figures depict the linearization error, illustrating its evaluation through Newton's method (top) and the L-scheme (bottom).

Figure 6 :

 6 Figure6: Illustration of a function ϕ and its regularization ϕ (left), along with the offline computation of the optimal L-scheme parameter. The optimal L is determined for a small regularization parameter = 10 -3 and a uniform coarse mesh with h = 1/8.

Figure 7 :

 7 Figure 7: Evolution of the error estimators as a function of nonlinear iterations with adaptive switching and stopping criteria (left). Evolution of the regularization and discretization estimators with respect to the regularization parameter (right).

Figure 8 :

 8 Figure 8: Spatial distribution of error estimators at the end of simulation: discretization (top left), regularization (top right), and linearization (bottom left). Convergence results comparing uniform and adaptive mesh refinements for p = 1 (bottom right).

Figure 9 :

 9 Figure 9: The ratio η k lin /η k disc as a function of nonlinear iterations (left) for Γsw = 0.1, 0.25, and 0.5. The results highlight the impact of Γsw on cumulative iterations (right).

Figure 10 :Figure 11 : 2 ,

 10112 Figure 10: Effectivity index as a function of model parameter m and mesh size h for the original formulation (left) and the transformed formulation (right). These results were obtained utilizing the L-scheme with an optimized L-parameter computed offline.

Figure 12 :

 12 Figure 12: Comparison of discretization errors (top left) and estimators (top right), along with linearization errors (bottom left) and estimators (bottom right) using the L-scheme for m = 0.5. The discretization errors remain identical when employing the Newton method (results not shown).

3 )Φv∈H 1 +

 31 is met. The arguments used for (B.2) togethet with (B.5) impliesv * * h + τ ,k * * -Π h ψ a β (u ,k h ) (ωa)\{0} g ,k a -ψ a r ,k h , v ωa -τ ,k a,disc , ∇v ωa

5 i=1 T i where T 1 :

 51 = (f, ψ a v) ωa -( β (u ,k h ), ψ a v) ωa -(ξ + ∇ ϕ (u ,k h ), ∇(ψ a v)) ωa , T 2 := (Π h (ψ a f ) -ψ a f, v) ωa , T 3 := (ξ -Π h ξ, v ∇ψ a ) ωa , T 4 := (ψ a ξ -Π RTN h (ψ a ξ), ∇v) ωa , T 5 := -r ,k h , ψ a v ωa .Next, we can easily verify that|T 2 | = K∈T a h (Π h (ψ a f ) -ψ a f, v -Π h v) K using the orthogonality of the L 2 -projection so that |T 2 | ≤ K∈T a h h K π ||Π h (ψ a f ) -ψ a f || K ||∇v|| K byusing the Cauchy-Schwarz and Poincaré-Friedrichs inequalities, as demonstrated in the proof of Theorem 4.10. Similarly, we can prove that |T 3 | ≤ ||ξ -Π h ξ|| ωa ||∇v|| ωa and |T 4 | ≤ ||ψ a ξ -Π h (ψ a ξ)|| ωa ||∇v|| ωa , so that using (4.45)

Table 1 :

 1 Comparison of exact and estimated errors and resulting effectivity index relative to mesh size h. Results are obtained with the L-scheme with Γ lin = 0.02.

	1/h error estimate exact error effectivity ind.
	8	0.03184	0.0296325	1.074
	16	0.01578	0.0149003	1.058
	32	0.00785	0.0074607	1.053
	64	0.00394	0.0037316	1.051

Table 2 :

 2 Comparison of standard and adaptive linearizations

Table 3 :

 3 Comprehensive comparison of nonlinear iterations, error estimates, exact errors, and effectivity indices across different values of m and mesh size h. The results for the L-scheme are obtained with an optimized Lparameter computed offline.

  which corresponds to an equilibration problem based on a pure diffusion problem. We also note that Φ ,k h,disc is undefined by(4.22) and is not needed in the estimates.By orthogonality of the L 2 -projection operatorΠ RTN h onto RTN p (T a h ) together with (4.24a), we can see that (σ ,k a,disc , v) ωa -(∇•v, r ,k a,disc ) ωa = (λ ,k a , v) ωa = -(τ ,k a,disc , v) ωa , for all v ∈ W a h and from (A.3) it's inferred that, for each a ∈ V h , Recall that (τ ,k a,disc , g ,k a -Π h (ψ a r ,k h )) ∈ RTN p (T a h ) × P p (T a h ), and the Neumann compatibility condition (4.26) holds true for all interior vertices. Then, using[START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF] Thm 7] in the case of two space dimensions and[START_REF] Ern | Stable broken h1 and h (div) polynomial extensions for polynomialdegree-robust potential and flux reconstruction in three space dimensions[END_REF] Thm 2.3] in the case of three space dimensions, we obtain, using (B.1) that||σ ,k a,disc + τ ,k a,disc || ωa = min v h ∈W a h , ∇•v h =g ,k a -Π h (ψar ,k h ) ||τ ,k a,disc + v h || ωa

	σ ,k a,disc = arg	min v h ∈W a h	||τ ,k a,disc + v h || ωa .	(B.1)
		∇•v h =g ,k a -Π h (ψar ,k h )
		max |∇v|=1	g ,k a -ψ a r ,k h , v
		v∈H 1 + (ωa)\{0}

ωa -τ ,k a,disc , ∇v ωa . (B.2)

  1) To bound T 5 , we use the Cauchy-Schwarz inequality and Poincaré-Friedrichs inequality to first get |T 5 | h ωa ||r ,k h || ωa ||∇v|| ωa . Now, we use property (2.b.), the triangle inequality, followed by an inverse inequality to get||r ,k h || K = ||∇•σ ,k h,lin + Φ ,k h,lin || K , ≤where the constant depends only on the mesh regularity. Now we are left with the bound of T 1 . As ψ Now, we proceed as in the proof of Theorem 4.10 by (1) using the triangle inequality and (4.1) to incorporate the regularization and quadrature estimators (2) applying the Cauchy-Schwarz and Poincaré-Friedrichs inequalities to bound each term (3) employing the mesh shape-regularity, yielding h ωa ≈ h K and the fact that ||∇(ψ a v)|| ωa ||∇v|| ωa :T 1 R u (u h ), ψ a v + η ,k reg,ωa + η ,k qd,ωa ||∇v|| ωa . u (u h ), ψ a v ||ψ a v|| ,ωa +where we used simply (4.42) for the last inequality for the case of δ = 0 and the fact that the ||ψ a v|| ,ωa (1+h ωa )||∇v|| ωa by combining (4.43)-(4.44) for δ = 1. The last inequality corresponds to (4.48) in which the constant in depends only the space dimension d, the shape-regularity constant of T h , and the polynomial degree p. Now, to prove the local efficiency (4.49), we use the partition of unity to get W replace the above inequality in the patchwise efficiency estimate (4.48), and proceeding by using the stopping criteria (3.8) and (3.9) with small enough user-parameters (Γ lin , Γ reg ) to bound η ,k

	By combining the above inequalities in (B.2), we find out that
					1
	||σ ,k a,disc + τ ,k a,disc || 2 ωa + δ r ε,k a,disc	2 ωa	2	max v∈H 1 0 (ωa)\{0}	R u (u h ), ψ a v |∇v|	+ η ,k lin,ωa + η ,k reg,ωa + η ,k qd,ωa + η ,k osc,ωa ,
						max v∈H 1 + (ωa)\{0}	R u (u h ), ψ a v ||ψ a v|| ,ωa	×	||ψ a v|| ,ωa ||∇v|| ωa	+	•	η ,k •,ωa ,
						max v∈H 1 0 (ωa)\{0}	•	η ,k •,ωa ,	(C.4)
	η ,k disc,K	2		||σ ,k h,disc + Π RTN
			a∈V K
		≤	a∈V K	||σ ,k a,disc + τ ,k a,disc || 2 K + δ r ε,k a,disc	2 K	,
		≤	a∈V K	||σ ,k a,disc + τ ,k a,disc || 2 ωa + δ r ε,k a,disc	2 ωa	.
								lin,	K + η ,k reg,	K
						||∇•σ ,k h,lin || K + ||Φ ,k h,lin || K ,
						K∈T a h
						h -1 K ||σ ,k h,lin || K + ||Φ ,k h,lin || K .	(C.2)
						K∈T a h
	That is we can conclude that
								
		T 5	  K∈T a h	||σ ,k h,lin || K +	h K π	||Φ ,k h,lin || K	  ||∇v|| ωa
								
						||σ ,k h,lin || K + κ ||Φ ,k h,lin || K		||∇v|| ωa ,	(C.3)
				 K∈T a h			
						η ,k lin,ωa

a v ∈ H 1 0 (Ω), we obtain first from (2.2) that

T 1 = β(u) -β (u ,k h ), ψ a v ωa + ∇ ϕ(u) -ϕ (u ,k h ) , ∇(ψ a v) ωa . R h (ψ a ξ) + ψ a ∇ ϕ (u ,k h )|| 2 K + δ ||r ,k a,disc || 2 K ,

Test case 2: A stefan-type model

In this case, the nonlinearity is represented by the Stefan-like function [START_REF] Nochetto | An adaptive finite element method for two-phase stefan problems in two space dimensions. i. stability and error estimates[END_REF], posing a greater challenge for the standard Newton's method:

This increased complexity arises from the potential for the solution u to exhibit discontinuities, as ϕ contains plateaux. This test case, taken from [START_REF] Droniou | Non-conforming Finite Elements on Polytopal Meshes[END_REF], introduces the interface given by s(x, y) = γ where s(x, y) = x+y √ 2 with γ = 1 3 , so that the exact solution

is discontinuous along γ. Note that ϕ(u) is, in fact, continuous since the discontinuity of u coincides with values where ϕ remains constant. Note that Assumptions (A1) to (A4) are satisfied for this case. We use a regularized function ϕ (u) so that its strictly increasing, where the regularization parameter will be set up

Appendix B. Proof of Proposition 4.11

We establish the stability result (4.47) using the approach presented in [START_REF] Smears | Simple and robust equilibrated flux a posteriori estimates for singularly perturbed reaction-diffusion problems[END_REF] and, more specifically, in [START_REF] Ern | Discrete p-robust H(div)-liftings and a posteriori estimates for elliptic problems with H -1 source terms[END_REF]. These references utilize stability bounds and quasi-minimizers for the patchwise discretization-flux equilibration problem (4.24).