Jean-Baptiste Clément
email: jean-baptiste.clement@univ-tln.fr

Frederic Golay
email: frederic.golay@univ-tln.fr

Mehmet Ersoy
email: mehmet.ersoy@univ-tln.fr

Damien Sous
email: damien.sous@mio.osupytheas.fr

Frédéric Golay

An adaptive strategy for discontinuous Galerkin simulations of Richards' equation: application to multi-materials dam wetting

Keywords: Unsaturated porous media, Weighted Discontinuous Galerkin method, Adaptive Mesh Refinement, A posteriori error estimation, Adaptive time stepping, BDF

Introduction

Predicting variably-saturated flows in porous media is a major issue for many fields in science and engineering. For example, such flows arise in soil physics, hydrogeology, environment, agriculture or oil industry for problems like subsurface contaminant transport, petroleum reservoir, water resources, surface ponding, etc.

In the present work, flows in variably-saturated porous media are described by Richards' equation. What makes Richards' equation very attractive is that it models the porous medium as a whole part including both saturated and unsaturated zones. Richards' equation is a nonlinear parabolic equation which can degenerate into an elliptic equation under complete saturation condition. Richards' equation has been extensively used for numerical simulations by the hydrogeology community, see e.g. [START_REF] Thoms | User's guide to the Variably Saturated Flow (VSF) process to MODFLOW[END_REF][START_REF] Szymkiewicz | Modelling Water Flow in Unsaturated Porous Media[END_REF]. Despite this research effort, many numerical challenges remain for solving Richards' equation with some simulations being still unreliable and/or expensive [START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF][START_REF] Zha | Review of numerical solution of Richardson-Richards equation for variably saturated flow in soils[END_REF]. This numerical complexity prevented a more general use of Richards-based model for a number of applications. Indeed, the solution of Richards' equation involves sharp wetting fronts which evolve both in space and time and are difficult to resolve. Besides, the simulation must treat simultaneously unsaturated/saturated regions of parabolic/elliptic natures, several porous media of different hydraulic properties and possibly fast-changing boundary conditions. These aspects can also lead to steep gradients acting like discontinuities. Considering the set of nonlinear functions for hydraulic properties, it is seldom possible to foresee the behaviour of Richards' equation and getting numerical solutions can be hard to achieve. Especially, numerical schemes often fail to converge or need fine discretization, making computation costly. Many numerical methods are dealing with Richards' equation to treat wetting fronts and heterogeneous soils in various contexts. Each of them shows advantages and drawbacks and that is the reason why research is still undergoing to push back the limits for the use of Richards' equation. For instance, to handle layered soils, some recent publications include a Kirchhoff integral transform-based method [START_REF] Suk | Numerical solution of the Kirchhoff-transformed Richards equation for simulating variably saturated flow in heterogeneous layered porous media[END_REF], a transversal method of lines [START_REF] Berardi | A mixed MoL-TMoL for the numerical solution of the 2D Richards' equation in layered soils[END_REF], a domain decomposition technique with L-scheme linearization [START_REF] Seus | A linear domain decomposition method for partially saturated flow in porous media[END_REF] and even upscaled models of Richards' equation for fractured porous media [START_REF] Kumar | Formal upscaling and numerical validation of unsaturated flow models in fractured porous media[END_REF].

In this paper, a discontinuous Galerkin (DG) method is chosen to solve Richards' equation. DG methods are based on a variational formulation in an element-wise fashion, sharing advantages both with finite elements and finite volumes methods. In particular, they are locally conservative which is crucial in fluid dynamics [START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF]. Moreover, the nature of DG formulation enables to work on non-conforming mesh and to change locally the degree of polynomial approximation. This is an important benefit since adaptive mesh refinement (AMR), the socalled h-adaptation, and high-order accuracy, the so-called p-adaptation, become possible [START_REF] Dolejší | Discontinuous Galerkin Method[END_REF], promoting a growing use of DG methods for transport phenomena in porous media such as two-phase flow problems [START_REF] Klieber | Adaptive simulations of two-phase flow by discontinuous Galerkin methods[END_REF][START_REF] Bastian | A fully-coupled discontinuous Galerkin method for two-phase flow in porous media with discontinuous capillary pressure[END_REF]. However, DG methods remain sparsely applied to Richards' equation. In 2007, Li et al. solved Richards' equation in 1D with a local discontinuous Galerkin method [START_REF] Li | Adaptive local discontinuous Galerkin approximation to Richards' equation[END_REF][START_REF] Li | Local discontinuous Galerkin approximations to Richards' equation[END_REF]. Sochala discretized Richards' equation through a SIPG mixed DG formulation in 2008 [START_REF] Sochala | Méthodes numériques pour les écoulements souterrains et couplage avec le ruissellement[END_REF][START_REF] Sochala | Mass conservative BDF-discontinuous Galerkin/explicit finite volume schemes for coupling subsurface and overland flows[END_REF]. Recently, in 2019, Dolejší et al. [START_REF] Dolejší | Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF] proposed a space-time DG method for solving Richards' equation.

To reach robustness and accuracy, an adaptive strategy is developed in this paper. Richards' equation is known to be a stiff differential equation with difficult convergence [START_REF] Paniconi | A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems[END_REF][START_REF] Lehmann | Comparison of Iterative Methods for Improved Solutions of the Fluid Flow Equation in Partially Saturated Porous Media[END_REF][START_REF] List | A study on iterative methods for solving Richards' equation[END_REF] so the time discretization and the nonlinear solver have to be addressed carefully. Adaptive mesh refinement will be used to capture moving wetting fronts thanks to a posteriori estimation. h-Adaptation is employed for two-phase flow in porous media [START_REF] Klieber | Adaptive simulations of two-phase flow by discontinuous Galerkin methods[END_REF] but also for Richards' equation [START_REF] Miller | A spatially and temporally adaptive solution of Richards' equation[END_REF][START_REF] Li | Local discontinuous Galerkin approximations to Richards' equation[END_REF]. Making the most of DG methods flexibility, adaptive mesh refinement is combined with a weighted discontinuous Galerkin (WDG) framework which allows discontinuity in the solution according to the nonlinear diffusivity. Following earlier works, this approach was formulated by Ern, Di Pietro and other collaborators [START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF][START_REF] Di Pietro | Discontinuous Galerkin Methods for Anisotropic Semidefinite Diffusion with Advection[END_REF] and by Proft and Rivière under the name improved and adapted discontinuous Galerkin methods [START_REF] Proft | Analytical and numerical study of diffusive fluxes for transport equations with near-degenerate coefficients[END_REF][START_REF] Proft | Discontinuous Galerkin methods for convection-diffusion equations with varying and vanishing diffusivity[END_REF]. Application of such strategy for Richards' equation, combining adaptive mesh refinement, a posteriori estimation and WDG framework, is the main novelty of the present study.

The Richards' equation framework is first recalled in Section 2. In Section 3, Richards' equation discretization is presented through a DG framework and then solved by a nonlinear iterative process aiming to robustness. Section 4 is dedicated to the adaptive strategy proposed to improve solving of Richards' equation. Several numerical experiments are discussed in Section 5. In particular, the modeling strategy is tested on a challenging benchmark case of full-scale wetting of a multi-materials dam. Finally, Section 6 is devoted to some concluding remarks.

Model problem 2.1 Richards' equation

Richards' equation is a classical nonlinear parabolic equation used to describe flow in unsaturated/saturated zones of an aquifer. Derivation of Richards' equation is described for instance in [START_REF] Szymkiewicz | Modelling Water Flow in Unsaturated Porous Media[END_REF]. Variables in play are θ the water content [-], ψ the pressure head [L] and K the hydraulic conductivity tensor [L • T -1 ]. The mixed formulation is selected in the present study owing to its versatility [START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF]:

∂ t θ(ψ) -∇ • (K(ψ)∇(ψ + z)) = 0. ( 1
)
Alternative pressure-based or saturation-based formulations have been discarded because they are either non-conservative or undefined for complete saturation and heterogeneous soils [START_REF] Celia | A general massconservative numerical solution for the unsaturated flow equation[END_REF].

The mixed formulation can be rewritten in the hydraulic head form which is more common in hydrology:

∂ t θ(h -z) -∇ • (K(h -z)∇h) = 0, (2)
where h = ψ + z is the hydraulic head [L]. One can also add a source/sink term Q [T -1 ] to Richards' equation to model various processes, like bacteria colony or plant roots uptake, or to couple Richards' equation with free surface flow.

Solving Eq. ( 1) requires two constitutive laws: one for hydraulic conductivity and one for water content. Several models have been proposed, depending on the hydraulic properties of the porous medium.

The hydraulic conductivity K is generally supposed to react to saturation identically for each space direction. This leads to write:

K(ψ) = K s K r (ψ), (3)
where K s the intrinsic or saturated hydraulic conductivity tensor [L • T -1 ] and K r the relative hydraulic conductivity [-].

For practical purposes, the water content is often described in terms of effective saturation S e [-]:

S e (ψ) = θ(ψ) -θ r θ s -θ r , (4)
where θ s denotes the saturated water content [-] and θ r the residual water content [-], corresponding to the maximal and minimal saturations, respectively. The hydraulic properties present two different behaviours depending on whether the porous media is saturated (ψ ≥ 0) or not (ψ < 0):

S e (ψ) = 1 if ψ ≥ 0, S e otherwise,
and

