Ali Raeisi

Vincent Fontaine

Hamid-Reza Ghafouri

Anis Younes

Marwan Fahs
email: fahs@unistra.fr

A fully interior penalty discontinuous Galerkin method for variable density groundwater flow problems

Keywords: Variable density flow, Compatible algorithms, Advanced model, seawater intrusion, Natural convection in porous enclosures 1. Introduction 1

Discontinuous Galerkin (DG) methods due to their robustness properties, e.g. local conservation, low numerical dispersion, and well-capturing strong shocks and physical discontinuities, are well-suited for the simulation of Variable Density Flow (VDF) in porous media. This paper aims at introducing, in a unified format, the general class of Interior Penalty DG (IPDG) methods to solve the VDF equations. A combination of symmetric, non-symmetric and incomplete IPDG methods is used to discretize both head and concentration variables. Compatibility analysis is performed to prevent the loss of accuracy of the IPDG methods in simulations of coupled flow and transport equations. An accurate technique is used for time integration, based on a non-iterative procedure and adaptive time stepping with embedded error control. Several benchmarks are investigated to validate the proposed DG scheme and to examine its performance in simulating VDF problems. The new DG scheme reproduces better the experimental data than the conventional SEAWAT model. Its results are in excellent agreement with a recent semianalytical solution of the Henry problem, dealing with seawater intrusion under convection-dominating conditions.

The performance of the DG scheme is examined by simulating the challenging problem of natural convection in porous enclosure. The method is compared against a finite element solution obtained with COMSOL multi-physics.

The numerical experiments indicate clearly that high-order DG method is much more appropriate than standard conforming Galerkin method in simulating VDF problems while at the same time, guaranteeing a better precision and high-fidelity solutions. The proposed numerical method can be extended to 3D problems.

of VDF problems requires accurate methods as numerical artifacts can affect the results, especially in unstable cases. This is also why benchmarking VDF numerical models is a common issue which has been widely discussed in the literature [START_REF] Fahs | A new benchmark semi-analytical solution for density-driven flow in porous media[END_REF][START_REF] Prasad | Using quantitative indicators to evaluate results from variabledensity groundwater flow models[END_REF][START_REF] Shao | A 3-D Semianalytical Solution for Density-Driven Flow in Porous Media[END_REF][START_REF] Simmons | Variable density groundwater flow: From current challenges to future possibilities[END_REF][START_REF] Stoeckl | A new numerical benchmark of a freshwater lens[END_REF][START_REF] Voss | Three-dimensional benchmark for variable-density flow and transport simulation: matching semi-analytic stability modes for steady unstable convection in an inclined porous box[END_REF]. Clearly, in a numerical simulation, the level of the solution accuracy profoundly depends on the applied numerical scheme.

In the literature, several studies investigated the different aspects of VDF numerical modelling. For instance, Mazzia and Putti [START_REF] Mazzia | Mixed-finite element and finite volume discretization for heavy brine simulations in groundwater[END_REF] validated the formulations resulting from the combination of two sets of dependent variables, i.e. head/concentration and pressure/mass-fraction, and concluded that the latter one is more accurate especially in long simulation periods. Comparison of conservative and non-conservative formulations [START_REF] Oltean | Coupled Groundwater Flow and Transport in Porous Media. A Conservative or Non-conservative Form?[END_REF] and validity of Oberbeck-Boussinesq based formulation were also explored [START_REF] Guevara Morel | Systematic investigation of non-Boussinesq effects in variable-density groundwater flow simulations[END_REF][START_REF] Johannsen | On the Validity of the Boussinesq Approximation for the Elder Problem[END_REF][START_REF] Kolditz | Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models[END_REF].

Younes et al. [START_REF] Younes | Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods[END_REF] studied the efficiency of high-order time integration schemes to solve the partial differential equations of the VDF model. Different time-stepping techniques have been also investigated in Hirthe and Graf [START_REF] Hirthe | Non-iterative adaptive time-stepping scheme with temporal truncation error control for simulating variable-density flow[END_REF], Younes and Ackerer [START_REF] Younes | Empirical versus time stepping with embedded error control for densitydriven flow in porous media[END_REF]. Consistent velocity approximation is one of the numerical challenges addressed by Albets-Chico and Kassinos [START_REF] Albets-Chico | A consistent velocity approximation for variable-density flow and transport in porous media[END_REF], Diersch and Kolditz [START_REF] Diersch | Variable-density flow and transport in porous media: approaches and challenges[END_REF], Frolkovič [START_REF] Frolkovič | Consistent velocity approximation for density driven flow and transport[END_REF], Herbert et al. [START_REF] Herbert | Coupled groundwater flow and solute transport with fluid density strongly dependent upon concentration[END_REF], Voss and Souza [START_REF] Voss | Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater-saltwater transition zone[END_REF]. Another common challenge that the authors have encountered is the solution of the coupled nonlinear systems obtained from numerical discretization. For example, Putti and Paniconi [START_REF] Putti | Picard and Newton linearization for the coupled model for saltwater intrusion in aquifers[END_REF] presented a partial Newton method in order to reduce the size of matrix systems equal to that of Picard. A new coupling procedure was also proposed by Ackerer [START_REF] Ackerer | A new coupling algorithm for density-driven flow in porous media[END_REF] to improve the rate of convergence in Picard method.

However, one of the most crucial questions for the numerical solution of the VDF model is the discretization of the spatial derivatives. Numerous studies have been devoted to the various spatial discretization methods such as the finite difference method [START_REF] Croucher | The Henry Problem for Saltwater Intrusion[END_REF][START_REF] Guo | User's Guide to SEAWAT: A computer program for simulation of three-dimensional variable-density ground-water flow[END_REF], the finite element method [START_REF] Diersch | Finite element modelling of recirculating density-driven saltwater intrusion processes in groundwater[END_REF][START_REF] Huyakorn | Saltwater intrusion in aquifers: Development and testing of a three-dimensional finite element model[END_REF][START_REF] Segol | A Galerkin-finite element technique for calculating the transient position of the saltwater front[END_REF], the finite volume method [START_REF] Frolkovič | Numerical modelling of convection dominated transport coupled with density driven flow in porous media[END_REF][START_REF] Hamidi | Modeling of 2D density-dependent flow and transport in porous media using finite volume method[END_REF], and the method of characteristics [START_REF] Oude Essink | MOC3D adapted to simulate 3D density-dependent groundwater flow[END_REF][START_REF] Pinder | A Numerical Technique for Calculating the Transient Position of the Saltwater Front[END_REF]. However, using the conventional discretization methods (especially standard Galerkin finite element and finite difference methods) for the transport equation leads to numerical dispersion [START_REF] Miller | Numerical simulation of water resources problems: Models, methods, and trends[END_REF][START_REF] Werner | Seawater intrusion processes, investigation and management: Recent advances and future challenges[END_REF], and using nonconservative methods for the flow equation in a coupled system can result in erroneous solutions [START_REF] Diersch | FEFLOW: Finite element modeling of flow, mass and heat transport in porous and fractured media[END_REF]. DG methods were firstly introduced in the middle of the seventies for the numerical approximation of hyperbolic problems, and independently, in the context of elliptic and parabolic problems [see e.g., [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF][START_REF] Arnold | Discontinuous Galerkin Methods for Elliptic Problems[END_REF][START_REF] Cockburn | The Development of Discontinuous Galerkin Methods[END_REF]. For diffusion problems, Arnold [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF] introduced a primal DG method inspired by the original work of Nitsche (1971) using Interior Penalty (IP) technique to weakly enforce some regularity requirements of the solution across the skeleton of the mesh. This derivation yields to the Symmetric Interior Penalty (SIP) method and constitutes a milestone in the development of primal DG methods. In the late 1990s, numerous variations of the SIP method have been proposed and studied in the literature. For instance, Oden et al. [START_REF] Oden | A discontinuous hp finite element method for diffusion problems[END_REF] introduced the Oden-Babuška-Baumann (OBB) method for pure diffusion problems and then extended it to convection-diffusion processes [START_REF] Baumann | A discontinuous hp finite element method for convectiondiffusion problems[END_REF]. Compared to SIP, the differences lie in the use of consistency and penalty terms at interfaces that are now skew-symmetric and null, respectively. Alternatively, penalty terms have been added to the OBB formulation leading to the Non-symmetric Interior Penalty (NIP) method analyzed by Rivière et al. [START_REF] Rivière | Improved energy estimates for interior penalty, constrained and discontinuous Galerkin methods for elliptic problems[END_REF]. Finally, we mention the Incomplete Interior Penalty (IIP) method as introduced first by Dawson et al. [START_REF] Dawson | Compatible algorithms for coupled flow and transport[END_REF]. Thus, symmetric, non-symmetric and incomplete interior Penalty methods are the three most famous variations of the primal DG method, and they have been judiciously combined to solve efficiently coupled single-phase flow and reactive transport problems [see e.g., [START_REF] Dawson | Compatible algorithms for coupled flow and transport[END_REF], 80, and the references therein]. The DG methods are numerically stable even for high Péclet numbers, able to capture physical discontinuities well and to handle nonconforming and unstructured meshes [START_REF] Abbà | Dynamic models for Large Eddy Simulation of compressible flows with a high order DG method[END_REF][START_REF] Franciolini | Efficient discontinuous Galerkin implementations and preconditioners for implicit unsteady compressible flow simulations[END_REF]. They are locally adaptive in mesh and polynomial degrees and are well suited for efficient parallel implementations [see e.g. 14, 57, 79, and references therein].