K r (ψ) = 1 if ψ ≥ 0, K r otherwise. (5)
S e and K r are monotonic increasing functions of pressure head ψ in the unsaturated zone. The water table corresponds to ψ = 0 by definition and is considered belonging to the saturated zone. The capillary fringe is the layer above water table where water is raised due to capillary actions. It belongs to the unsaturated zone and there is no standard definition for its upper limit. Throughout this paper, several constitutive laws will be used to model hydraulic properties in the unsaturated zone. They are compiled in Tab. 1.

Name Expression Parameters

Gardner-Irmay relations (1958) [START_REF] Irmay | On the hydraulic conductivity of unsaturated soils[END_REF][START_REF] Gardner | Some Steady-State Solutions of the Unsaturated Moisture Flow Equation with Application to Evaporation from a Water Table[END_REF] S e = e Van Genuchten-Mualem relations (1980) [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF][START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils[END_REF]]

S e = (1 + (α|ψ|) n ) -m K r = S l e 1 -1 -S 1 m
e m 2 l = 0.5 for Mualem [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF] 1 for Burdine [START_REF] Burdine | Relative Permeability Calculations From Pore Size Distribution Data[END_REF] : pore connectivity [-] α: parameter linked to air entry pressure inverse for particular set of parameters, when ψ → 0 -, constitutive laws may exhibit very steep gradients.

[L -1 ] n > 1: pore-size distribution [-] m = 1 - 1 n : pore-size distribution [-]
These constitutive laws are mainly responsible for the numerical challenges of Richards' equation, such as nonlinearities, degenaracies and instabilities, which are often observed in the presence of material heterogeneities or dynamic boundary conditions like seepage boundary condition. Special treatment may be done to regularize the constitutive relations [START_REF] Jäger | Solution of doubly nonlinear and degenerate parabolic problems by relaxation schemes[END_REF][START_REF] Sorin | Error Estimates for a Time Discretization Method for the Richards' Equation[END_REF][START_REF] Florin | Newton-Type Methods for the Mixed Finite Element Discretization of Some Degenerate Parabolic Equations[END_REF]. For instance, Dolejší et al. [START_REF] Dolejší | Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF] modify slightly the functions to avoid some types of degeneracy mentioned above and improve the convergence properties. It can cause difficulties for numerical analysis but, in practice, we find that the convergence of our solving algorithm described in Section 3 can be ensured choosing carefully the parameters of the computations. As a result, no modification on the constitutive relations is done in this study.

Seepage boundary condition

The seepage boundary condition is specific to subsurface model [START_REF] Scudeler | Examination of the seepage face boundary condition in subsurface and coupled surface/subsurface hydrological models[END_REF]. This condition is used to model the interface between a porous medium and the atmosphere. If the porous medium is saturated and an outflow occurs, then water pours out at atmospheric pressure: ψ = 0 ⇐⇒ h = z. Otherwise, the interface acts as an impervious boundary and there is no flux. This condition mimics an outflow condition. Figure 1 depicts the situation. The treatment of the seepage boundary condition is difficult because the length of seepage face is unknown a priori : it depends on the sought solution. There are many ways to express and realize the seepage boundary condition. A simple and common approach is to refer to its base definition by considering a switch between a Dirichlet ψ D and a Neumann q N boundary condition [START_REF] Richard | Some new procedures for numerical solution of variably saturated flow problems[END_REF][START_REF] Beaugendre | A seepage face model for the interaction of shallow water tables with the ground surface: Application of the obstacle-type method[END_REF][START_REF] Šimůnek | The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and MultipleSolutes in Variably Saturated Media[END_REF][START_REF] Scudeler | Examination of the seepage face boundary condition in subsurface and coupled surface/subsurface hydrological models[END_REF]:

h = z if h ≥ z and -K(h -z)∇h • n > 0, -K(h -z)∇h • n = 0 otherwise. ( 6
)
The realization of such condition should be done inside the nonlinear iterative process according to the previous solution guess at a local level. However, it is not always the case because of the involved numerical methods and discretization schemes. Then, the exit point is searched iteratively until a convergence criterion is reached according to different techniques. Alternatively, the seepage boundary condition may be interpreted as a nonlinear Robin boundary condition:

1 S (h)(h -z) -(1 -1 S (h))K(h -z)∇h • n = 0, (7)
where the seepage indicator function is:

1 S : Γ S → {0, 1} (8)
h → 1 if h ≥ z and -K(h -z)∇h • n > 0, 0 otherwise. ( 9
)
This compact formulation will be used here for the discretization to incorporate directly the seepage boundary condition into the weak formulation, keeping in mind that it encompasses a mix of a Dirichlet and a Neumann boundary conditions whose actual realization is done such as in Eq. ( 6).

3 Numerical methods

Discontinuous Galerkin discretization

More careful developments about discontinuous Galerkin (DG) methods can be found in Rivière [START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF] or Dolejší and Feistauer [START_REF] Dolejší | Discontinuous Galerkin Method[END_REF]. Let d ∈ {1, 2, 3} be the dimension. The porous medium is represented by the computational domain Ω ⊂ R d of boundary ∂Ω and the final time is

T ∈ R * + .
The boundary ∂Ω is subdivided into three mutually disjoint boundaries, ∂Ω = Γ D ∪ Γ N ∪ Γ S , corresponding to the Dirichlet, Neumann and seepage boundary conditions respectively. The following problem is considered:

find h(x, t) : Ω × (0, T ) -→ R such that                ∂ t θ(h -z) -∇ • (K(h -z)∇h) = 0, in Ω × (0, T ), h = h 0 , in Ω × {0}, h = h D , on Γ D × (0, T ), -K(h -z)∇h • n = q N , on Γ N × (0, T ), 1 S (h)(h -z) -(1 -1 S (h))K(h -z)∇h • n = 0, on Γ S × (0, T ). ( 10
)
The time duration (0, T ) is subdivided into N time intervals such that 0 = t

0 < t 1 < • • • < t N = T . Let n ∈ N, 0 < n ≤ N : if the time interval T n = [t n-1 , t n ] is considered, the 120 corresponding time step is τ n = t n -t n-1 .
Ω is subdivided into N E mutually disjoint polygonal elements E forming a mesh. The mesh at the time partition T n is denoted by

E n h = {E n } E⊂Ω . So, Ω = E∈E n h E.
Moreover, the boundary of one element E ∈ E n h is ∂E while its diameter h E is defined by:

h E := sup x,y∈E |x -y|. (11)
For every mesh E n h , the space step is set as:

h n := max E∈E n h (h E ). (12)
The set of all open faces of all elements E ∈ E n h is denoted by F n h . Furthermore, these definitions stand:

F n h,D := F ∈Γ D F, F n h,N := F ∈Γ N F, F n h,S := F ∈Γ S F (13)
F n h,B := F n h,D ∪ F n h,N ∪ F n h,S , F n h,I := F n h \ F n h,B . (14)
For each face F ∈ F n h , there is a unit normal vector n F . For F ∈ F n h,I , its orientation is arbitrary but kept fixed and, for F ∈ F n h,B , it is oriented outward. Let two neighbouring elements E l and E r sharing one face F ∈ F n h,I . There are two traces of a function u along F denoted by u l and u r :

x ∈ F, u l (x) := lim ε→0 + u(x + εn F ), x ∈ F, u r (x) := lim ε→0 - u(x + εn F ). (15)
Formally, the jump and the mean of the function u across a face F ∈ F n h,I are defined respectively by:

u := u l -u r , ⦃u⦄ := 1 2 (u l + u r ). ( 16
)
In the case that F ∈ F n h,B , u = ⦃u⦄ = u l where u l is the trace of the function u from the element which holds F . Throughout this paper, the quantity h F denotes the length of F ∈ F n h and p F the mean of polynomial degree of neighbouring elements:

h F := |F |, p F := ⦃p 2 ⦄. (17)
In the following derivation of the DG variational formulation, several interior penalty methods are considered. Their numerical properties are slightly different. Further information are available in [START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF][START_REF] Dolejší | Discontinuous Galerkin Method[END_REF] for example. They use discontinuous approximations over the mesh E n h defining a so-called broken Sobolev space H s E n h with s > 3/2. The solution is sought in a subspace of this broken Sobolev space, taken to be:

S p E n h := v ∈ L 2 (Ω) : v E ∈ P p (E), ∀E ∈ E n h , (18)
where P p (E) denotes the space of polynomial functions on E of degree less than p ∈ N. Equation [START_REF] Casulli | A Nested Newton-Type Algorithm for Finite Volume Methods Solving Richards' Equation in Mixed Form[END_REF] On the other hand, two penalty terms are considered. One mimics the solution continuity by constraining the interior solution jump to be zero, and one enforces the Dirichlet boundary condition, including the remaining part of the seepage boundary condition. If the solution is a continuous function satisfying the Dirichlet boundary conditions, the two penalty terms are vanishing so they can be added to the formulation:

J I (h, v) = F ∈F I h F I F h v dF, J D (h, v) = F ∈F D h F D F (h -h D )v dF. (19)
The interior and Dirichlet penalty weights are set as:

I F = σ I F γ F µ F and D F = σ D F γ F µ F . ( 20
)
where the face measure is defined as µ F = h F p F 2 and where σ I F , σ D F and γ F are suitably chosen positive constants. σ I F and σ D F are user-defined parameters which must be above a threshold value to provide coercivity and guarantee stability for some of the DG methods [START_REF] Epshteyn | Estimation of penalty parameters for symmetric interior penalty Galerkin methods[END_REF]. γ is a diffusion penalty coefficient set to one generally whereas some studies consider other values, see e.g. [START_REF] Schötzau | A robust a-posteriori error estimator for discontinuous Galerkin methods for convection-diffusion equations[END_REF].

A first weak formulation rises but it is not is not symmetric between the unknown h and the test function v. Making the formulation symmetrical can be useful to preserve the natural symmetry in the discrete diffusion operator, to use appropriate solvers or to enhance numerical properties of the formulation. This can be done by adding the corresponding symmetric term which vanishes because h = 0 if h is a sufficiently smooth function.

Finally, the space semidiscretization reads:

find h ∈ S p E n h such that ∀v ∈ S p E n h , m h,n (∂ t θ(h -z), v) + a h,n (h, v) = l h,n (v), (21)
where the DG bilinear forms, m and a, and linear form l are:

m h,n (∂ t θ(h -z), v) = E∈E h E ∂ t θ(h -z)v dE, (22)
a h,n (h, v) = E∈E h E K(h -z)∇h • ∇v dE - F ∈F I h F ⦃K(h -z)∇h • n F ⦄ v dF - F ∈F D h F K(h -z)∇h • n F v dF + F ∈F D h F D F hv dF - F ∈F S h F 1 S (h)K(h -z)∇h • n F v dF + F ∈F I h F I F h v dF + F ∈F S h F D F 1 S (h)hv dF -Θ F ∈F I h F ⦃K(h -z)∇v • n F ⦄ h dF -Θ F ∈F D h F K(h -z)∇v • n F h dF -Θ F ∈F S h F 1 S (h)K(h -z)∇v • n F h dF, (23)
l h,n (v) = F ∈F D h F D F h D v dF + F ∈F S h F D F 1 S (h)zv dF - F ∈F N h F q N v dF -Θ F ∈F D h F K(h -z)∇v • n F h D dF -Θ F ∈F S h F 1 S (h)K(h -z)∇v • n F z dF. (24)
Θ = {-1, 0, 1} is a constant leading to different symmetric versions of the DG formulation listed in Tab. 2.

Symmetrization Θ = -1 Θ = 0 Θ = 1 Penalization ∀F ∈ F h , σ I F = σ D F = 0 OBB method - global element method ∀F ∈ F h , σ I F = 0, σ D F = 0 NIPG IIPG SIPG
NIPG: non-symmetric interior penalty Galerkin SIPG: symmetric interior penalty Galerkin IIPG: incomplete interior penalty Galerkin OBB method: Oden-Baumann-Babuška method Table 2: Different types of DG methods.

The seepage boundary condition is treated numerically thanks to the indicator function 1 S which alternates dynamically between the Dirichlet part when 1 S = 1 and the Neumann part for 1 S = 0, directly in the usual weak framework of the DG methods. This treatment can model many physical situations, like multiple seepage faces simultaneously, because it does not require any assumption about the mesh or seepage face [START_REF] Scudeler | Examination of the seepage face boundary condition in subsurface and coupled surface/subsurface hydrological models[END_REF]. This treatment was also done by Dolejší et al. in [START_REF] Dolejší | Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF] where they pointed out this is a great benefit compared to conforming FE methods which need to modify the FE space or use numerical tricks, resulting in difficult implementation. However, the indicator function 1 S is not differentiable which causes difficulties for the numerical analysis. A regularization of this function can be done [START_REF] Schweizer | Regularization of outflow problems in unsaturated porous media with dry regions[END_REF][START_REF] Sorin | Regularization Schemes For Degenerate Richards Equations And Outflow Conditions[END_REF]. For example, Dolejší et al. [START_REF] Dolejší | Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF] smoothed the indicator function 1 F to avoid a brutal switch between the values 0 and 1. In this work, no regularization is performed because numerical experiments do not show major convergence solver problems. Besides, DG methods handle already discontinuities and weak penalty naturally. This makes the numerical scheme even more simple.

Time discretization

To get the fully discrete DG formulation, the time derivative has now to be discretized. Such a procedure which takes care of space discretization before time discretization is called the method of lines. It is used commonly by others studies for Richards' equation, see e.g. [START_REF] Miller | A spatially and temporally adaptive solution of Richards' equation[END_REF][START_REF] Li | Adaptive local discontinuous Galerkin approximation to Richards' equation[END_REF][START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF]. Backward differentiation formula (BDF) are implicit numerical integration methods especially used for stiff differential equations thanks to their wide region of stability. These linear multistep methods use q already computed solutions to produce a method of order q ≤ 6. Analysis can be found in [START_REF] Süli | An Introduction to Numerical Analysis[END_REF][START_REF] Dolejší | Analysis of a BDF-DGFE scheme for nonlinear convection-diffusion problems[END_REF]. The 1-step BDF is the backward Euler scheme. For this study, BDF are interesting because they provide high-order accuracy in time which is suitable if one wants to reach high-order in space without loss of gain. The following notation will be used for any function u ∈ L 2 (0, T ; Ω):

∀n ∈ N + , u n := u(x, t n ). ( 25
)
The DG formulation in Eq. ( 21) is integrated over T n and the time integral is approximated by the BDF methods of 0 < q < 7 steps:

find a sequence of (h n ) n∈N + ∈ S p E n h such that      h 0 = h 0 , ∀v ∈ S p E n h , m h,n q k=0 α q,k τ n θ h n+1-k -z , v + a h,n h n+1 , v = l h,n v, t n+1 , (26)
where BDF coefficients are given in Tab. 3.

q 1 1 -1 2 3 /2 -2 1 /2 3 11 /6 -3 3 /2 -1 /3 4 25 /12 -4 3 -4 /3 1 /4 5 137 /60 -5 5 -10 /3 5 /4 -1 /5 6 49 /20 -6 15 /2 -20 /3 15 /4 -6 /5 1 /6 α s,k 0 1 2 3 4 5 6 k
Table 3: BDF coefficients α q,k for all methods.

However, stability of BDF methods has some restrictions. Only BDF methods of order 1 and 2 are A-stable according to the second Dahlquist barrier for implicit linear multistep methods as defined and proved in [START_REF] Germund | A special stability problem for linear multistep methods[END_REF]. BDF methods of order 3 to 6 are conditionally stable or A(α)-stable where α decreases with the order, see e.g. [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems[END_REF][START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF]. Then, time stepsizes have to be small enough to make the method stable. BDF methods of order q > 6 are unconditionally unstable so they cannot be used.

BDF methods presented here are based on fixed time steps so initialization and a change of time stepsize has to be treated carefully. For a scheme with q-step BDF, the first time step is subdivided into as many times as necessary to compute the sub-time steps with an increasing sequence of lower-order BDF methods. The remaining next macro-time steps are still subdivided while the number of computed solution, stored gradually, is not sufficient to go on with the qstep BDF method. Working with BDF methods of variable time stepsize is possible using divided differences-based formula but with some restrictions on time step variations to keep stability. This is done for Richards' equation in [START_REF] Clément | Numerical simulation of flows in unsaturated porous media by an adaptive discontinuous Galerkin method: application to sandy beaches[END_REF] following the method of Hay et al. [START_REF] Hay | hp-Adaptive time integration based on the BDF for viscous flows[END_REF]. For this study, the backward Euler scheme will be mainly used because it provides good stability properties and remains simple to work with adaptive time stepping. For Tracy's benchmark, the time step is kept unchanged to deal with high-order BDF methods which are necessary to decrease time error regarding space error. This study uses the mixed form of Richards' equation (physically mass conservative) and DG methods (numerically local mass conservative). Then, mass conservation should be very good which is checked, globally at machine precision and locally with small discrepancies, in [START_REF] Clément | Numerical simulation of flows in unsaturated porous media by an adaptive discontinuous Galerkin method: application to sandy beaches[END_REF].

Linearization

Richards' equation is a nonlinear equation usually solved by an iterative procedure such as fixed-point iteration (also called Picard's) or Newton-Raphson method whose choice is determining for computation time performances and convergence. Lots of derived methods exist: modified Newton-Raphson, mixed Picard-Newton, quasi-Newton, L-scheme, etc. Studies have been carried out to compare these methods, and others, for solving Richards' equation for various problems [START_REF] Paniconi | A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems[END_REF][START_REF] Lehmann | Comparison of Iterative Methods for Improved Solutions of the Fluid Flow Equation in Partially Saturated Porous Media[END_REF][START_REF] List | A study on iterative methods for solving Richards' equation[END_REF]. They emphasize that fixed-point iteration and even Newton-Raphson scheme are very sensitive and do not converge systematically according to Richards' equation's formulations, initial and boundary conditions and because of the nonlinearities of constitutive laws. In particular, the nonlinear iterative solver can oscillate between two solutions. There are still recent and significant works on this topic like Anderson acceleration [START_REF] Lott | An accelerated Picard method for nonlinear systems related to variably saturated flow[END_REF], nested Newton's type algorithm [START_REF] Casulli | A Nested Newton-Type Algorithm for Finite Volume Methods Solving Richards' Equation in Mixed Form[END_REF] or modified L-schemes [START_REF] Mitra | A modified L-scheme to solve nonlinear diffusion problems[END_REF]. For this study, a damped Newton-Raphson method and a fixed-point method were implemented. Let m be the number of the nonlinear iteration. The residual of Eq. ( 26) is expressed as:

r h,n h n+1 , v := m h,n q k=0 α s,k τ n θ h n+1-k -z , v + a h,n h n+1 , v -l h,n v; t n+1 . (27)
The difference between two successive iterations is written δ n+1,m h := h n+1,m+1 -h n+1,m so the iterative procedure reads:

   dr h,n h n+1,m , v dh n+1,m δ n+1,m h = -r h,n h n+1,m , v , h n+1,m+1 = h n+1,m + δ n+1,m h . ( 28
)
Underlying ideas are already described in [START_REF] Lehmann | Comparison of Iterative Methods for Improved Solutions of the Fluid Flow Equation in Partially Saturated Porous Media[END_REF] and have just been adapted here to the DG system. Roughly speaking, a fixed-point method is obtained from the Newton-Raphson method by avoiding the first derivative terms inside dr h,n h n+1,m , v dh n+1,m . The damped Newton-Raphson method relaxes the increment δ n+1,m h while the new residual is greater than the previous one. Such a procedure is described in [START_REF] Dolejší | Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF]. One important choice for nonlinear iterative process is the stopping criterion which is set for this study as:

r h,n h n+1,m , v L 2 (Ω) a h,n (h n+1 , v) L 2 (Ω)
< ε 1 and

δ n+1,m h L 2 (Ω) h n+1,m L 2 (Ω) < ε 2 , (29)
where ε 1 , ε 2 are a user-defined tolerances. These two criteria are relative in order to be independent from the characteristic quantities of the problem. The first criterion is residual and the second one is incremental. The numerical tests of Section 5 were carried out with the fixed-point method because it was more robust concerning initialization.

For each nonlinear iteration, the linear system has to be constructed and solved. We use an in-house direct solver based on LU decomposition and working with skyline matrix storage.

Adaptive time stepping

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

For this study, the time step is adjusted heuristically according to the previous number of iterations N it from the nonlinear solver such as [START_REF] Bergamaschi | Mixed finite elements and Newton-type linearizations for the solution of Richards' equation[END_REF][START_REF] Thoms | User's guide to the Variably Saturated Flow (VSF) process to MODFLOW[END_REF]. Indeed, the time stepsize has a direct effect on the convergence of the solver. The simulations begin with a time step τ 0 . Then, the future time steps are calculated accordingly to the following rule: each time step is kept if convergence at the previous time iteration is achieved between m it and W it nonlinear iterations, it is increased by an amplification factor λ amp if the convergence requires fewer than m it nonlinear iterations and it is decreased by a reduction factor λ red if the convergence requires more than M it nonlinear iterations. In the case where the convergence is not fulfilled because of a solver failure (poor initial guess, bad condition number) or because the nonlinear iterations go over a prescribed maximum bound W it , the time step is started again (back-stepping) using a stepsize reduced by λ red . The time stepping scheme to choose the next time step τ n+1 is calculated from the previous one τ n can be summarised with:

           τ n+1 =      λ amp τ n if N it ≤ m it , τ n if m it < N it ≤ M it , λ red τ n if M it < N it ≤ W it , τ n = λ red τ n if W it < N it or if the solver has failed (time step is started again). (30)
The factors 0 < λ red < 1 and 1 < λ amp as well as the threshold values 0 < m it ≤ M it < W it are prescribed by the user and found empirically. Default values are λ amp = 2, λ red = 0.5, m it = 3, M it = 7 and W it = 10. A minimum time step can be defined to avoid excessive small time steps. With this approach, the nonlinear solver is more robust because the time step is adjusted until success of convergence independently on τ 0 . Thanks to the amplification/reduction coefficients, the time step is adjusted smoothly. Nevertheless, the method depends on fixed empirical parameters and does not provide an optimal time step. The resulting loss/gain in computational time is difficult to assess in regards to balance between nonlinear iterations number and time steps length.

4 Adaptive strategy

Adaptive Mesh Refinement

The Adaptive Mesh Refinement techniques (AMR) are now widely used and have since proven their efficiency, whether on 2D or 3D mesh, structured or unstructured mesh, conforming or not conforming mesh, with domain decomposition or not, see e.g. [START_REF] Min | A second order accurate level set method on non-graded adaptive cartesian grids[END_REF][START_REF] Losasso | Simulating water and smoke with an octree data structure[END_REF][START_REF] Fuster | Numerical simulation of droplets, bubbles and waves: state of the art[END_REF][START_REF] Coupez | Solution of high-Reynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing[END_REF].

In this work, applications with complex geometries are aimed, or even later couplings of models on hybrid mesh. This is why, as in a previous work [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF], a non conforming unstructured mesh is used. Therefore, although computing time may be lost compared to fully structured Cartesian codes, Block-Based Adaptive Mesh Refinement strategy is adopted. In accordance with the simulations presented in this paper, only the bidimensional case is presented. The strategy adopted is then as follows and illustrated in the Fig. 2. The unstructured mesh is composed of quadrilateral and/or triangular elements, where each element, as in [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF], defines a root element or block, see Fig. 2(a). Then, a mesh refinement level is defined for each block, initially by the user and thereafter in accordance with the chosen mesh refinement criterion, see Fig. 2(b). Coarsening and refinement are decided according to threshold values 0 < β c < β r . If η denotes the block-level error estimate, three situations are encountered:

For β c < η < β r , the block remains unchanged and so for the related elements; For β r < η, the block is refined so its elements are split into four isotropic subelements; For η < β c , the block is coarsened so four neighbouring subelements are merged into one element.

Limitations of such selection are described in [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF] as well as a procedure to overcome them. The mesh refinement level is adapted according to the rule that the ratio of mesh refinement should not exceed 2 between two neighbouring blocks. As illustrated in [START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF], this constraint allows smooth transitions between refined and unrefined regions. The mesh is then built using a quadtree graph (available for quadrangle and triangle) in order to define easily and precisely the neighbouring elements of each face. This procedure is pursued until the desired level of mesh refinement is reached as shown in Fig. 2(c). In the case of evolutionary problems, this implies frequent remeshing to follow the phenomenon studied. This is why we prefer to widen the stencil with a coarse root mesh, namely blocks, and therefore to use the whole quadtree graph, in order to remesh less often. It is illustrated for quadrangles and triangles in Fig. 3. Finally, a space filing curve using Morton numbering (Z-order) is built in order to number easily the degrees of freedom line in Fig. 2(d).

During a simulation, if a mesh refinement or coarsening process takes place, a new evaluation of the unknown field must be carried out. The prolongation (refinement) and restriction (coarsening) process are facilitated by the fact that the quadtree graph connect explicitly the "mother" element to its four "daughter" elements. To do that, one may use local interpolation/extrapolation. This is not straightforward in practice because it requires to identify the number of points and their position on the element to construct a right-determined system giving a fairly close solution. In the present study, the projection is made by solving a local DG problem. This technique gives good degrees of freedom by performing calculations at quadrature points. The element-wise weak formulation corresponding to h Enew = h E old reads:

Enew h Enew v dE new = Enew h E old v dE old . (31)
Then, mass matrix can be used to project the solution on the new mesh. This method is also used to enforce the initial condition h = h 0 .

Figure 3: Mesh refinement for quadrangle and triangle using quadtree.

A posteriori error estimation

An a posteriori error indicator is employed to show how a simple quantity linked to problem physics is able to provide information about refinement. This error indicator is inspired from Miller et al. [START_REF] Miller | A spatially and temporally adaptive solution of Richards' equation[END_REF] where they use changes in effective saturation. Here, the error indicator is based directly on changes in the solution that is to say the gradient of hydraulic head h. This allows to work both for the saturated and unsaturated regions. So, it reads:

C E = 1 max E∈E h 1 |E| ∇h L 2 (E) 1 |E| ∇h L 2 (E) . (32)
The error indicator is normalized so that 0 ≤ C E ≤ 1. Another a posteriori error indicator is introduced but it is based on error estimation. Indeed, Verfürth developed a posteriori estimation-based error indicators based on residual and derived thanks to suitable norms in the context of finite elements methods for many equations [START_REF] Verfürth | A Posteriori Error Estimation Techniques for Finite Element Methods[END_REF], in particular for nonlinear parabolic problems such as Richards' equation. For this kind of problem, Verfürth used some variants to prove reliability as well as global and local efficiency for implicit Runge-Kutta schemes with some restrictions [START_REF] Verfürth | A posteriori error estimates for nonlinear problems: L r (0, T ; W 1,ρ (Ω))error estimates for finite element discretizations of parabolic equations[END_REF] and more standard time discretization in [START_REF] Verfürth | A posteriori error estimates for finite element discretizations of the heat equation[END_REF][START_REF] Verfürth | A posteriori error estimatesfor non-linear parabolic equations[END_REF]. Later on, this type of residual-based derivation was adapted by Melenk and Wohlmuth to hp-adaptive finite elements methods for elliptic problems [START_REF] Melenk | On residual-based a posteriori error estimation in hp-FEM[END_REF]. Afterwards, Houston et al. [START_REF] Houston | Energy Norm A Posteriori Error Estimation of hp-Adaptive Discontinuous Galerkin Methods for Elliptic Problems[END_REF] as well as Schötzau and Zhu [START_REF] Schötzau | A robust a-posteriori error estimator for discontinuous Galerkin methods for convection-diffusion equations[END_REF][START_REF] Zhu | A robust a posteriori error estimate for hp-adaptive DG methods for convection-diffusion equations[END_REF] developed residual-based a posteriori error estimation for hp-adaptive interior penalty DG methods applied to elliptic problems and convection-diffusion equations respectively.