Since its introduction, and due to their advantages, DG methods have benefited from intensive research and development, and they have been applied to a variety of physical issues and situations. In the context of porous media, the DG method has been used to discretize the hyperbolic terms of the governing equations in applications involving contaminant transfer in aquifers [START_REF] Younes | Solving the advection-dispersion equation with discontinuous Galerkin and multipoint flux approximation methods on unstructured meshes[END_REF], two phase flow nonfractured domains [START_REF] Fabien | A hybridizable discontinuous Galerkin method for two-phase flow in heterogeneous porous media[END_REF][START_REF] Fučík | Discontinous Galerkin and mixed-hybrid finite element approach to two-phase flow in heterogeneous porous media with different capillary pressures[END_REF][START_REF] Kummer | Extended discontinuous Galerkin methods for two-phase flows: the spatial discretization[END_REF][START_REF] Zidane | An implicit numerical model for multicomponent compressible two-phase flow in porous media[END_REF], two phase flow in fractured domains [START_REF] Hoteit | An efficient numerical model for incompressible two-phase flow in fractured media[END_REF][START_REF] Hoteit | Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures[END_REF][START_REF] Moortgat | Adaptive implicit finite element methods for multicomponent compressible flow in heterogeneous and fractured porous media[END_REF][START_REF] Moortgat | Implicit finite volume and discontinuous Galerkin methods for multicomponent flow in unstructured 3D fractured porous media[END_REF][START_REF] Zidane | An efficient numerical model for multicomponent compressible flow in fractured porous media[END_REF] and diffusion and natural convection in fractured domains [START_REF] Hoteit | Modeling of multicomponent diffusions and natural convection in unfractured and fractured media by discontinuous Galerkin and mixed methods[END_REF]. In this context, the combination of DG method and mixed finite element (MFE) method has received particular attention. The former has been used to discretize the convective term of the transport equation while the later has been used for flow. Sun et al. [START_REF] Sun | A Combined Mixed Finite Element and Discontinuous Galerkin Method for Miscible Displacement Problem in Porous Media[END_REF] used this combination with a cut-off operator in DG method to make the method converge. Li and Riviere [START_REF] Li | Numerical solutions of the incompressible miscible displacement equations in heterogeneous media[END_REF] developed a discretization method based on the combination of MFE-DG without any slope limiter into heterogeneous media which uses a high-order Runge-Kutta approximation for time.

The DG method has been also used to entirely solve the flow or transport (both hyperbolic and parabolic terms) equations in porous media. Li and Riviere [START_REF] Li | High order discontinuous Galerkin method for simulating miscible flooding in porous media[END_REF] used a weighted version of high-order IPDG methods for simulation of miscible displacement problems. Besides, Rivière and Wheeler [START_REF] Rivière | Discontinuous Galerkin methods for flow and transport problems in porous media[END_REF] described a full DG method with slope limiter for miscible displacement. For a coupled system of flow and reactive transport, a fully primal DG using a cut-off operator was developed by Sun and Wheeler [START_REF] Sun | Discontinuous Galerkin methods for coupled flow and reactive transport problems[END_REF]. Full primal DG method was used to simulate two-phase flow in Arbogast et al. [START_REF] Arbogast | A discontinuous Galerkin method for two-phase flow in a porous medium enforcing H(div) velocityand continuous capillary pressure[END_REF], Bastian [START_REF] Bastian | A fully-coupled discontinuous Galerkin method for two-phase flow in porous media with discontinuous capillary pressure[END_REF], Epshteyn and Rivière [START_REF] Epshteyn | On the solution of incompressible two-phase flow by a p-version discontinuous Galerkin method[END_REF][START_REF] Epshteyn | Fully implicit discontinuous finite element methods for two-phase flow[END_REF][START_REF] Epshteyn | Convergence of high order methods for miscible displacement[END_REF], Ern et al. [START_REF] Ern | Discontinuous Galerkin approximation of two-phase flows in heterogeneous porous media with discontinuous capillary pressures[END_REF], Jamei et al. [START_REF] Jamei | A Lax-Wendroff-IMPES scheme for a twophase flow in porous media using interior penalty discontinuous Galerkin method[END_REF], Kou and Sun [START_REF] Kou | Upwind discontinuous Galerkin methods with mass conservation of both phases for incompressible two-phase flow in porous media[END_REF], Mozolevski and Schuh [START_REF] Mozolevski | Numerical simulation of two-phase immiscible incompressible flows in heterogeneous porous media with capillary barriers[END_REF]. The three primal DG methods, namely IIP, SIP, and NIP methods were discussed in these works. Different DG methods have been also used to investigate variably-saturated flow expressed by Richards equation [START_REF] De Maet | A fully-explicit discontinuous Galerkin hydrodynamic model for variably-saturated porous media[END_REF][START_REF] Dolejší | Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF][START_REF] Li | Local discontinuous Galerkin approximations to Richards' equation[END_REF][START_REF] Li | Adaptive local discontinuous Galerkin approximation to Richards' equation[END_REF] and reactive transport in porous media [START_REF] Sun | Discontinuous Galerkin methods for coupled flow and reactive transport problems[END_REF][START_REF] Sun | Symmetric and Nonsymmetric Discontinuous Galerkin Methods for Reactive Transport in Porous Media[END_REF][START_REF] Sun | A posteriori error estimation and dynamic adaptivity for symmetric discontinuous Galerkin approximations of reactive transport problems[END_REF][START_REF] Sun | Analysis of Discontinuous Galerkin Methods for Multicomponent Reactive Transport Problems[END_REF][START_REF] Sun | Anisotropic and dynamic mesh adaptation for discontinuous Galerkin methods applied to reactive transport[END_REF][START_REF] Wheeler | Discontinuous Galerkin Method for Modeling Flow and Reactive Transport in Porous Media[END_REF]. This brief review shows an increasing interest in the development and utilization of the DG methods in several applications related to flow and mass transfer in fractured/unfractured porous media. The review paper by Miller et al. [START_REF] Miller | Numerical simulation of water resources problems: Models, methods, and trends[END_REF] estimated that, with the advancement of computing technology, DG methods have great potential to be foundational technique for future simulators of multi-physical processes in porous media and more generally for future solvers of physical partial differential equations. The dynamic adaptivity properties and localized nature of DG methods, especially, are well-adapted for many VDF situations where the dense physics/high spatial gradients occur only in a small space may vary continuously over time. Thereby, a comparable memory and processing time can be saved. However, despite their advantages, wide range of applications, increasing popularity and promising results, applications of DG methods in VDF are limited and performance of these methods in such a configuration has been never evaluated. For VDF, the DG methods have been only used for discretizing the convective terms of the transport equation and has been coupled with other methods (mainly MFE) for discretizing the flow and parabolic term of the transport equation [START_REF] Ackerer | Efficient approximations for the simulation of density driven flow in porous media[END_REF][START_REF] Ackerer | Modeling variable density flow and solute transport in porous medium: 1. Numerical model and verification[END_REF][START_REF] Buès | Numerical simulations for saltwater intrusion by the mixed hybrid finite element method and discontinuous finite element method[END_REF][START_REF] Younes | Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods[END_REF]. DG methods have never been used to solve the full VDF model while this could be very advantageous.

In a full DG model, one can efficiently solve pure diffusion, pure advection or any mixed situation with one and only one mathematical formalism, no splitting operator techniques, no need to introduce different variables for the diffusive and convective parts, and solution procedure can be efficiently accelerated using parallel computing. This avoids splitting errors and allows for treating mixed boundary conditions in accurate manner. Thus, the main aim of the present work is to show how DG methods can be used to solve the full VDF model and to evaluate the performance and benefits of this class of methods in simulating VDF problems. We develop a unified and coherent format of the general class of IPDG. We test symmetric, non-symmetric and incomplete IPDG methods to discretize flow and transport variables (pressure head and concentration). Flow and transport equations are solved sequentially. We implement an efficient scheme for time integration based on adaptive time stepping with error control. The new developed numerical scheme is validated against experimental data and semi-analytical solutions for problems involving VDF driven by salinity gradients. A highly convective thermal-driven case is considered to evaluate the performance of the developed DG scheme, by comparing it against standard finite element solution obtained with COMSOL Multiphysics.

The paper is organized as follows: In section 2, we describe the governing equations. Section 3 is devoted to the discretization of the governing equations using the IPDG method and to the numerical scheme implemented to solve the discretized equations. Section 4 aims at validating and verifying the new developed scheme by comparison against semi-analytical solutions (Henry and thermal porous-cavity problems) and experimental data. To highlight the performance of the developed DG scheme, a convergence analysis is presented in section 4. Finally, we end in section 5 with general conclusions.