Greatly inspired by these works and supported by an heuristic analysis, the following a posteriori estimation-based error indicator is used:

(η n E ) 2 = η n E,R 2 + η n E,F 2 + η n E,J 2 , ( 33
)
where η n E,R , η n E,F and η n E,J are respectively the element residual, the face residual and the solution jump sub-estimates. I F and D F were defined in Eq. [START_REF] Dolejší | Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF]. The quantities λ m (K) and λ M (K)

stand respectively for the lowest and largest eigenvalue of K on E. For a face F of normal n, the minimum of normal component from neighbouring K is chosen: κ m := min(κ l , κ r ) where

κ = n • K u n+1 h
• n. So sub-estimates write:

η n E,R 2 = h 2 E p 2 E λ m (K) θ u n+1 h -θ u n h τ n -∇ • K u n+1 h ∇u n+1 h 2 L 2 (E) , (34)
η n E,F 2 = F ∈∂E∩F I h h F 2p F κ m K u n+1 h ∇u n+1 h • n 2 L 2 (F ) + F ∈∂E∩F N h h F p F κ l q N -K u n+1 h ∇u n+1 h • n 2 L 2 (F ) , (35)
η n E,J 2 = F ∈∂E∩F I h 1 2 I F + h F p F κ m u n+1 h 2 L 2 (F ) + F ∈∂E∩F D h D F + h F p F κ l u D -u n+1 h 2 L 2 (F ) . (36)
Here, the time discretization is based on the implicit Euler scheme but the a posteriori estimationbased error indicator can be determined for other time schemes such as in [START_REF] Verfürth | A posteriori error estimates for nonlinear problems: L r (0, T ; W 1,ρ (Ω))error estimates for finite element discretizations of parabolic equations[END_REF]. Global estimates are obtained with:

(η n ) 2 = E∈E n h (η n E ) 2 . (37)
Unlike the aforementioned papers, there is no rigorous mathematical proof for this residualbased energy norm a posteriori error estimation. Extension to the case of nonlinear and timedependant parabolic equation solved by hp-adaptive DG methods remains an open problem. However, some remarks can be drawn. Firstly, the numerical treatment of the time dependence and of the nonlinearity resorts to consider a sequence of linear and steady problems. Secondly, as noted by Verfürth [START_REF] Verfürth | A Posteriori Error Estimation Techniques for Finite Element Methods[END_REF], the element residual term is related to the residual of the numerical solution with respect to the strong form of the equation which may be viewed as the error from the nonlinear process. The face residual term is related to the boundary operator which is associated with the strong and weak forms of the differential equation. It reflects, on one hand, that the numerical solution gradient -the flux -is discontinuous and on the other hand, that Neumann boundary conditions may not be satisfied. The solution jump term is related to the penalization which are associated with the weak form of the equation. It reflects that the numerical solution is discontinuous at interior faces in the DG framework and Dirichlet boundary conditions may be slightly violated according to the penalty coefficients. Moreover, the second term and third parts of the estimation-based error indicators measure how valid the hypothesis on the seek solution are to derive the DG formulation.

Weighted discontinuous Galerkin framework

In the context of convection-diffusion equations or coupling of parabolic/hyperbolic domains, the solution may hold spurious oscillations at internal layers for vanishing or varying diffusion [START_REF] Ern | Multi-algorithmic methods for coupled hyperbolicparabolic problems[END_REF]. This situation is typical of porous media problems where there are material heterogeneity and degenerate hydraulic properties. One interpretation is the flow stays continuous but the solution mimics a discontinuity. Then, the numerical scheme is unable to recognize the sharp internal layer leading to numerical instabilities.

Even though this phenomenon may be resolved by mesh refinement, the corresponding computational cost is substantial to maintain a physically acceptable solution. Slope limiters are another possibility to control these undershoots/overshoots but further computational developments are needed to cope with the geometry of elements in high-dimension [START_REF] Sochala | Méthodes numériques pour les écoulements souterrains et couplage avec le ruissellement[END_REF][START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF]. A satisfactory approach would be to design a DG method that can handle internal layers in an automated fashion. That is the reason why the weighted discontinuous Galerkin (WDG) framework was introduced [START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF][START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF]. These methods work thanks to two key ingredients which can be seen as a way to incorporate into the definition of the scheme some partial knowledge of the solution. The first one is the use of weighted averages instead of the standard arithmetic average inside the discontinuous formulation. Then, the amount of diffusivity flux for each side of a face is controlled. The second ingredient is to relax the face penalization used for continuity constraint by a coefficient depending on the harmonic mean of the diffusivity of the neighbouring elements. Such a penalty strategy turns out to tune automatically the amount of local penalty to regulate the degree of smoothness of the approximate solution. Analysis and details can be found in [START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF][START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF] and related.

Internal layers depend on the spectral structure of the diffusivity. The choice is to take its normal component for face evaluation [START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF][START_REF] Burman | A Domain Decomposition Method Based on Weighted Interior Penalties for Advection-Diffusion-Reaction Problems[END_REF]. Weights across a face are positive numbers such that ω l + ω r := 1. The weighted average and the conjugate weighted average are respectively defined as:

⦃u⦄ ω := ω l u l + ω r u r , ⦃u⦄ ω := ω r u l + ω l u r . (38)
Then, the relation ∀u, v ∈ S p E n h , uv = ⦃u⦄ ω v + u ⦃v⦄ ω is replacing Section 3.1 in the weak formulation derivation. The weighted diffusion penalty coefficient is taken to be the harmonic mean of the neighbouring normal hydraulic conductivities across the face [START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF]: [START_REF] Li | Adaptive local discontinuous Galerkin approximation to Richards' equation[END_REF]

γ F = 2κ l κ r κ l + κ r . ( 40
)
Standard DG methods assume ω l = ω r = 1 2 , which reduces the weighted average to the arithmetic average, and γ F = 1. The WDG method chooses other definitions. Following definitions in [START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF], let the weights be:

   ω l = κ r κ l + κ r , ω r = κ l κ l + κ r if κ l + κ r = 0, ω l = ω r = 1 2 otherwise. (41)
With this definition choice for weights, γ F results to be equivalent to the weighted average of diffusivities:

γ F = 2κ l κ r κ l + κ r = ⦃κ⦄ ω . (42)
WDG methods assume that diffusivity discontinuities fit the mesh [START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF]. This is generally the case for heterogeneous media but it is not true for nonlinear diffusivity such as in Richards' equation. Then, sharp internal layers may occur inside an element so the weighted framework is not expected to work in this situation. Nevertheless, a suitable adaptive mesh refinement can be used as a capturing technique. This idea is one of the main novelties proposed in this paper. The weighted framework and mesh adaptation are working dynamically in synergy: the former changes the smoothness nature of the numerical solution while the latter tracks the internal layer linked to nonlinearity thanks to refinement driven by the a posteriori estimation-based error indicator Eq. ( 33).

Numerical results

Simulations are performed to highlight some numerical behaviour of Richards' equation and the ability of the numerical methods to treat them. The IIPG method is employed in every simulation thereafter.

Polmann's test-case

This test-case considers a soil from New Mexico whose hydraulic properties are provided by Polmann et al. [START_REF] Polmann | Application of Stochastic Methods to the Simulation of Large-scale Unsaturated Flow and Transport[END_REF]. A 1D vertical sample of this soil subject to downward infiltration was simulated by Celia et al. [START_REF] Celia | A general massconservative numerical solution for the unsaturated flow equation[END_REF], Manzini and Ferraris [START_REF] Manzini | Mass-conservative finite volume methods on 2-D unstructured grids for the Richards' equation[END_REF] as well as Sochala [START_REF] Sochala | Méthodes numériques pour les écoulements souterrains et couplage avec le ruissellement[END_REF]. This test-case can be challenging because, on one hand, the prescribed pressure head on both sides of the 100 cm soil column has a difference of 925 cm resulting in steep solicitation, and, on the other hand, For the mesh M100, nonlinear iterations have difficulties to reach the requirements of stopping criteria and the solution holds spurious oscillations, in particular an undershoot ahead of the sharp wetting front as showed in Fig. 4. This behaviour was already noticed by Celia et al. [START_REF] Celia | A general massconservative numerical solution for the unsaturated flow equation[END_REF] and Sochala [START_REF] Sochala | Méthodes numériques pour les écoulements souterrains et couplage avec le ruissellement[END_REF]. The latter decides to eliminate this undershoot by adding a slope limiter which causes a small late compared to the non-limited solution. Better results are obtained with the mesh M1000, the oscillations in the solution vanish, see Fig. 5. The M1000-results agree with those from Celia et al. [START_REF] Celia | A general massconservative numerical solution for the unsaturated flow equation[END_REF] as well as from Manzini and Ferraris [43]. This shows that a good discretization is necessary for Richards' equation to get quality solution. Such requirement can be fulfilled at optimal cost by a suitable adaptive local refinement. The error indicator from Eq. ( 32) is employed with refinement threshold values β c = β r = 50. Adaptation is performed every 5 time steps, starting with the M100 mesh. In Fig. 5, results show that refinement is able to capture the wetting front dynamically and eliminate the undershoot. The average number of elements is 211. More frequent adaptations allow to follow perfectly the front while less frequent adaptations are less costly.

Figure 5: Comparison of pressure head at t = 24 h for Polmann's test-case: IIPG on mesh with 100 elements (M100), IIPG on mesh with 1000 elements (M1000), IIPG with adaptive mesh refinement (AMR), WDG with 100 elements (WDG).

In addition, the WDG method from Section 4.3 is tried without any adaptation. The mesh M100 is used. Results in Fig. 5 show that the undershoot is still present but reduced to one single element. It is relevant to point out that the solution holds jumps where the wetting front is steep. They prevent oscillations to propagate and make solution catch the front steepness. The remaining undershoot depends on the wetting front localization in relation with the mesh geometry. This is due to the nonlinear nature of diffusivity and that is why it is important to associate mesh adaptation.

Tracy's benchmark

This test-case has an analytical solution given by Tracy [START_REF] Tracy | Clean two-and three-dimensional analytical solutions of Richards' equation for testing numerical solvers[END_REF] for 2D and 3D problems which can be used as a benchmark for Richards' equation [START_REF] Tracy | Three-dimensional analytical solutions of Richards' equation for a box-shaped soil sample with piecewise-constant head boundary conditions on the top[END_REF][START_REF] Tracy | Testing Computational Algorithms for Unsaturated Flow[END_REF]. Tracy's benchmark is particularly relevant because it is transient, has a simple parameter to vary nonlinearities for Richards' equation and holds differentiated steep regions. It can serve to check the nonlinear solver robustness, convergence properties and adaptive mesh refinement. Then, it is very convenient to assess the solution quality and compare codes: Šolín and Kuraz [START_REF] Solin | Solving the nonstationary Richards equation with adaptive hp-FEM[END_REF] and Dolejší et al. [START_REF] Dolejší | Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF] used it to evaluate the performance of their respective method.

Tracy's benchmark employs Gardner-Irmay relations with K s = 1.1 m•d -1 , θ s = 0.5, θ r = 0, α = 0.1 m -1 and m = 1. It is solved for pressure head ψ. The residual pressure head ψ r is a parameter. The computational domain Ω is a square (0, a) × (0, b) m. Here, a = b = 1 m and ψ r = -10 m. A specific Dirichlet boundary condition is prescribed on the top edge with the function:

ψ top = 1 α log e αψr + 1 -e αψr sin πx a . ( 43
)
While for other edges, a constant Dirichlet boundary condition is imposed: ψ = ψ r . The initial condition is ψ 0 = ψ r . Tracy gives the 2D exact solution ψ ex in [START_REF] Tracy | Clean two-and three-dimensional analytical solutions of Richards' equation for testing numerical solvers[END_REF]. Figure 6 shows the solution obtained for a fixed grid with τ = 10 -4 day and T = 10 -2 day. The treatment of initial condition is not straightforward. First, because the initial condition does not respect the top boundary condition: they are inconsistent. This issue can be troublesome for discretizations like finite elements methods [START_REF] Solin | Solving the nonstationary Richards equation with adaptive hp-FEM[END_REF] but it is natural for DG methods where boundary conditions are usually enforced by weak penalization [START_REF] Dolejší | Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF]. Second, because the top boundary condition condition is particularly stiff compared to the initial state. Therefore, the very first moments are the most interesting to simulate in order to address the time-boundary inconsistency and the early stiffness. Nonlinearities are controlled by the parameter α in Eq. ( 43). More α is large, more the relative hydraulic conductivity drops quickly to near zero, and so, nonlinearities are increased for Richards' equation [START_REF] Tracy | Three-dimensional analytical solutions of Richards' equation for a box-shaped soil sample with piecewise-constant head boundary conditions on the top[END_REF]. Moreover, errors are linked to the top corners since the top boundary condition induces a steep slope at these corners [START_REF] Tracy | Three-dimensional analytical solutions of Richards' equation for a box-shaped soil sample with piecewise-constant head boundary conditions on the top[END_REF][START_REF] Tracy | Testing Computational Algorithms for Unsaturated Flow[END_REF]. As Tracy stated, these features make this problem a good option to investigate adaptive mesh refinement.

The classic L 2 -norm and, following Rivière [START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF], Schötzau and Zhu [START_REF] Schötzau | A robust a-posteriori error estimator for discontinuous Galerkin methods for convection-diffusion equations[END_REF][START_REF] Zhu | A robust a posteriori error estimate for hp-adaptive DG methods for convection-diffusion equations[END_REF] or Dolejší et al. [START_REF] Dolejší | hp$-Adaptation Driven by Polynomial-Degree-Robust A Posteriori Error Estimates for Elliptic Problems[END_REF], the energy norm, also called the DG intrinsic norm, are introduced for a scalar function u : R d -→ R:

u 2 L 2 (E) := E u 2 dE |||u||| 2 E(E) := u 2 R(E) + F ∈∂E u 2 J(F ) , (44)
where the residual seminorm and the jump seminorm are respectively:

u 2 R(E) := E K(u)∇u • ∇u dF, ( 45
)
u 2 J(F ) := F F u 2 dF. (46)
Errors between numerical and exact solution are computed for a time step such that:

e X = ψ ex -ψ X with X = L 2 , E, R, J . (47)
To measure how effective the estimation-based error indicator is, the effectivity indices for the estimation-based error indicator [START_REF] Irmay | On the hydraulic conductivity of unsaturated soils[END_REF], the element and face residuals [START_REF] Jäger | Solution of doubly nonlinear and degenerate parabolic problems by relaxation schemes[END_REF][START_REF] Klieber | Adaptive simulations of two-phase flow by discontinuous Galerkin methods[END_REF], and the solution jump (36) will be used:

I eff E := η n E |||e||| E(E) , I eff R := η n E,R 2 + η n E,F 2
1 2 e R(E) , I eff J := η n E,J e J(F ) . (48)
Due to the lack of rigorous derivation, the effectivity index is not expected to represent properly the true error. The estimator efficiency is rather appreciated to set up a mesh refinement strategy.

First of all, several simulations are performed to show convergence properties. Because of symmetry, one half vertical plan of the domain is considered to speed up the simulation. Duration is set to T = 10 -4 day to focus on small times. The computations are carried out with the 6-step BDF and time steps are set to τ = 10 -6 s. This is needed to reduce time discretization errors expecting that they will be negligible compared to space discretization errors in order to observe space convergence. Indeed, according to theoretical error estimates [START_REF] Dolejší | Discontinuous Galerkin Method[END_REF], error in L 2 -norm and H 1 -seminorm behave such that:

e L 2 ≈ c h h p + c τ τ q , ( 49
)
e E ≈ c h h p+1 + c τ τ q , ( 50
)
where c h , c h , c τ and c τ are constants independent of h and τ . Numerical experiments involve four meshes of 50, 200, 400 and 3200 quadrilateral elements. Order approximation varies from one to four. Figure 7 shows the convergence of the error e h for the L 2 -and energy norms with respect to the number of degrees of freedom. It is interesting to note that space error are so small for the most dense mesh with p = 4 that time error saturates the convergence. In Fig. 8, effectivity indices are plotted. The error measured with the jump seminorm is independent from the exact solution and so is equivalent to the jump indicator [START_REF] Kumar | Formal upscaling and numerical validation of unsaturated flow models in fractured porous media[END_REF]. Then, the scaling factor of the jump indicator was chosen to be the penalty weight in order to be equal to jump seminorm of the error. As a consequence, I eff J controls only the numerical representation of exact solution (interpolation error) and must be equal to one as observed in Fig. 8e. On the contrary, I eff E and I eff R are not representative of the true error in term of magnitude like observed on Fig. 8a and8c. They are blowing up, at least for P1 and P2. It is probably due to a wrong choice for the scaling factors. The factors come from from Schötzau and Zhu [START_REF] Schötzau | A robust a-posteriori error estimator for discontinuous Galerkin methods for convection-diffusion equations[END_REF][START_REF] Zhu | A robust a posteriori error estimate for hp-adaptive DG methods for convection-diffusion equations[END_REF] who considered another problem. They were kept because they scale properly the sub-indicators to drive mesh adaptation as a capturing technique like showed later on. Moreover, the scaling factors always overestimate the true error which is good for mesh refinement monitoring. True a posteriori error analysis for the nonlinear problem of Richards' equation is known to be difficult and remains completely open.

The behaviour of the three contributions with respect to the number degrees of freedom is shown in Fig. 8 by varying the order of approximation and the mesh size. Since they do not represent properly the true error, observations are difficult to relate with the true error. However, some remarks can be drawn. Firstly, these diagrams show the convergence up to a certain extent. Secondly, one can see that the volume contribution -the residual indicator from Eq. ( 34) -behaves quite differently compared to the surface contributions -the flux and jump indicators from Eq. ( 34) and [START_REF] Klieber | Adaptive simulations of two-phase flow by discontinuous Galerkin methods[END_REF]. The convergence is globally maximised with p-refinement for the surface contributions while h-refinement seems to improve convergence better for the volume contribution. This observation should be assessed in more details but could drive the hp-decision making. Error estimation is evaluated in terms of ability to drive the mesh adaptation. The error in L 2 -norm is above the error measured with energy norm as stated in Fig. 9. It is worth to notice that the estimation-based error indicator Error estimation is evaluated in terms of ability to drive the mesh adaptation. The distribution of e L 2 (E) with no mesh adaptation is given in Fig. 9 for reference to localize the zones of interest. It is above the error measured with energy norm. In Fig. 10, it is worth to notice that the estimation-based error indicator η n E is able to give the shape (but not the magnitude) of the true error measured with the energy norm, particularly for steady state, even though it was not expected because there is no theoretical proof. Besides, during transient state, the estimation-based error indicator spots the error in the top corner where Tracy stated there were problems for this test-case, see Fig. 9 at t = 10 -4 d. The gradient-based error indicator C E was used successfully in a previous paper [START_REF] Clément | Adaptive Discontinuous Galerkin Method for Richards Equation[END_REF] to capture the moving wetting front on a simple 1D infiltration case. Yet, as shown in Fig. 11