Governing equations

The governing equations of variable density groundwater flow and solute transport are described following [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Huyakorn | Saltwater intrusion in aquifers: Development and testing of a three-dimensional finite element model[END_REF][START_REF] Kolditz | Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models[END_REF]. The set of partial differential equations include (i) the generalized Darcy' law, (ii) continuity equations of fluid and solute mass, and (iii) state equations for the bulk fluid density and/or viscosity. We will propose an equivalent compact formulation of both original problems based on the introduction of suitable variables. In the rest of the paper, we assume that previous coupled physical processes are given in a bounded domain Ω ⊂ R 2 with a Lipschitz boundary ∂Ω, and in time interval [0, T ], with T > 0.

Variable density flow equation

The generalized Darcy' law can be written in term of the reference head variable ϕ = p/(ρ 0 g) + z where p corresponds to the dynamic pressure [M L -1 T -2 ], ρ 0 the reference fluid density [M L -3 ], z the vertical elevation above some datum level [L], and g the gravity acceleration [LT -2 ]. Namely,

q = -K ∇ϕ + ρ -ρ 0 ρ 0 ∇z in Ω × (0, T ], (1)
where q denotes the Darcy's velocity [LT -1 ] and ρ is the fluid density [M L -3 ]. Here, K = ρ0gk µ is defined as hydraulic conductivity where k is the permeability tensor of the porous material [L 2 ], and µ the bulk

fluid viscosity [M L -1 T -1 ]
. Following Huyakorn et al. [START_REF] Huyakorn | Saltwater intrusion in aquifers: Development and testing of a three-dimensional finite element model[END_REF], the time-dependent mass balance equation of the fluid in a porous medium is given by:

ρs∂ t ϕ + φ∂ c ρ∂ t c + ∇• (ρq) = 0, in Ω × (0, T ], (2)
where s is the specific storativity of the porous medium related to head change given by:

ϕ = ϕ D , on ∂Ω D × (0, T ], (3a)
ρq • n = q N , on ∂Ω N × (0, T ], (3b)
ϕ (•, 0) = ϕ 0 , in Ω × {0}, ( 3c
)
where n is the outward unit normal vector on ∂Ω, ϕ D and q N are prescribed functions on ∂Ω D and ∂Ω N , respectively.

Solute transport equation

The solute transport process is governed by a time-dependent advection-dispersion equation,

∂ t (φρc) + ∇• (ρqc -ρD∇c) = 0, in Ω × (0, T ]. (4)
Here D corresponds to the hydrodynamic dispersion tensor [L 2 T -1 ] given by,

D = (α T |q| + φD m ) I + (α L -α T ) q ⊗ q |q| , (5)
where ⊗ denotes the dyadic product,

|•| the Euclidean norm, D m the molecular diffusion coefficient [L 2 T -1 ],
I the identity tensor, and α L and α T the longitudinal and transverse dispersivity, respectively [L]. The domain boundary is now splitted into the inflow part, ∂Ω in := {x ∈ ∂Ω : q • n < 0}, and the outflow part, ∂Ω out := ∂Ω \ ∂Ω in , such that ∂Ω in ∪ ∂Ω out = ∂Ω and ∂Ω in ∩ ∂Ω out = ∅. Thus, the boundary and initial conditions for the solute transport process are given by:

(ρcq -ρD∇c) • n = ρ in c in q • n, on ∂Ω in × (0, T ], (6a)
-D∇c • n = 0, on ∂Ω out × (0, T ], (6b)
c (•, 0) = c 0 , in Ω × {0}, (6c)
where ρ in and c in are the specified inflow density and mass fraction at the inflow boundary ∂Ω in , respectively.

State equation for the bulk fluid density

The state equations represent fundamental thermodynamic relationships reflecting some properties and physical characteristics of a fluid using state variables. In practice, the mathematical dependencies of the bulk fluid density and viscosity quantities are usually derived by using a total derivative representation with respect to the set of state variables. Here, we assume that the bulk fluid density and viscosity depend only on the mass fraction variable, i.e., ρ = ρ(c) and µ = µ(c), on the whole domain. We consider the simplest linearized model for the density :

ρ = ρ 0 (1 + β 0 c) , in Ω × (0, T ], (7)
where

β 0 = ρ -1 0 ∂ c ρ

Mathematical model

In this section, we propose to rewrite the set of partial differential equations governing variable-density flow with the corresponding initial and boundary conditions in a compact form. To this aim, we introduce suitable variables for the description of each of both processes. Thus, we use the following notations:

(Flow) s * = ρs, β * 0 = ρ 0 φβ 0 , q * = ρq and K * = ρK, (8a)
(Transport) φ * = ρφ, D * = ρD and c * in = ρ in c in ρ . ( 8b
)
Respecting the previous notations, the flow problem is now described by the simplified set of equations:

s * ∂ t ϕ + β * 0 ∂ t c + ∇• q * = 0, in Ω × (0, T ], (9a)
q * = -K * (∇ϕ + β 0 c∇z) , in Ω × (0, T ], (9b)
where the initial and Dirichlet boundary conditions (3a)-(3c) remain unchanged, and the Neumann condition (3b) becomes:

q * • n = q N , on ∂Ω N × (0, T ]. (10a)
The time-dependent advection-dispersion equation now can be written as follows:

γ∂ t (c) + ∇• (q * c -D * ∇c) = 0, in Ω × (0, T ], (11)
where γ(c) = φ * + β * 0 c, with the following boundary conditions

(cq * -D * ∇c) • n = c * in q * • n, on ∂Ω in × (0, T ], (12a)
-D * ∇c • n = 0, on ∂Ω out × (0, T ]. (12b)
Let us underline that the initial condition (6c) remains unchanged as well as the definitions of inflow ∂Ω in and outflow ∂Ω out boundary parts since ρ > 0.

Interior penalty discontinuous Galerkin methods

The class of Interior Penalty discontinuous Galerkin (IPDG) methods is now derived to solve the coupled non-linear system of partial differential equations. Precisely, we consider the Incomplete Interior Penalty (IIP), the Symmetric Interior Penalty (SIP) and the Non-symmetric Interior Penalty (NIP) variants for each physical process leading to nine possible combinations schemes. The main difference between these three IP variants lies in the symmetrization term used in the discrete bilinear form [START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF]. Before establishing the weak formulations, let us introduce some conventional notations that will be used throughout the remainder of the paper.

Some notations

We denote by 0 = t 0 < t 1 < t 2 < . . . < t N = T a non-uniform partition of the simulation time

interval (0, T ]. Let E h = {E i } N h be
h E = diam (E).
Let N H loc and N c loc be the number of degrees of freedom (DOFs) per element for head and mass fraction, respectively. Let F h be the set of all edges comprising interior and boundary edges such

that F h = F I h ∪ F B h .
Let Γ be an interior edge Γ ∈ F I h with size |Γ|, then there exist two adjacent elements

E ± Γ such that Γ = ∂E - Γ ∩ ∂E + Γ .
Further, we denote n ± Γ the unit normal vector on Γ pointing exterior to E ± Γ , respectively (see e.g., Fig. 1), and n Γ is coincident with n + Γ . Let Γ be a boundary edge Γ ∈ F B h , then there exist an element E ∈ E h such that Γ = ∂E ∩ ∂Ω. We define the set of Dirichlet, Neumann, inflow and outflow edges such that

F B h = F D h ∪ F N h = F in h ∪ F out h .
For clarity purposes, we denote by

F 0 h := F I h ∪ F D h
for the flow problem, and by

F + h := F I h ∪ F out h
for the solute transport problem. The DG method is based on the use of discontinuous approximations of discrete variables. To this aim, let us introduce the broken polynomial space V p h which is a finite dimensional space of discontinuous piecewise polynomial functions belonging to the broken Sobolev space H s (E h ) with s ≥ 1,

V p h := v ∈ L 2 (Ω) : v| E ∈ P p (E) ∀E ∈ E h ,
where P p (E) corresponds to the space of all polynomials of total degree ≤ p on E. Let ψ ± Γ be the traces of function ψ on the edge Γ between E ± Γ . Therefore, p ϕ and p c are taken to be the total degrees of discontinuous polynomials for approximating reference head and mass fraction, respectively. The jump and average of ψ on Γ can be now defined on the skeleton F h as follows,

E + Γ E - Γ Γ n Γ ∂Ω Ω = ∂Ω ∪ Ω Ei F I h F B h
ψ :=      ψ - Γ n - Γ + ψ + Γ n + Γ , if Γ ∈ F I h , ψn, if Γ ∈ F B h , (13a)
{ {ψ} } :=      ψ - Γ + ψ + Γ /2, if Γ ∈ F I h , ψ, if Γ ∈ F B h . (13b)
Let us precise that previous definitions of trace operators can be easily extended to any vector-valued functions ψ in the same way. For all x := (x 1 , . . . , x p ) ∈ R p with p ≥ 1, we define the uniform norm of x by x ∞ := max 1≤i≤p |x i | .

Spatial discretization of the flow process

The IPDG weak formulations for flow and transport are derived by multiplying the corresponding governing equation by appropriate test functions, integrating by parts over an element E ∈ E h (via the Green theorem), summing-up over all elements of the partition, and finally by imposing initial and boundary conditions. By performing the afore-mentioned operations on the flow problem ( 9), the weak form (in continuous time) is achieved and consists to seek, for any compatible c h (•, t) ∈ V pc h , the discrete variable

ϕ h (•, t) ∈ V pϕ h such that, (s * ∂ t ϕ h , v h ) E h + (β * 0 ∂ t c h , v h ) E h + a ( f ) h (ϕ h , v h ; c h ) = l ( f ) f (v h ; c h ) ∀v h ∈ V pϕ h , (14)
where f ∈ {0, ±1} is a symmetrization parameter, and the discrete operators a

( f ) h and l ( f ) f
are given by:

a ( f ) h (ϕ h , v h ; c h ) = (K * (c h )∇ϕ h , ∇v h ) E h -{ {K * (c h )∇ϕ h } }, v h F 0 h -f { {K * (c h )∇v h } }, ϕ h F 0 h + σ f Γ ϕ h , v h F 0 h , (15) l
( f ) f (v h ; c h ) = -(K * (c h )β 0 c h ∇z, ∇v h ) E h + { {K * (c h )β 0 c h ∇z} }, v h F 0 h -f K * (c h )∇v h • n Γ , ϕ D F D h + σ f Γ ϕ D , v h F D h -q N , v h F N h . (16)
The second, third and fourth terms on the right-hand side of ( 15) are respectively called consistency, symmetry and penalty terms. The role of the parameter f consists of regulating the symmetry term impact inside the bilinear form. Here, σ f Γ corresponds to the penalty parameter on Γ, and it is defined on the mesh skeleton as follows:

σ f Γ =        2τ - Γ τ + Γ τ - Γ + τ + Γ if Γ ∈ F I h τ E,Γ if Γ ∈ F D h ( 17
)
where τ ± Γ = τ E ± ,Γ denotes the transmissibility coefficient on Γ = ∂E -∩ ∂E + . For all E ∈ E h , and for any Γ ∈ ∂E, this parameter is given by,

τ E,Γ = σ f 0 |Γ| |E| (p ϕ + 1)(p ϕ + 2)κ E,Γ , (18)
where κ E,Γ = n Γ K * E n Γ represents the normal diffusivity on Γ and σ f 0 ≥ 0 is a user-dependent parameter.

The unified DG formalism of the variable density flow problem ( 14) includes the Incomplete Interior Penalty (IIP), the Symmetric Interior Penalty (SIP), the Non-symmetric Interior Penalty (NIP), and the Oden-Babuška-Baumann (OBB) methods. All these variants are deduced by choosing precisely both parameters f and σ f 0 : IIP ( f = 0 and σ f 0 > 0), SIP ( f = 1 and σ f 0 > 0), NIP ( f = -1 and σ f 0 > 0), and OBB ( f = -1 and σ f 0 = 0). Regardless of the choice of IP methods, it is a fundamental aspect at this stage to reconstruct the Darcy's velocity field q * . To this aim, we shall use the favourable local conservation properties associated with IPDG methods. Let us denotes by q * h the discrete approximation of the velocity field q * on the partition E h . Two description levels are used to approximate q * (i) at the element-level E of the mesh E h , and (ii) at the edge-level Γ of the mesh skeleton F h .

q * h =      q * h , for all E ∈ E h , q * h for all Γ ∈ F h , (19)
where q * h = -K * (c h ) (∇ϕ h + β 0 c h ∇z) and the numerical flux q * h is defined as follow:

q * h =            { {q * h } } + σ f Γ ϕ h , if Γ ∈ F I h q * h + σ f Γ (ϕ h -ϕ D ) n Γ , if Γ ∈ F D h q N n Γ , if Γ ∈ F N h (20)

Spatial discretization of the transport process

This section is devoted to the presentation of the IPDG discretization of the transport equation respecting the imposed boundary conditions. We use a similar approach to the previous section to handle the dispersive part that we combine with a traditional upwind-strategy to control the advective part. Thus, the spatial discretization of the transport problem (in continuous time) consists to seek, for any reconstructed

q * h (•, t), the discrete variable c h (•, t) ∈ V pc h such that, (γ∂ t c h , w h ) E h + b ( t) h (c h , w h ; q * h ) = l ( t) t (w h ; q * h ) , ∀w h ∈ V pc h , (21)
where t ∈ {0, ±1}, the discrete bilinear form b

( t) h is given by, b ( t) h (c h , w h ; q * h ) = (D * (q * h )∇c h -c h q * h , ∇w h ) E h + q * h { {c h } }, w h F + h -{ {D * (q * h )∇c h } }, w h F I h -t { {D * (q * h )∇w h } }, c h F I h + σ t Γ c h , w h F I h , (22a)
and the linear form l

( t) t by, l Γ τ - Γ + τ + Γ + 1 2 |q * h • n Γ | , ∀Γ ∈ F I h , (23)
where the parameters τ ± Γ = τ E ± ,Γ are now derived using the hydrodynamic dispersion tensor D * and the polynomial order p c of the discrete variable c h . For all E ∈ E h , and for any Γ ∈ ∂E, it is given by,

τ E,Γ = σ t 0 |Γ| |E| (p c + 1)(p c + 2)κ E,Γ , (24)
where κ E,Γ = n Γ D * E n Γ represents the normal dispersion on Γ and σ t 0 ≥ 0 is a user-dependent parameter.

Despite the natural skew-symmetry of the advection-dispersion operator, we have included the additional symmetry term of the dispersive operator which is controlled here by the parameter t . The resulting DG method is also called IIP, SIP, and NIP discretizations of the solute transport problem if t is equal to 0, 1, and -1, respectively [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF][START_REF] Klieber | Adaptive simulations of two-phase flow by discontinuous Galerkin methods[END_REF][START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF].

Time-discretization of coupled flow-transport problem

The IPDG approach for solving the coupled flow-transport process yields to a set of differential and algebraic equations. To this aim, let us consider the discrete variables

ϕ h (•, t) ∈ V pϕ h and c h (•, t) ∈ V pc h that
we decompose as follows:

ϕ h = N f i=1 ϕ i (t)v i and c h = Nc i=1 c i (t)w i ( 25
)
where

N f = dim(V pϕ h ) and N c = dim(V pc h ).
Here, (v i ) i=1,...,N f and (w i ) i=1,...,Nc denote the set of trial functions of approximation spaces V pϕ h and V pc h , respectively. The matrix form of the global discrete problem can be recast as follows:

  S C 0 T     φh ċh   +   A ( f ) 0 0 B ( t)     ϕ h c h   =   L f L t   ( 26
)
where the dot symbol denotes the partial time-derivative operator, and

ϕ t h = [ϕ 1 , . . . , ϕ N f ] and c t h = [c 1 , . . . , c Nc ],
the vector of degrees of freedom of the discrete variable ϕ h and c h , respectively. The blockmatrices associated with [START_REF] Ern | Discontinuous Galerkin approximation of two-phase flows in heterogeneous porous media with discontinuous capillary pressures[END_REF] are given by,

S |ij = (s * (c h )v j , v i ) E h , (27a)
C |ik = (β * 0 (c h )w k , v i ) E h , (27b)
T |kl = (γ(c h )w l , w k ) E h , (27c)
A ( f ) |ij = a ( f ) h (v j , v i ; c h ) , (27d)
B ( t) |kl = b ( t) h (w l , w k ; q * h ) , (27e)
L f|i = l ( f ) f (v i ; c h ) , (27f)
L t|k = l ( t) t (w k ; q * h ) . ( 27g
)
where indexes i, j = 1, . . . , N f , k, l = 1, . . . , N c , and block-matrices

A, S ∈ R N f × R N f , B, T ∈ R Nc × R Nc ,
and

C ∈ R N f × R Nc .
Due to the discontinuous nature of V pϕ h and V pc h , all computations (matrix assembly) can be done locally at the element level. We underline that the IPDG method leads to a very compact discretization stencil consist of a given mesh element and its direct neighbors. A popular way for solving the coupled differential-algebraic system (26) consists in splitting the numerical treatment of both physical sub-processes (flow and transport) and to treat them sequentially via an adapted flux reconstruction.

Usually, for each time step, an iterative procedure is applied between the flow and transport operators until convergence. This approach could be more efficient in CPU time than the fully coupled approach that proceeds by solving both flow and transport operators simultaneously. However, the slow convergence of the iterative procedure between both operators affects significantly the performance of the solution.

To avoid this problem, we use the non-iterative approach suggested by Younes and Ackerer [START_REF] Younes | Empirical versus time stepping with embedded error control for densitydriven flow in porous media[END_REF]. In this paper, the time discretization is performed by a first-order Backward Differentiation Formula (BDF 1 ) and a direct solver using an unsymmetric-pattern multifrontal method and a direct sparse LU factorization (UMFPACK) is implemented for solving linear systems. Accuracy of the time discretization is insured using a proper time step management where the time step length is controlled by the temporal truncation error. Note that high order BDF can also be used and allow large time steps which may further improve the efficiency of the simulations [START_REF] Younes | Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods[END_REF]. The main steps of the non-iterative procedure used to solve the coupled flow and transport operators are described here. We denote by c n h and ϕ n h the discrete approximation of c(•, t n ) and ϕ(•, t n ) at the instant t n with 0 ≤ n ≤ N , respectively. Similarly, we denote by ϕ n h and c n h the corresponding vector of degrees of freedom. We describe below the staggered procedure for switching from t n to t n+1 and calculating the (n + 1)-th time step ∆t n+1 .

• Step 1. Solve the transport problem: Given c n h and q * h n , seek c n+1 h solving the equation,

T(c n h ) ∆t n + B ( t) (q * h n ) c n+1 h = L t (c n h ) + T(c n h ) ∆t n c n h , (28)
where ∆t n = t n+1t n .

• Step 2. Solve the flow problem: Given initial values c n+1 h , seek ϕ n+1 h solving the equation,

S(c n+1 h ) ∆t n + A ( f ) (c n+1 h ) ϕ n+1 h = L f (c n+1 h ) + S(c n+1 h ) ∆t n ϕ n h - C(c n+1 h ) ∆t n (c n+1 h -c n h ). (29)
• Step 3. Flux reconstruction: Given c n+1 h and ϕ n+1 h , we reconstruct the global velocity field q *

h n+1 = [q * h n+1 , q * h n+1
] on the whole domain as follows :

-Mesh element: For all A ∈ E h , we compute interior approximations as follows,

q * h n+1 = -K * (c n+1 h ) ∇ϕ n+1 h + β 0 c n+1 h ∇z . (30)
-Mesh skeleton: For all Γ ∈ F 0 h , we compute boundary approximations as follows,

q * h n+1 =      { {q * h n+1 } } + σ f Γ ϕ n+1 h , if Γ ∈ F I h , q * h n+1 + σ f Γ ϕ n+1 h -ϕ n+1 D n Γ , if Γ ∈ F D h . (31)
• Step 4. Update the time step: The truncation error e n+1 measures the difference between first-and second-order temporal approximations of the concentration [START_REF] Younes | Empirical versus time stepping with embedded error control for densitydriven flow in porous media[END_REF]:

e n+1 := 1 2 c n+1 h -c n h + ∆t n ∆t n-1 (c n h -c n-1 h ) . (32)
The uniform norm e n+1 ∞ is used to accept or reject the time step with respect to the user tolerance γ tol (set at 10 -3 in this work).

-If e n+1 ∞ ≤ γ tol , the time step is accepted and the next one is estimated using,

∆t n+1 = ∆t n min κ γ tol e n+1 ∞ , r max , (33)
where κ(= 0.95) is a safety factor and r max (= 2) is the maximum allowed time step.

-Else, the time step is repeated with a smaller step size using the latest error estimate:

∆t n j+1 = ∆t n j max κ γ tol e n+1 j ∞ , r min , (34)
where j indexes the consecutive time step estimates and r min (= 0.1) is the minimum allowed time step.

• Step 5. Update the initial guess:

n ← n + 1, c n h ← c n+1 h , ϕ n h ← ϕ n+1 h
, ∆t n ← ∆t n+1 and back to

Step 1.

Compatibility analysis & Frolkovič-Knabner procedure

Traditional algorithms employ operator-splitting to treat the coupled flow and transport processes sequentially and separately. However, when solving coupled processes, the discretization method employed to solve the flow is crucial since the approximate velocity field can strongly influence the positive features of the transport method. In 2003, Dawson et al. [START_REF] Dawson | Compatible algorithms for coupled flow and transport[END_REF] analyzed compatible DG-algorithms for the coupled flow and transport problem. Specifically, they established some minimal requirements on the flow DG schemes to maintain optimal accuracy and conservation properties of the DG algorithms used for transport. Firstly, the authors proved that DG methods employed for the transport process must imperatively respect the zero-order accuracy requirement. This latter condition measures the ability of the DG method to reproduce a constant field by replacing the true velocity field q * by its discrete approximation q * h in the transport process. Particularly, Sun and Wheeler [START_REF] Sun | Discontinuous Galerkin methods for coupled flow and reactive transport problems[END_REF] analyzed this criterion in the context of primal IPDG methods, and they established that all schemes namely, IIP-, NIP-and SIP-DG variants, verify it. However, it also imposes that the discrete velocity field verifies a compatibility requirement, implying that the DG method employed for the resolution of flow can not be chosen arbitrarily i.e.,

(q * h , ∇w h ) E h + q * h , w h F h = q * • n Γ , w h F in h ∀w h ∈ V pc h . (35)
Let us note that the compatibility condition ( 35) corresponds here to a stronger form of the local conservative principle since it holds for any polynomial functions in the space V pc h , and not only for piecewise constants.

Thus, in virtue of [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF], the IIP-DG scheme method for flow is the only one compatible with primal IPDG methods for transport in the sense defined in [START_REF] Dawson | Compatible algorithms for coupled flow and transport[END_REF]. In compliance with Sun and Wheeler [START_REF] Sun | Discontinuous Galerkin methods for coupled flow and reactive transport problems[END_REF], the total number of permissible combinations to solve VDF equations is now reduced from nine to three compatible ones: IIP-IIP, IIP-SIP, and IIP-NIP. Any other combination of (non-compatible) methods can cause severe struggles such as a loss of accuracy or conservative properties, and may not even converge. For solving the VDF equations, we can choose related polynomial degrees of approximation for flow and transport i.e., p ϕ ∼ p c , however, optimal or nearly optimal solution convergence can be attained when p ϕ = p c , in coupled flow and transport processes [START_REF] Sun | Discontinuous Galerkin methods for coupled flow and reactive transport problems[END_REF].

Besides, in buoyancy-driven problems, existence of additional gravity term in Darcy's velocity make a constraint on use of equal degrees of approximation for flow and transport due to the confliction between the head gradient and gravity terms of Darcy in parts of the domain where the velocity is zero or nearly zero

[see e.g. [START_REF] Diersch | Variable-density flow and transport in porous media: approaches and challenges[END_REF][START_REF] Frolkovič | Consistent velocity approximation for density driven flow and transport[END_REF][START_REF] Voss | Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater-saltwater transition zone[END_REF]. To overcome this problem, a consistent velocity approximation based on Frolkovič-Knabner (FK) approximation [START_REF] Frolkovič | Consistent velocity approximation for density driven flow and transport[END_REF] was developed for the IPDGs presented here for different polynomial degrees. The FK allows us to approximate flow and transport with the same degree of polynomials and also reduces the computational effort, considerably.

Results and discussion

A new numerical scheme is developed to solve the VDF equations with the DG method. This scheme is implemented in a numerical code. The first goal of this section is to validate the new developed scheme and to verify the correctness of the developed code. To do so, we use common benchmarks dealing with different configurations and applications of VDF and we compare the results of the new code against experimental results and semi-analytical solutions. The first benchmark is the Goswami-Clement laboratory experimental problem [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF] that deals with saltwater intrusion in a rectangular experimental tank. This problem is used to test the new developed code in the case of specified-head boundary conditions (Dirichlet type). The comparison against experimental data would give an overall assessment of the robustness of the mathematical model in reproducing physical processes but cannot confirm the accuracy of the developed numerical solution and the correctness of the developed numerical code. The later can be well assessed by comparing numerical results against exact-analytical solutions. As there is no analytical solution for the VDF model, we compare the results of the new developed code against the semi-analytical solution of the Henry problem [START_REF] Henry | Effects of dispersion on salt encroachment in coastal aquifers: sea water in coastal aquifers[END_REF], which is a common benchmark for seawater intrusion. Several semi-analytical solutions of the Henry problem have been suggested in the literature [START_REF] Henry | Effects of dispersion on salt encroachment in coastal aquifers: sea water in coastal aquifers[END_REF][START_REF] Segol | No Title. Classic Groundwater Simulations: Proving and Improving Numerical Models[END_REF][START_REF] Simpson | Theoretical analysis of the worthiness of Henry and Elder problems as benchmarks of density-dependent groundwater flow models[END_REF][START_REF] Younes | A Semi-Analytical Solution for the Reactive Henry Saltwater Intrusion Problem[END_REF][START_REF] Zidane | The Henry semianalytical solution for saltwater intrusion with reduced dispersion[END_REF], here we use the recent solution developed by Fahs et al. [START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF] as, in contrast to the previous solutions, it involves velocitydependent dispersion. The Henry problem allows not only comparison against semi-analytical solution but also checking the model in a case dealing with specified-flux boundary conditions (Neumann type).

Both Goswami-Clement and Henry problems deal with mixed convection as the flow is generated by a head gradient or forced conditions. In order to investigate the new developed code in a case involving purely natural convection, we consider the problem of natural convection in porous square cavity with vertically heated walls. In this case, VDF is thermally driven. The flow is generated by the temperature gradient. This problem is commonly used to represent natural convection in porous media in either cases of vertical (horizontally heated walls) or horizontal temperature gradient (vertically heated walls). We do not consider horizontally heated wall as the problem could be unstable. Unstable problems are not suitable for benchmarking as they could have multiple solutions [START_REF] Prasad | Using quantitative indicators to evaluate results from variabledensity groundwater flow models[END_REF]. Our new developed code is compared against the stable semi-analytical solution developed by [START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF]. For validation, we consider cases dealing with relatively small Rayleigh number. Beside verification, the problem of porous square cavity would also highlight the flexibility of the new developed model in treating either solute or thermal problems and even coupled thermohaline problems.

The Results section aims also at highlighting the performance and accuracy of the developed DG scheme in solving VDF equations. To do so, we consider a case of the porous cavity problem with high Rayleigh number and we compare the new developed model against a commercial model based on standard finite element method (COMSOL Multiphysics). Fahs et al. [START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF] have shown that such a case is computationally challenging as accurate solution of this case is beyond the capacity of current models. Standard numerical techniques used in current models could lead to spurious oscillations or numerical diffusion.

Validation: The Goswami-Clement experimental problem

Herein, to validate our model we simulate the laboratory experiments originally performed by Goswami and Clement [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF]. These experiments have been developed with the purpose of evaluating the densitydependent models. The experimental setup composed of a rectangular flow tank comprising three distinct chambers. The central chamber contains an unconfined, homogeneous, porous medium and two constanthead chambers on the left and right sides containing saltwater and freshwater, respectively. Goswami and Clement [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF] completed their experiments by recording the data in both transient and steady state conditions. As a result, in addition to the steady state validation, the transient test of numerical models is also possible. The saltwater intrusion experiments were conducted under three transient phases ended

with steady state conditions. The first phase was established by setting the right freshwater head equal to 26.7 cm and fixing the left equivalent freshwater head proportional to the saltwater head of 25.5 cm. Under these conditions, the first steady state condition (SS 1 ) was obtained for the salt wedge (initial phase). Then the freshwater head was instantaneously reduced to 26.2 cm and maintained until the second steady state condition (SS 2 ) is reached (advancing phase). Finally, the freshwater head was increased to 26.55 cm which forces the salt wedge to recede and the conditions was kept until the third steady state condition (SS 3 ) was established (receding phase).

In the present study, we tested our DG model by comparing the simulation results against experiments data in both transient and steady state conditions. For numerical simulations, the unconfined porous medium was assumed to be confined with the boundary conditions depicted in Fig. 2, [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF]. The domain was discretized using 44K structured right-angled triangles with equal grid spacing and linear polynomial approximations were used for both variables i.e., head and mass fraction. The initial mass fraction and freshwater head were set to zero and 26.7 cm, respectively. Table 1 illustrates the physical parameters used for simulation of the experiments.

Confined aquifer

Saltwater level: 25.5 cm

Groundwater discharge:

D * ∇c • n = 0
Saltwater inflow to SS 3 . In addition to the experimental data from measurements of salt wedge location, the numerical results of SEAWAT [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF] were superimposed for comparison. As can be seen, our results are in excellent agreement with the experimental data. The maximum difference between the results and experimental data respects to SS 2 . This is while for the SEAWAT model, the difference exists for all the three steady state conditions. In general, our model can provide a closer prediction of salt wedge, particularly in the toe, for SS 1 to SS 3 . Note that the results from IIP-IIP, IIP-SIP and IIP-NIP were closely similar [START_REF] Isa-Abadi | Modeling of seawater intrusion into heterogeneous coastal aquifers using discontinuous Galerkin method[END_REF].

(cq * -D * ∇c) • n = c * in q * • n
Flood: DG solution (present study) Dotted: Experimental data [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF] Dashed: SEAWAT solution [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF] Fig. 3. 10%, 50% and 90% isochlors of DG solution in comparison with Goswami experimental data [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF] and 50% SEAWAT isochlor [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF] in the (a) first steady state (SS 1 ), (b) second steady state (SS 2 ), and (c) third steady state (SS 3 ) conditions

Here, for the first transient test, the advancing phase of salt wedge was also considered. Again, the most difference is observed when the salt wedge is closer to SS 2 position (Fig. 5). This condition is also true for SEAWAT [START_REF] Isa-Abadi | Modeling of seawater intrusion into heterogeneous coastal aquifers using discontinuous Galerkin method[END_REF]. Finally, a flux test was also presented to compare the freshwater fluxes predicted by the model against the steady state flux measurements. Table 2 reports the measured fluxes and the results from our DG code and SEAWAT. For the DG, the most difference is happened in SS 3 , though, this is not significant and is less than 2 percent. Generally, the fluxes estimated by our model compared to SEAWAT have closer agreement with the experimental measurements.

Verification against semi-analytical solution: The Henry problem

Henry [START_REF] Henry | Effects of dispersion on salt encroachment in coastal aquifers: sea water in coastal aquifers[END_REF] was the first to model seawater intrusion phenomenon in a variable density system. His simplified rectangular model is constructed based on a confined, homogeneous, isotropic 2D aquifer in a vertical section (Fig. 6). As shown, the left landward boundary is recharged by a constant rate of freshwater (Q) and the right boundary is ended to the sea such that the top aquifer boundary is just leveled at the sea surface. Saltwater intrudes from the seaside until an equilibrium is established with the landward inflowing freshwater. Henry developed a semi-analytical solution based on the double Fourier series for the problem under the steady state condition. Due to the availability of the analytical solution, for many years, various attempts have been made [see the review in 76, 95] aimed at utilizing the Henry problem as benchmark for variable density flow models. The first semi-analytical solution, developed by Henry [START_REF] Henry | Effects of dispersion on salt encroachment in coastal aquifers: sea water in coastal aquifers[END_REF], has been limited to high diffusion for which the buoyancy processes are dominated by diffusion. Voss et al. [START_REF] Voss | Three-dimensional benchmark for variable-density flow and transport simulation: matching semi-analytic stability modes for steady unstable convection in an inclined porous box[END_REF] showed that this solution is a necessary but not a sufficient test for a simulator to present variable-density physics. Several further studies suggested new semi-analytical solutions that are more sensitive to the variable-density flow physics by reducing the imposed inland flow or the diffusion coefficient [START_REF] Simpson | Improving the worthiness of the Henry problem as a benchmark for density-dependent groundwater flow models[END_REF][START_REF] Younes | A semi-analytical solution for saltwater intrusion with a very narrow transition zone[END_REF][START_REF] Zidane | The Henry semianalytical solution for saltwater intrusion with reduced dispersion[END_REF].

Here we use the semi-analytical solution developed by Fahs et al. [START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF] that includes velocity-depended dispersion. This solution is more realistic than previous solutions in representing seawater intrusion processes. It is more suitable for code verification as it includes velocity dependent dispersion. Here for verification purposes, we consider the three cases presented in Fahs et al. [START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF]. The first one is similar to the standard solution suggested by Henry [START_REF] Henry | Effects of dispersion on salt encroachment in coastal aquifers: sea water in coastal aquifers[END_REF]. It deals with pure molecular diffusion. The second one in-cludes velocity-dependent dispersion, but with exaggerated dispersion coefficients leading to a wide mixing zone. These two cases are relatively simple from computational point of view as they do not involve sharp solutions. They are useful in verifying the correctness of the developed code. The third case deals with small dispersion coefficients, leading to a very narrow mixing zone. As shown in Fahs et al. [START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF], this case is computationally challenging because of its high sensitivity to the numerical scheme. Fahs et al. [START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF] showed that, for this case, standard finite element method leads to inaccurate results due to numerical diffusion.

The physical parameters used in the three test cases are given in Table 3. Note that the three cases can be derived from the standard Henry problem by changing one or more of the following parameters: domain length (l), molecular diffusion (D m ) and dispersion coefficients (α L , α T ).

For the purpose of numerical solution, following Fahs et al. [START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF], the domain was discretized using structured meshes comprising 7500, 12288 and 76800 triangles for test cases 1 to 3, respectively. For all the test cases except the third one which utilizes quadratic polynomial for head, the linear polynomial approximations were applied for both variables. The imposed boundary conditions are depicted in Fig. 6.

All combinations of DG schemes are used for the simulation of Henry problem, however, we observed no significant difference between the results. Main isochlors (10%, 50% and 90%) as well as the velocity field for the three test cases are represented on Fig. 7. It can be observed that there is an excellent agreement between generated DG isochlors and those of semi-analytical solution throughout the total aquifer thickness.

The third case highlights the capacity of the developed DG solution in reproducing a narrow mixing zone that could be very sensitive to numerical dispersion [START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF]. A more quantitative comparison is performed using scalar indicators representing the seawater intrusion metrics. The results are summarized in Table 4.

For test case 1, the difference between the DG and semi-analytical solutions is less than 1.5%. The highest difference is observed for the width of the mixing zone. The post treatment procedure used for the evaluation of this zone depends on the computational mesh. For the second and third test cases, the discrepancy between the solutions is less than 6%, except for the mixing zone and the total flux. The small discrepancies between the solutions regarding L s and L toe confirms that DG avoids the problem of numerical dispersion. In fact, numerical dispersion could overestimate the diffusion processes and leads to an overestimation of L toe and L s .

Verification against semi-analytical solution:

The problem of natural convection in a porous enclosure An important application of the VDF is in natural convection in porous enclosure. In this case, the density variation is related to a thermal gradient. This problem can be modeled by a system of equations equivalent to the solute VDF problems. For thermal problems, the porosity is analogue to the equivalent specific heat of the porous domain and the molecular diffusion is analogue to the equivalent thermal diffusivity of the porous domain. In thermal application, the expansion coefficient in [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF], should be negative as the density decreases with the increase of temperature. Thus the developed DG scheme can be applied to solve thermal density-driven problems. The previous test cases deal with mixed convection as the flow Ltoe is the dimensionless distance to the seaside boundary that seawater intrudes along the aquifer bottom (measured based on 50% isochlor), Ls is the dimensionless distance between 10 and 90% isochlors along the aquifer bottom, Wmz is the vertically dimensionless average width of the mixing zone, Z I is the dimensionless vertical coordinate of the point at which flow direction at seaside boundary alters adversely, Qs is the saltwater inflow flux from seaside boundary divided by freshwater recharge (Q) .

is either imposed or generated by an imposed pressure gradient. The main goal of this section is to test the new developed scheme in cases involving natural convection, where flow is induced by density gradient. Thus, we consider the well-known benchmark of natural convection in a porous cavity (Fig. 8). The domain of the problem is a unit square filled with homogeneous porous medium. All sides are impervious respective to fluid. Right and left sides of the domain have different temperatures with 1 and 0 dimensionless temperatures, respectively. Other boundary conditions and data are depicted in Fig. 8. The main factor of fluid density changes is the temperature differences due to heating the fluid from the boundaries. The convective flow and heat transfer processes are governed by the dimensionless Rayleigh number (Ra = Klβ 0 ∆c/φD m , where l is the domain length) expressing the ratio of the buoyancy forces to the viscousdiffusion effects.

Fahs et al. [START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF] developed a semi analytical solution for this problem which is used here to verify the proposed DG scheme. We consider cases dealing with low and average Rayleigh numbers (Ra = 100 and Ra = 1000). The physical parameters used in the simulations are given in Table 5.

In the previous test for Henry problem, DG was evaluated for linear approximation. Here, higher order Flood: DG (present study) Lines: Semi-analytical [START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF] Flood: DG (present study) Lines: Semi-analytical [START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF] Flood: DG (present study) Lines: Semi-analytical [START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF] (a) (b) (c) This section aims at examining the performance of the new developed DG scheme against conventional methods. Thus, we consider a computationally challenging case and we simulate this case with the DG scheme and with a standard finite element solution obtained using COMSOL. We simulate the case on natural convection in porous enclosure under high Rayleigh regime (Ra = 10000). Fahs et al. [START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF] showed that in this case, conventional numerical methods suffer from numerical dispersion and/or nonphysical oscillations. They showed that accurate simulation of this test case is beyond the capacity of conventional models. In their work, the author proposed a reference solution based on the Fourier series method. We simulate this case using the DG scheme and COMSOL. For the DG scheme we use a computational mesh of 10K elements and we test different polynomial orders from 1 to 4. In COMOS, we use two levels of mesh refinement consisting of about 100K and 600K elements. COMSOL simulations are performed with first polynomial order. We use the average Nusselt number to evaluate the performance of both models.

The average Nusselt number represents the dimensionless heat flux to the domain. It is calculated by:

Nu = ∂Ω|x=0 Nu = ∂Ω|x=0 ∇c • n [29].
Table 6 summarizes the Nusselt number for the DG scheme in comparison with those of semi-analytical

Flood: DG (present study)

Lines: Semi-analytical [START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF] Flood: DG (present study)

Lines: Semi-analytical [START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF] (a) (b) solution and COMSOL-Multiphysics. We should mention that, by default in COMSOL, second polynomial order is used for flow and first order is used for transport. As shown, the Nu achieved for DG approaches to that of analytical as the order of polynomial approximation (the number of DOFs) is increased until at order 4 including totally 300, 000 DOFs, Nu is extremely close to analytical value. Fig. 10 pinpoints how the relative error on Nusselt number decreases as the number of DOFs is increased. COMSOL-Multiphysics is a standard finite element software and its results were also included to demonstrate how DG is superior to standard finite element method. We observe for high order DG solution, Nu is more close to the analytical value, though, less than a quarter of the number of COMSOL-Multiphysics' DOFs was used.

The DG scheme intrinsically has less numerical dispersion and this is why a more accurate solution is achieved for, even with a much less number of DOFs. Fig. 11 depicts the isotherms for both models against the semi-analytical solution. As illustrated, DG isotherms (for order 4) have closer agreement with the semi-analytical solution than COMSOL.

The results of this test case show clearly a significant improvement in the accuracy of the solution due to the DG method. This is in agreement with several previous studies on flow and transport in porous media, confirming the high accuracy of the DG methods. However, in general the DG method is seen to be more consuming in CPU time and memory than standard methods, due to the high number of DOFs per element of the computational mesh. The results here show that for VDF problems, this is not totally true.

In fact, a reliable comparison should respect the balance between computational efficiency and solutions accuracy. In this logic, Table 6 confirms the superiority of the DG methods for the simulation of VDF problems.

Conclusion

Several recent studies have highlighted the advantages of the DG methods in solving partial differential equations governing flow or transport processes in porous media. This method has been never used to solve the full system of VDF equations. Yet, simulations of VDF requires high accurate numerical methods as the accuracy of the numerical solutions is highly sensitive to the numerical scheme. The first objective of this work is to show how the equations of VDF can be solved with the DG method. This requires extension of the DG method to solve nonlinear equations coupling flow and convective/dispersive transport equations under variable density. We develop, in a unified format, the general class of Interior Penalty DG (IPDG) methods to solve VDF equations. We test symmetric, non-symmetric and incomplete IPDG methods to discretize both head and concentration variables. Numerical experiments are performed to validate the developed DG scheme and to examine its performance in solving VDF problems.

The developed numerical scheme is implemented in a numerical code which is used to simulate common benchmarks dealing with VDF. The simulations are performed to validate and verify the developed DG scheme. The benchmarks are selected to cover different types of boundary conditions and convective flow processes (either, mixed, natural, solute or thermal). Validation is performed by comparing the results of the developed DG scheme with experimental data, based on the Goswami-Clement experimental problem.

Good agreement has been found. The new developed DG scheme reproduces the experimental data better than the finite element SEAWAT code. Verifications are performed by comparing the code against the semi- with semi-analytical solutions for cases involving wide mixing zone or small Rayleigh numbers. This confirms the correctness of the developed code and DG scheme. Furthermore, close agreement has been obtained for Henry problem in the case of narrow mixing, which illustrates the promising results of the DG method in solving VDF with sharp fronts. These numerical experiments show that, whole algorithm combinations, i.e. IIP-IIP, IIP-SIP, and IIP-NIP have no significant differences in the solutions.

To highlight the performance of the DG method in solving the equations of VDF, we simulate the problem of natural convection at a high Rayleigh number (10,000). Accurate solution of this problem is out of the capacity of current models based on conventional methods. We compared the DG results to a semi-analytical solution and to a finite element solution obtained using COMOSL multi-physics. The results prove the promising performance of the DG method at high polynomial order as it provides higher accuracy than standard finite element method with reduced number of degree of freedom. The DG method outperforms the conventional method that we tested.

This study is the first step towards application of DG methods in the area of VDFs. Future work should target applying other nice properties of DG such as adaptivity in mesh and polynomial approximation

  [L -1 ], c the solute mass fraction [-], φ the kinematic porosity [-], and ∂ x the partial derivative operator with respect to the xvariable. The boundary ∂Ω is divided into a Dirichlet part, ∂Ω D , and a Neumann part, ∂Ω N , such that ∂Ω D ∪ ∂Ω N = ∂Ω and ∂Ω D ∩ ∂Ω N = ∅. Thus, the boundary and initial conditions for the flow process are

  corresponds to the expansivity coefficient resulting from a change of the solute mass fraction [-]. Different mathematical models of state equations have been proposed in the literature, including the thermal effect and based on a polynomial expansion of the state variables. The viscosity is considered as constant in the numerical simulations developed in this work, but the developed DG model can handle variable viscosity.

Fig. 1 .

 1 Fig. 1. Schematic discretization of Ω and representation of elements, interior and boundary edges

Freshwater 26 .Fig. 2 .

 262 Fig. 2. Geometry and boundary conditions for the Goswami-Clement experimental setup

Fig. 4

 4 compares the results of numerical simulations with experimental measurements at times 5, 15, 55 min after starting the second phase. Similar to the steady state test, the results match well with the experimental data at the different times. The results show that the more salt wedge advances over time, the less agreement appears due to approaching to SS 2 position. The results of SEAWAT model were also presented in the figure for comparison. Similarly, for SEAWAT, the differences was increased in the later times. As a second part of transient test, we performed a comparison of numerical results with experimental data in receding phase.

Flood:

  DG solution (present study) Dotted: Experimental data[START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF] Dashed: SEAWAT solution[START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF]

Fig. 4 .

 4 Fig. 4. 10%, 50% and 90% isochlors of DG solution in comparison with Goswami experimental data [35] and 50% SEAWAT isochlor [35] at (a) 5 min, (b)15 min, and (c) 55 min of advancing phase (moving from SS 1 toward SS 2 )

Fig. 5 .

 5 Fig. 5. 10%, 50% and 90% isochlors of DG solution in comparison with Goswami experimental data [35] and 50% SEAWAT isochlor [35] at (a) 10 min, (b)15 min, and (c) 25 min of advancing phase (moving from SS 2 toward SS 3 )

Fig. 6 .

 6 Fig.6. Simplification of seawater intrusion by Henry[START_REF] Henry | Effects of dispersion on salt encroachment in coastal aquifers: sea water in coastal aquifers[END_REF] and the implemented boundary conditions

Fig. 7 .Fig. 8 . 4 . 4 .

 7844 Fig. 7. Main iscohlors and velocity field for the test case 1 (a), test case 2 (b), and test case 3 (c) of the Henry problem: comparison between the DG and semi-analytical solutions

Fig. 9 .

 9 Fig. 9. Comparison of DG solution with that of semi-analytical for (a) Ra = 100 and Ra = 1000

Fig. 10 .

 10 Fig. 10. Relative error on Nu versus number of DOFs

Flood:Fig. 11 .

 11 Fig. 11. Isotherms concluded from DG compared to that of analytical and COMSOL-Multiphysics

Table 1 .

 1 Physical parameters for the Goswami-Clement experimental problem

	Parameter	Value	Unit
	K	0.0122	ms -1
	Dm α L	1.0 × 10 -9 0.001	m 2 s -1 m
	α T s φ ρ0	0.0001 1.0 × 10 -5 0.385 1000	m m -1 -kg m -3
	ρs	1026	kg m -3
	µ0	0.001	kg m -1 s -1
	βµ β0	0 0.026	--

Table 2 .

 2 Numerical fluxes (cm 3 • s -1 ) in comparison with the experimental data and SEAWAT for the three steady state

	conditions of Goswami-Clement problem
		SS1	SS2	SS3
	Experimental measurement [35]	1.42	0.59	1.19
	DG (present study)	1.41	0.6	1.17
	SEAWAT [35]	1.46	0.59	1.13

Flood: DG solution (present study) Dotted: Experimental data

[START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF]

Dashed: SEAWAT solution

[START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF]

Table 3 .

 3 Physical parameters used for the simulations of the three test cases of the Henry problem

	Parameter		Value		Unit
		Case 1	Case 2	Case 3
	K Dm α L	0.01 × I 18.86 × 10 -6 0.0	0.01 × I 9.43 × 10 -8 0.1	0.01 × I 9.43 × 10 -8 0.001	ms -1 m 2 s -1 m
	α T s	0.0 0	0.01 0	0.0001 0	m m -1
	φ ρ0	0.35 1000	0.35 1000	0.35 1000	-kg m -3
	ρs	1025	1025	1025	kg m -3
	Q µ	6.6 × 10 -5 0.001	6.6 × 10 -5 0.001	6.6 × 10 -5 0.001	ms -1 kg m -1 s -1
	β0	0.025	0.025	0.025	-

Table 4 .

 4 Computed seawater intrusion metrics for Henry test cases[START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF]

	Metrics *		Test case 1		Test case 2		Test case 3
		DG	Semi-analytical	DG	Semi-analytical	DG	Semi-analytical
	L T oe	0.624	0.624	1.246	1.256	1.604	1.594
	Ls	0.752	0.751	0.393	0.368	0.087	0.088
	Wmz	0.756	0.757	0.323	0.295	0.088	0.09
	Z1	0.425	0.419	0.51	0.527	0.691	0.684
	Qs	1.066	1.068	1.069	1.061	0.35	0.3

* Adapted from

[START_REF] Fahs | The Henry problem: New semianalytical solution for velocity-dependent dispersion[END_REF]

:

Table 6 .

 6 Nusselt number achieved for different number of DOFs per each variable

				DG		COMSOL-Multiphysics	Semi-analytical
		Order 1	Order 2	Order 3	Order 4	Order 1	Order 1
	Number of elements Number of DOFs	10, 000 60, 000	10, 000 120, 000	10, 000 200, 000	10, 000 300, 000	100, 110 201, 022	655, 752 1, 313, 506	--
	Nu	33.14	50.45	45.24	46.15	38.07	44.9	46.14

( t) t (w h ; q * h ) = -(q * h • n Γ )c in , w h F in h .(22b)Here, σ t Γ ≥ 0 corresponds to the penalty parameter for the advective-dispersive problem and we assume that it can decomposed into a dispersive and advective parts,σ t Γ = 2τ - Γ τ +

(hp-version) to even more reduce the computational cost. The proposed method is generic and can be extended to 3D, which is an ongoing research project of this work.

Acknowledgement

The authors would like to appreciate the use of the computational clusters of the High Performance Computing Center (ShahreKord University, Iran), to complete this work. The data used in this work are available on the GitHub repository: https://github.com/fahs-LHYGES


MD Antonio Cubisino
email: antoniocubisino@gmail.com

MD ; F Guillon

MD ; J M Fabre

Laparoscopic right partial adrenalectomy (with video) Titre: Surrénalectomie partielle droite laparoscopique (avec vidéo

TEXT:

The laparoscopic approach has become the gold standard approach for the management of most adrenal gland pathologies. In 1992, Gagner et al. described the first laparoscopic adrenalectomy, and further reports showed how minimally invasive access was able to reduce complications and hospital stay [START_REF] Gagner | Laparoscopic adrenalectomy in Cushing's syndrome and pheochromocytoma[END_REF]. Partial adrenalectomy is an alternative to total adrenalectomy when normal adrenal tissue has to be preserved in order to avoid side effects of adrenal insufficiency. In selected cases such as patients with small benign functional adrenal lesions, small pheochromocytomas, patients with hereditary disease predisposing to multiple adrenal masses, and in patients with prior contralateral adrenalectomy or at high risk for developing contralateral lesions, partial adrenalectomy is a valid alternative to total adrenalectomy. Tumoral diameter (< 3cm) and tumor localization within the adrenal gland are also two criteria impacting the feasibility of partial adrenalectomy. The main goal of partial adrenalectomy is the reduction of patients who would become lifetime dependent of hormonal replacement, if not performed. This video shows a laparoscopic right partial adrenalectomy in a fifty-eight year-old patient with hypertension, muscular weakness and episodic headaches due to primary hyperaldosteronism caused by aldosterone-producing adenoma (Conn's Syndrome).

Biological workup revealed a severe hypokalemia, (K + 2.6 mmol/L) with an aldosterone/renin ratio of 307 ng/dL per ng/(mL/hour), well above normal upper limit value.

An abdominal MRI showed a 18mm nodule of the right adrenal gland and a CT scan confirmed the presence of an enhancing nodule of 9 HU measuring 19mm x 15mm on the right adrenal gland. No lesions were found on the left adrenal gland. An adrenal venous sampling confirmed an unilateral abnormal secretion of aldosterone (2102,4 ng/mL) confirming the diagnosis of unilateral adrenal adenoma. In consideration of the morphological characteristics of the adrenal lesion (lateral localization and <2cm) and of the advantages associated with an adrenal-sparing surgery, we decided to perform a laparoscopic partial adrenalectomy in this patient. This video showed different steps necessary to perform a right partial laparoscopic adrenalectomy [START_REF] Nagaraja | Recurrence and functional outcomes of partial adrenalectomy: a systematic review and meta-analysis[END_REF]. This didactic video will be useful for all surgeons performing adrenal surgical procedures [START_REF]Pheochromocytoma and paraganglioma: an endocrine society clinical practice guideline[END_REF].

Disclosure of interest

The authors declare that they have no competing interest.

Figure 1 .

 1 Figure 1. Patient positioning and port placement for the right partial adrenalectomy.

Figure 2 .

 2 Figure 2. Incision of the posterior peritoneum.

Figure 3 .

 3 Figure 3. Right adrenal gland section.

Figure 4 .

 4 Figure 4. Analysis of the surgical specimen.