Wetting of the La Verne dam

Setting

The model is here used to simulate the full-scale wetting of a multi-materials dam: the La Verne dam. La Verne dam was constructed in 1991 to supply Gulf of Saint-Tropez (south of France) with drinkable water. La Verne dam is 42 m high and peaks at 90 m above sea level. Its width ranges from -102.3 m to 110.5 m. The reservoir is on the left in Fig. 13.

The genuine La Verne dam inner section shows a complicated heterogeneous material patterns, partially described in [START_REF] Fleureau | Validation des modèles de couplage sur ouvrages réels[END_REF][START_REF] Bonelli | Résolution auto-adaptative par éléments finis de problèmes de diffusion fortement non linéaires[END_REF]. The dam is an earth-filled embankment dam. A clay core allows for impermeability while outer zones are semi-permeable alluvial shells. The dam is protected by rip-rap upstream. Filters and drains of very permeable gravel-like materials secure the dam inner saturation. The dam foundation lies on an impervious rock stratum thanks to concrete injections. A fine loam layer covers the downstream dam part. The filling of La Verne dam was controlled with instruments [START_REF] Fleureau | Validation des modèles de couplage sur ouvrages réels[END_REF] during 40 days, which provides a unique and challenging dataset to test the model. In particular, the reservoir height evolution is available as well as three groundwater hydraulic head measurements from the core, see Fig. 13 For the simulation, the reservoir height is used as a forced boundary condition in the form of a Dirichlet boundary condition monitored by a function based on experimental data (see Fig. 14). Since comparison focusses on the core of the dam, it is assumed that outer high permeable and/or fine layers materials do not have a significant impact on inner groundwater flows and saturation. The downstream thin loam layer, the filter-isolated toe drain and the lower submerged rip-rap are therefore not represented by the simulation and assimilated to the adjacent materials (see Tab. 4 for description of numerical dam structure). A Dirichlet boundary condition for hydraulic head is prescribed at the downstream boundary. The upper rip-rap and the cofferdam are simulated as such to assess the method robustness owing to their direct contact with the upstream dynamic boundary condition (forcing) and their hydraulic properties are steep. Rock/concrete foundations are supposed to be perfectly impervious which means a zero-valued Neumann boundary condition can be prescribed. Finally, seepage boundary condition are prescribed everywhere else because outflow can drain from the exposed upstream and downstream shell slopes. The actual simulated configuration of La Verne dam is sketched in Fig. 13. Five different materials are simulated (Tab. 4), with hydraulic properties provided by Bonelli et al. [START_REF] Bonelli | Résolution auto-adaptative par éléments finis de problèmes de diffusion fortement non linéaires[END_REF]. The dam wetting simulation involves each tool presented in this paper for the adaptive strategy: the WDG method in combination with the AMR technique driven by the estimation-based error indicator. Mesh adaptation is done every five time steps. Fixed-point iteration is used together with adaptive time stepping. Order of approximation is quadratic for space and one for time (implicit Euler) in order to have a robust and cost-effective computation.

Penalty parameters are σ I = σ B = 100. The simulation of La Verne dam case is a very challenging benchmark for numerical models, combining heterogeneous materials, steep constitutive laws and dynamic boundary conditions. However, despite of its interest, the La Verne dam benchmark remains only partly documented by in-situ instruments, precluding a comprehensive quantitative confrontation between observations and model. In particular, the hydraulic properties, which have a drastic effect on flow dynamics, remain approximative due to the absence of direct in-situ characterization. In addition, note that no initial in-situ values for water table elevation and saturation are available inside the dam, which can also strongly affect the subsequent evolution. For the present simulation, the initial water table is imposed at 67 m to fit the initial values of experimental hydraulic head. Finally, Bonelli et al. pointed out that hydromechanical coupling should be considered to be fully consistent with the real case, see [START_REF] Fleureau | Validation des modèles de couplage sur ouvrages réels[END_REF]. Figure 14: Time evolution of hydraulic head h, mesh and water table (white line) during La Verne dam wetting.

Results

Figure 14 depicts hydraulic head (color contours), water table (white line) and adaptive mesh for both simulations every 10 days. Time evolution of simulated hydraulic head and water content are displayed in Fig. 15 at the position of the three experimental sensors (S1 to S3) and for two additional relevant points called numerical sensors (S4 and S5). It is recalled that initial reservoir free surface and water table elevation are at 57 and 67 m, respectively. Note that the water table is not displayed in Fig. 14 at some points due to insufficient resolution and oscillations w.r.t the front sharpness, which are detailed in the Discussion section. The numerical model is observed to provide a good overall description of the dam wetting process. Figure 14 shows that the propagation of the wetting wave inside the dam is strongly dependent on the inner materials properties. The first eleven days are characterized by an adjustement phase during which the reservoir surface remains below the initial water table elevation. Hydraulic head and water content at numerical sensors S4 and S5 show a very constrasting behavior (Fig. 15b andc). This reveals that, even if both m2-and m3-materials show the same hydraulic conductivity at saturation, the difference in constitutive laws induces a radically different dynamics of the capillary fringe. At S4, the head remains nearly constant during the adjustement phase while the water content is very low. The m3-material is here nearly desaturated, well above a thin capillary fringe. By contrast, S5 is within the large capillary fringe associated to the m2material, providing a much higher water content. The hydraulic head shows first a small decay, probably due to the drain proximity, before starting to rise in response to the reservoir filling. Around day 14, both S4 and S5 sensors show a simultaneous regular rising, slightly lagged from the reservoir level. A strong jump is observed at S4 (Fig. 15c), indicating the crossing of a sharp capillary fringe. While rising, the water table remains nearly horizontal in the upstream m2and m3-materials zone but a sharp front is observed to develop across the left inner drain (m5material), see Fig. 14. Small numerical oscillations in hydraulic head are observed in Fig. 15 but do not affect the global dynamics. Further downstream and lower within the dam, the hydraulic head simulated at S1, S2 and S3 shows the slow wetting wave propagation, with attenuation and damping (Fig. 15b). At S1, S2, S3, the water content remains constant during the whole simulation, in accordance to the sensors position below water table.

When compared to hydraulic head experimental data (see Fig. 15a), the overall time evolution and order of magnitude are well represented. Some discrepancies are observed, in particular the late increase of head which is experimentally observed on each sensor around 30 days, which tends to be underestimated by the model. However, recalling the lack of experimental control evoked before, more detailed comparison should be made with caution. For the sake of comparison, the previous numerical results presented by Bonelli et al. [START_REF] Bonelli | Résolution auto-adaptative par éléments finis de problèmes de diffusion fortement non linéaires[END_REF] are recalled in Fig. 15 in green lines. While the overall time evolution is rather similar between models, a better agreement with experimental data is obtained with the present simulation.

Discussion

505

The proposed modeling strategy for Richard's equation showed very satisfactory performance in reproducing analytical and laboratory test cases and in simulating complex full-scale exper-iment. In this latter case, the model is able to capture automatically wetting fronts which are moving dynamically thanks to AMR. In particular, the fronts associated with heterogeneity are well resolved by the WDG technique. This is illustrated by water content distribution in Fig. 16a. Moreover, the simulation is robust since it is able to compute completely without user intervention and handle many features: degeneracy for saturated/unsaturated zones, dynamic forcing boundary condition and seepage. However, the simulation struggles to capture wetting fronts associated with nonlinear varying diffusivity resulting in severe overshoots/undershoots for the transition between saturated and unsaturated zones in very permeable materials. Further insight on this behaviour can be discussed using the Péclet number framework. Richards' equation is an elliptic-parabolic equation predominantly diffusive, it shares properties with convection-diffusion equation and can behave like them. For example, sharp wetting fronts and internal layers can be reinterpreted within this context. The Péclet number represents the ratio of the rate of advection to the rate of diffusion. This dimensionless number is useful to analyse locally the nature of the flow regime and determine whether advection or diffusion is dominant. For this study, the Péclet number is extracted from the saturation-based formulation of Richards' equation which is similar to a convection-diffusion equation with the water content θ as single unknown. This formulation is undefined in saturated zones because hydraulic diffusivity D tending to infinity but this zone is not of interest for the Péclet number since Richards' equation is purely diffusive. Moreover, it is assumed that θ = 0. This assumption is not restrictive because if water content drops to zero, nothing happens and, in practice, there is often a residual water content. The total flux can be separated into a diffusive flux defined by hydraulic diffusivity and driven by capillarity, and an advective flux defined by hydraulic conductivity/water content and driven by gravity. The saturation-based formulation can be rewritten as: In order to overcome these limitations, an additional ad-hoc configuration, named augmented simulation, has been implemented. This one does not use WDG framework and other difference lies in the mesh and tolerances for nonlinear convergence criterion. Finer discretization is allowed in permeable material that is to say the m1-, m2-and m3-materials. The mesh is also refined around water table and both for the gradient-based and estimation-based error indicators. Instead of the reference simulation, 4-order BDF is employed and hydraulic properties are relaxed with a numerical minimal value to avoid complete degeneracy. The augmented computation allowed to eliminate spurious oscillations from the solution thanks to a more refined mesh, which explains why WDG method is not needed. As the mesh is extremely refined, the augmented simulation takes approximately 13.5 times longer to compute (42 h 37 min) than the reference simulation (3 h 11 min). One core Intel(R) Xeon(R) CPU 5-2630 v3 with 2.40 GHz was used. Additional information and results can be found in [START_REF] Clément | Numerical simulation of flows in unsaturated porous media by an adaptive discontinuous Galerkin method: application to sandy beaches[END_REF].

∂ t θ -∇ • D(θ)∇θ + K(θ)∇z θ θ = 0, (51)
The simulation of La Verne dam filling shows many difficult points which lead to two unwanted effects. The first one is the development of non-physical oscillations around steep wetting fronts associated with nonlinear varying diffusivity which affect solution accuracy. The other one is nonlinear convergence possible only with very small time steps so that the calculation duration is high. These challenges illustrate numerical issues of Richards' equation for heterogeneous media with high-varying diffusivity. The improvements offered by the augmented simulation support the idea that a series of heuristic tools (threshold values, maximum refinement level, refinement frequency, adaptive time steps) should be investigated to improve accuracy while maintaining cost-effectiveness.

Conclusion

In this work, a new strategy has been formulated for the solution of Richards' equation on the basis of discontinuous Galerkin method and adaptive mesh refinement. The approximation can reach high-order both in space and time. Issues outlined in the introduction have been tackled and further discussed throughout the paper so the following concluding remarks can be drawn:

The seepage boundary condition is directly incorporated into the DG weak formulation which makes its treatment natural and simple.

The adaptive time stepping allows the nonlinear iterative solver to converge giving robustness.

Mesh adaptation is employed to monitor spatial errors of the Richards' equation by capturing the moving wetting fronts thanks to an a posteriori error indicator.

The adaptive mesh refinement is based on a block structure to surround large regions where wetting fronts move aiming to save computational time by avoiding systematic refinement. It also prepares a future parallelization.

Oscillations in sharp wetting fronts of Richards' equation have been interpreted within the context of convection-diffusion equation. Wetting fronts are considered as internal layers due to the nonlinear varying, and possibly vanishing, diffusivity as well as material heterogeneity and anisotropy. A local Péclet number is exhibited to highlight troublesome regions. A weighted discontinuous Galerkin method is used to allow for jumps in the solution which better approximate the sharp internal layers. Mesh adaptation and the weighted framework work in synergy to capture and resolve sharp wetting fronts through the proposed a posteriori error indicator.

The strategy leads to mass conservative, efficient and robust solution of Richards' equation.

However, the method holds heuristic parameters so that questions remain open regarding the numerical analysis, the hp-decision making, the nonlinear convergence or the error balance. All these aspects have great potential of improvements and will be investigated in future studies to optimize the solving of Richards' equation.

A

  + |ψ| BA, C: empirical shape parameters [L B;D ] B, D: empirical shape parameters[-]

Figure 1 :

 1 Figure 1: Seepage modelling by the boundary condition.

Figure 2 :

 2 Figure 2: AMR Strategy: (a) Block mesh; (b) Level of mesh refinement; (c) mesh generation; (d) Morton numbering.

  hydraulic conductivity shows strong variations under the set of values taken by pressure head. The Polmann's test-case employs Van Genuchten-Mualem relations with K s = 9.22 × 10 -3 cm•s -1 , θ s = 0.368, θ r = 0.102, α = 3, 35 × 10 -2 cm -1 , n = 2, m = 0.5 and l = 0.5. The computational domain Ω is a rectangle (0, 20) × (0, 100) cm. The test-case is solved for pressure head ψ during T = 172800 s (48 h) with a constant time step τ = 120 s and the initial condition is ψ 0 = -1000 cm. Throughout the examples, implicit Euler scheme, p = 1, σ I = σ B = 100 are used. Computation is done for two meshes: M100 is a coarse mesh of 100 elements and M1000 is a refined mesh of 1000 elements.

Figure 4 :

 4 Figure 4: Pressure head along the vertical for Polmann's test-case with the mesh M100.

Figure 6 :

 6 Figure 6: Hydraulic head for Tracy's benchmark with p = 1 and 6-step BDF.

( a )

 a Errors in the L 2 -norm (b) Errors in the energy norm

Figure 7 :

 7 Figure 7: Numerical convergence for Tracy's benchmark. P1, P2, P3 and P4 stands respectively for linear, quadratic, cubic and quartic approximations.

Figure 8 :

 8 Figure 8: Effectivity indices for different indicators (left) and convergence diagrams for the three parts of the indicator η n E (right). P1, P2, P3 and P4 stands respectively for linear, quadratic, cubic and quartic approximations.

  and 12, it performs badly for the Tracy's benchmark where steep gradients arise on a large part of the domain. The refinement thresholds are difficult to tune which leads to over-refinement for large regions at the top. They were set up to β c = β r = 0.5 in order to give e L 2 (Ω) ≈ 0.167 at t = 10 -2 d. Then, the computation is very costly since it takes around 227 min. For the estimation-based error indicator, the ability to drive the mesh adaptation is more convincing. β c = β r = 0.85 in order to target e L 2 (Ω) = 0.137 at t = 10 -2 d. The mesh is adapted dynamically without inducing over-refinement. Elements are refined where it is needed because adaptation occurs locally while the error is controlled. Then, the cost is reasonable with a computation around 35 min.

Figure 9 :

 9 Figure 9: True error distribution in the L 2 -norm on Tracy's benchmark.

( a )

 a Distribution at t = 10 -4 days (b) Distribution at t = 10 -2 days

Figure 10 :

 10 Figure 10: Comparison between distributions of true error in energy norm with the indicator η n E on Tracy's benchmark.

( a )

 a Mesh and error distribution in L 2 -norm for CE-driven AMR (b) Mesh and error distribution in L 2 -norm for η n E -driven AMR

Figure 11 :

 11 Figure 11: Comparison with L 2 -norm between C E -based adaptation (left) and η n E -based adaptation on Tracy's benchmark.

Figure 12 :

 12 Figure 12: Comparison with energy norm between C E -based adaptation (left) and η n E -based adaptation on Tracy's benchmark.

  and 15. The three sensors are at z = 55 m and their x-coordinates are -10.5 m, -3 m and 5 m.

Figure 13 : 9 m4 64 A 64 A

 1396464 Figure 13: Geometry, materials and boundary conditions of La Verne dam for the numerical case. S1, S2, S3 are numerical sensors corresponding to the experimental sensors while S4 and S5 are additional purely numerical sensors.

Figure 15 :

 15 Figure 15: Time evolution of inner dam hydraulic parameters. a) Hydraulic head at sensors S1 (short-dashed lines), S2 (dash-dotted lines) and S3 (long-dashed lines) from the present simulation (red), the experimental data (light blue) and the previous results of Bonneli et al. [7]. The reservoir height is depicted in dark solid line. b) Hydraulic head at numerical sensors S4 (orange) and S5 (yellow). c) Water content at S1 to S5 sensors (note that S1, S2 and S3 show similar values).

  with the hydraulic diffusivity D(θ) = K(θ) dψ(θ) dθ [L 2• T -1 ]. The Péclet number can therefore be defined as: capacity C = dθ(ψ) dψ [L -1 ] and L the characteristic linear dimension assimilated to an element size h E . High Péclet numbers are clearly reached in the sharp-gradient areas in Fig. 16b, with m1-, m2-and m3-materials (i.e. drains, rip-rap and upstream shell) showing relatively high Péclet numbers, from 1 to more than 5. Consequently, water table is

510

  not straight anymore and a delay in elevation is observed upstream, see Fig.14.

( a )

 a Distribution of water content at t = 40 days (b) Distribution of Péclet number at t = 40 days

Figure 16 :

 16 Figure 16: Heterogeneity and nonlinear varying diffusivity highlighted respectively by water content and Péclet number for La Verne dam simulation with the reference simulation.

Table 1 :

 1 Hydraulic relations used in this paper.

	It is worth noting that:
	under complete saturation, hydraulic properties become constant and Richards' equation

degenerates into an elliptic equation characterised by fast diffusion; under almost complete unsaturation, hydraulic properties get very near-zero values which stop diffusion and may be inconvenient numerically;

  is multiplied by a test function v ∈ S p E n

	and then integrated on each h . Green's theorem is used element-wise before summing over all elements in h element E ∈ E n E n h . The Neumann boundary condition has appeared naturally within the formulation and has
	been substituted. So it stands for the homogeneous Neumann part of the seepage boundary
	condition.
	On one hand, knowing ∀u, v ∈ S p E n h , uv = ⦃u⦄ v +⦃v⦄ u and assuming K(h -z)∇h • n F =
	0 because this quantity is smooth enough, product jump is decomposed and simplified where it
	is possible.

Table 4 :

 4 Materials of the numerical La Verne dam.

Acknowledgement

The first author is partially supported by a grant from Provence-Alpes-Côte d'Azur region, France.


Samia Moumade
email: samia.moumade@iae-aix.com

Aurélie Hemonnet-Goujot

Pierre Valette-Florence

Cognitive or Affective? A Dual Perspective of Consumer

Keywords: Brand Activism, Consumer-Brand Relationship, Brand Relationship Quality Track: Product and Brand Management

come   L'archive ouverte pluridisciplinaire