L Nouveau
email: leo.nouveau@duke.edu

M Ricchiuto
email: mario.ricchiuto@inria.fr

G Scovazzi

High-Order Gradients with the Shifted Boundary Method: An Embedded Enriched Mixed Formulation for Elliptic PDEs

Keywords: Darcy flow, embedded boundary, finite element method, high-order approximation, stabilized methods, computational fluid dynamics

We propose an extension of the embedded boundary method known as "shifted boundary method" to elliptic diffusion equations in mixed form (e.g., Darcy flow, heat diffusion problems with rough coefficients, etc.). Our aim is to obtain an improved formulation that, for linear finite elements, is at least second-order accurate for both flux and primary variable, when either Dirichlet or Neumann boundary conditions are applied. Following previous work of Nishikawa and Mazaheri in the context of residual distribution methods, we consider the mixed form of the diffusion equation (i.e., with Darcy-type operators), and introduce an enrichment of the primary variable. This enrichment is obtained exploiting the relation between the primary variable and the flux variable, which is explicitly available at nodes in the mixed formulation. The proposed enrichment mimics a formally quadratic pressure approximation, although only nodal unknowns are stored, similar to a linear finite element approximation. We consider both continuous and discontinuous finite element approximations and present two approaches: a non-symmetric enrichment, which, as in the original references, only improves the consistency of the overall method; and a symmetric enrichment, which enables a full error analysis in the classical finite element context. Combined with the shifted boundary method, these two approaches are extended to high-order embedded computations, and enable the approximation of both primary and flux (gradient) variables with second-order accuracy, independently on the type of boundary conditions applied. We also show that the the primary variable is third-order accurate, when pure Dirichlet boundary conditions are embedded.

Introduction

The proposed work stems from the Shifted Boundary Method (SBM), an embedded method that has been recently developed and applied in the context of the Laplace and Stokes problems [START_REF] Main | The shifted boundary method for embedded domain computations. part i: Poisson and stokes problems[END_REF], the Navier-Stokes equations [START_REF] Main | The shifted boundary method for embedded domain computations. part ii: Linear advection-diffusion and incompressible navier-stokes equations[END_REF] and hyperbolic systems [START_REF] Song | The shifted boundary method for hyperbolic systems: Embedded domain computations of linear waves and shallow water flows[END_REF]. The SBM leverages a surrogate/approximate interface representation, by which the location where boundary conditions are applied is shifted from the true to the surrogate interface, and, most importantly, a Taylor expansion is used to modify the value of boundary conditions, with the goal of preventing a reduction in the convergence rates of the overall formulation. The appropriate (modified) boundary conditions are applied weakly, using, for example, a Nitsche strategy (other weak enforcement strategies are, in principle, also possible).

Immersed and embedded boundary methods present several advantages with respect to conformal methods especially regarding mesh generation. If complex geometries are involved, generating the mesh can require considerable time and effort compared to the whole simulation process. Introduced by Peskin in 1972 [START_REF] Charles S Peskin | Flow patterns around heart valves: A numerical method[END_REF], the Immersed Boundary Method (IBM) proposes to mesh the entire computational domain, independently of the geometry of the problem. Bodies or boundaries present in the domain are described implicitly in the mesh by means of a distance function, or similar indicator function. The main idea behind the IBM is to solve the model equations on the entire domain, and to impose the boundary conditionss via a forcing term. The method was originally designed on Cartesian grids and its development has been an active field of research for some decades now. Two exhaustive reviews were written by Mittal and Iaccarino in 2005 [START_REF] Mittal | Immersed boundary methods[END_REF], and Sotiropoulos and Yang in 2014 [START_REF] Sotiropoulos | Immersed boundary methods for simulating fluid-structure interaction[END_REF]. The main drawback of the IBM is the accuracy of boundary conditions, which are most often only first-order accurate. Some strategies have been proposed in the past to increase the IBM accuracy, at the expense of more involved implementations [START_REF] Griffith | On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems[END_REF][START_REF] Griffith | An adaptive, formally second order accurate version of the immersed boundary method[END_REF][START_REF] Leveque | Immersed interface methods for stokes flow with elastic boundaries or surface tension[END_REF][START_REF] Linnick | A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains[END_REF], or to compensate low-order accuracy by means of mesh adaptation [START_REF] Nouveau | An adaptive, residual based, splitting approach for the penalized navier stokes equations[END_REF][START_REF] Billon | Anisotropic boundary layer mesh generation for immersed complex geometries[END_REF][START_REF] Hachem | Immersed stress method for fluid-structure interaction using anisotropic mesh adaptation[END_REF].

A similar objective is pursued by Embedded Boundary Methods, which still use an immersed description of bodies and boundaries, but solve the model equations only in the regions of interest. Within a finite element context, the cut-cell and cutFEM approach [START_REF] Parvizian | Finite cell method[END_REF][START_REF] Düster | The finite cell method for three-dimensional problems of solid mechanics[END_REF][START_REF] Hansbo | An unfitted finite element method, based on nitsche's method, for elliptic interface problems[END_REF][START_REF] Clarke | Euler calculations for multielement airfoils using cartesian grids[END_REF][START_REF] Kirkpatrick | A representation of curved boundaries for the solution of the navier-stokes equations on a staggered three-dimensional cartesian grid[END_REF][START_REF] Udaykumar | Elafint: A mixed eulerian-lagrangian method for fluid flows with complex and moving boundaries[END_REF] are most commonly utilized. These techniques usually combine a weak enforcement of the boundary conditions with a XFEM strategy. For example, in CutFEMs, the boundary is reconstructed as the intersection between the boundary surface and the embedding (underlaying) grid. Approaches of this kind often require complex implementations and may suffer from the so-called small-cut cell issue, unless specific stabilization mechanisms are put in place [START_REF] Burman | Ghost penalty[END_REF][START_REF] Burman | Fictitious domain methods using cut elements: Iii. a stabilized nitsche method for stokes' problem[END_REF][START_REF] Schott | A new face-oriented stabilized xfem approach for 2d and 3d incompressible navier-stokes equations[END_REF][START_REF] Schott | A face-oriented stabilized nitsche-type extended variational multiscale method for incompressible two-phase flow[END_REF]. In fact, the intersection of the grid with the surface boundary may produce elements of extremely small sizes, with negative consequences on the stability of the variational formulation and the magnitude of the condition number of the associated algebraic system.

The SBM is effective in overcoming all these issues, and proved robust and accurate in the case of Dirichlet conditions, while in the case of Neumann conditions, a standard SBM implementation may not be optimally accurate, since the Taylor expansion of the gradient of the solution is required, and this is not feasible for the simple linear finite element approximation. In [START_REF] Main | The shifted boundary method for embedded domain computations. part i: Poisson and stokes problems[END_REF], for the Laplace problem, the use of a reconstructed gradient has been used to overcome this difficulty.

A more robust strategy to achieve accurate gradients is pursued here, and consists in developing a mixed formulation, in which the Poisson operator is replaced by the Darcy operator, and to design advanced schemes for high-order flux approximations. In this context, we would like to mention recent contributions [START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme -i: Basic principle[END_REF][START_REF] Mazaheri | Improved second-order hyperbolic residual-distribution scheme and its extension to third-order on arbitrary triangular grids[END_REF][START_REF] Mazaheri | Efficient high-order discontinuous galerkin schemes with first-order hyperbolic advection-diffusion system approach[END_REF][START_REF] Mazaheri | A first-order hyperbolic system approach for dispersion[END_REF] for body-fitted computations, which were inspired by the earlier work of Caraeni [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for navier-stokes simulations[END_REF]. In these works, mixed forms of the Poisson, diffusion, or advection diffusion operators were proposed, in which the solution derivatives are adopted as main unknowns, together with the solution itself. This procedure enhances the accuracy of the solution, and allows to obtain uniformly second-order accurate derivatives of the solution. An additional pressure enrichment was also proposed in these contributions to improve the order of accuracy of the solution and the associated gradients.

The context and scope of the present work is within this idea, which is applied to mixed formulations of the Poisson problem. Mixed formulations of the Poisson problem (e.g., the Darcy flow equations) encompass a wide range of applications in the areas of thermo-mechanics, heat transfer, geotechnical and petroleum engineering. We will consider both continuous Galerkin (CG) and discontinuous Galerkin (DG) formulations, to demonstrate the flexibility of the proposed approach. For the sake of simplicity, we focus on piecewise linear finite element approximations. The main contributions of this paper can be summarized as follows:

1. We elaborate on the enrichment approach used in [START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme -i: Basic principle[END_REF][START_REF] Mazaheri | Improved second-order hyperbolic residual-distribution scheme and its extension to third-order on arbitrary triangular grids[END_REF][START_REF] Mazaheri | Efficient high-order discontinuous galerkin schemes with first-order hyperbolic advection-diffusion system approach[END_REF][START_REF] Mazaheri | A first-order hyperbolic system approach for dispersion[END_REF] in the context of body fitted stabilized finite elements in mixed form. We propose a symmetric enrichment method that can be rigorously proved to be stable and third-/second-order accurate for the primary variable/gradient, respectively.

2. We extend these enriched schemes to embedded computations using the SBM. Higher-order extrapolations are introduced to preserve the same consistency of the enriched formulations, thus providing a genuinely secondorder embedded method for both the primary variable and its gradient, and a third-order method for the primary variable when only Dirichlet boundary conditions are embedded.

3. We provide an extensive numerical investigation of the properties of the methods proposed, confirming the expected orders of accuracy and stability properties. In particular, we show that quadratic exact solutions can be recovered within machine precision with linear elements, and that additional stabilization terms are required with the non-symmetric enrichment, in order to obtain stability properties similar to those of the symmetric enrichment.

In conclusion, the proposed constructions allow to overcome one of the limitations of the original shifted boundary method, providing a genuinely high-order formulation.

The present article is organized as follows: the strong formulation of the problem of interest is presented in Section 2; the variational formulations for conformal grids are given in Section 3; the shifted boundary method is presented in Section 4; Section 5 is devoted to the primary variable enrichment and the construction of numerical schemes of improved accuracy; finally, an extensive validation of the proposed methods is detailed in Section 6; and conclusions are summarized in Section 7.

Problem statement: Darcy flow

Darcy's flow equations provide a homogenized macroscopic transport model through porous media. The Darcy operator can also be interpreted as a mixed form of the Poisson operator, in which the gradient (i.e., the flux) is solved as a variable together with the primary variable (i.e., the pressure). Darcy-type operators are very important in a number of fields that do not necessarily involve flow in porous media, such as heat transfer for systems with rough distributions of conductivities, etc. The Darcy flow problem is governed by the following set of equations:

                   Λ -1 β + ∇p = 0 in Ω ∇ • β = φ in Ω p = p D on Γ D β • n = h N on Γ N (1)
where β denotes the flux, p the pressure, Λ the permeability tensor and φ a source term. Note that the permeability is in general spatially dependent, that is Λ = Λ(x), and may even be discontinuous. The boundary of the computational domain Ω is denoted ∂Ω, and partitioned into Γ D and Γ N , on which Dirichlet and Neumann boundary conditions are imposed, respectively. In particular, we have the following non-overlapping decomposition of the boundary:

∂Ω = Γ N ∪ Γ D , Γ N ∩ Γ D = ∅.

Conformal mixed variational formulations

We start our discussion by introducing the base variational formulations on conformal (i.e., body-fitted) grids, and for that purpose, by introducing some general notation.

Generalities

Let T h denote a tessellation of the domain Ω, composed of non-overlapping regular elements K, with T h = ∪K. For a given element K, we denote by |K| its area, by ∂K its boundary, by γ K a generic face/edge, and by h K its reference length. We denote the reference mesh size by h = max K∈T h h K . The set of internal edges is denoted by E i , and for an interior edge γ K , K + and K -are the elements on each side of it. For a given face/edge γ K of size |γ K |, the characteristic length in a direction perpendicular to it is denoted by h ⊥ , and computed as:

h ⊥ =                |K + | + |K -| 2|γ K | if γ K ∩ ∂Ω = ∅ |K| 2|γ K | if γ K ∩ ∂Ω ∅ (2)
The mesh is assumed to verify classical regularity assumptions. In particular it is assumed that there exist bounded constants

C 1 , C 2 ∈ R + such that: C 1 h ≤ h ⊥ ≤ C 2 h (3)
Let us now introduce the function spaces that will be used to build the discontinuous Galerkin (DG) and continuous Galerkin (CG) schemes. The discontinuous Galerkin spaces are defined as:

S l d (Ω) = p ∈ L 2 (Ω) : p| K ∈ P l (K), ∀K ∈ T h V k d (Ω) = β ∈ L 2 (Ω) d , ∇ • β ∈ L 2 (Ω) : β| K ∈ P k (K) d , ∀K ∈ T h (4)
where P n (K) is the space of polynomial functions of order at most n on K. From these definitions, the continuous Galerkin spaces are defined as

S l c (Ω) = S l d (Ω) ∩ C 0 (Ω) V k c (Ω) = V k d (Ω) ∩ (C 0 (Ω)) d (5)
where C 0 (Ω) denotes the set of continuous functions on Ω. We also introduce the tensor product spaces

M d (Ω) = V k d (Ω) × S l d (Ω) and M c (Ω) = V k c (Ω) × S l c (Ω)
. The L 2 (K) inner products over an element and its boundary are defined as:

(v, w) K = K vw and (v, w) K = K v • w u, v ∂K - ∂K uv and u, v ∂K = K u • v (6)
Next, we introduce the average and jump operators {•} and [[•]] for discontinuous scalar and vector fields as follows:

{p} = 1 2 (p + + p -) [[p]] = p + n + + p -n - {β} = 1 2 (β + + β -) [[β]] = β + • n + + β -• n - (7)
where the + andsigns indicate the two sides of an edge/face. By a reordering argument, the following useful identity can be proved:

K w, v • n ∂K = {w}, [[v]] E i + [[w]], {v} E i + w, v • n ∂Ω (8)

Continuous and discontinuous Galerkin mixed formulations

We briefly recall here the construction of discontinuous Galerkin scheme and refer to [START_REF] Cockburn | Discontinuous Galerkin Methods: Theory, Computation and Applications[END_REF][START_REF] Arnold | Unified analysis of discontinuous galerkin methods for elliptic problems[END_REF][START_REF] Douglas | Interior penalty procedures for elliptic and parabolic galerkin methods[END_REF] and references therein for more details. The corresponding continuous formulations will be obtained by enforcing continuity of the discrete approximations across element interfaces.

Discontinuous Galerkin schemes have undergone considerable development from early works [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF][START_REF] Douglas | Interior penalty procedures for elliptic and parabolic galerkin methods[END_REF][START_REF] Cockburn | Discontinuous Galerkin Methods: Theory, Computation and Applications[END_REF] and allow for a robust treatment of rough distributions of permeabilities/conductivities [START_REF] Arnold | Unified analysis of discontinuous galerkin methods for elliptic problems[END_REF][START_REF] Douglas | Interior penalty procedures for elliptic and parabolic galerkin methods[END_REF][START_REF] Huang | A high-order, fully coupled, upwind, compact discontinuous galerkin method for modeling of viscous fingering in compressible porous media[END_REF][START_REF] Scovazzi | A discontinuous galerkin method for gravity-driven viscous fingering instabilities in porous media[END_REF][START_REF] Douglas | Self-adaptive finite element simulation of miscible displacement in porous media[END_REF]. The main drawback of DG methods is the large number of degrees-of-freedom, which relates to the overall computational cost and has spurred the development of the hybridizable DG method [START_REF] Cockburn | A hybridizable discontinuous galerkin method for steady-state convectiondiffusion-reaction problems[END_REF], the enriched Galerkin approach [START_REF] Lee | A locally conservative enriched galerkin approximation and efficient solver for elliptic and parabolic problems[END_REF], or the variational multiscale DG approach [START_REF] Thomas | A multiscale discontinuous galerkin method with the computational structure of a continuous galerkin method[END_REF][START_REF] Wang | Dual-scale galerkin methods for darcy flow[END_REF]). These methods are generally more complex to implement, and for this reason we also consider stabilized CG methods [START_REF] Masud | A stabilized mixed finite element method for darcy flow[END_REF], which may become more appropriate and efficient when the variation of the permeabilities/conductivities is smooth.

Specifically, we consider the stabilized CG and DG schemes of Masud et. al. [START_REF] Masud | A stabilized mixed finite element method for darcy flow[END_REF][START_REF] Thomas | A stabilized mixed discontinuous galerkin method for darcy flow[END_REF] and the interior penalty DG scheme of [START_REF] Douglas | Interior penalty procedures for elliptic and parabolic galerkin methods[END_REF][START_REF] Arnold | Unified analysis of discontinuous galerkin methods for elliptic problems[END_REF]. As a starting point, each equation of the Darcy problem (1) is multiplied by the test function pair (w, q) ∈ M d (Ω) and the result is integrated over a generic element:

(w, Λ -1 β) K + (w, ∇p) K = (w, Λ -1 β) K -(∇ • w, p) K + w • n, p ∂K = 0 (q, ∇ • β) K = -(∇q, β) K + q, β • n ∂K = (q, φ) K (9)
The boundary traces are then replaced by the numerical fluxes p, β and an integration by parts is performed again on the second equation before summing on all the elements of the mesh:

(w, Λ -1 β) Ω -(∇ • w, p) Ω + K w • n, p ∂K = 0 (q, ∇ • β) Ω + K q, ( β -β) • n ∂K = (q, φ) Ω (10)
The numerical fluxes link the discontinuous polynomial approximations of the solution between neighboring elements. Using the definition of average and jump operators along with equality (8), the set of equations [START_REF] Clarke | Euler calculations for multielement airfoils using cartesian grids[END_REF] becomes:

(w, Λ -1 β) Ω -(∇ • w, p) Ω + [[w]], { p} E i + {w}, [[ p]] E i + w • n, p ∂Ω = 0 (q, ∇ • β) Ω + [[q]], { β -β} E i + {q}, [[ β -β]] E i + q, ( β -β) • n ∂Ω = (q, φ) Ω (11)
The definition of the numerical fluxes is what differentiates the many possible DG formulations. In this work, we use the Interior Penalty (IP) method [START_REF] Douglas | Interior penalty procedures for elliptic and parabolic galerkin methods[END_REF][START_REF] Fanett | An elliptic collocation-finite element method with interior penalties[END_REF][START_REF] Douglas | An interior penalty finite element method with discontinuous elements[END_REF][START_REF] Arnold | Unified analysis of discontinuous galerkin methods for elliptic problems[END_REF], and the numerical fluxes are defined as

             p = {p} + α p [[β]], β = {β} + α β [[p]], on E i p = p D , β = β + α β (p -p D )n, on Γ D p = p, β • n = h N , on Γ N ( 12
)
where α p = αp h ⊥ /||Λ|| and α β = αβ /(h ⊥ )||Λ|| are penalty coefficients, with αp and αβ real positive constants. As proposed in [START_REF] Thomas | A stabilized mixed discontinuous galerkin method for darcy flow[END_REF], a stabilizing term similar to the one employed in the continuous formulation [START_REF] Masud | A stabilized mixed finite element method for darcy flow[END_REF] can be added and the final formulation yields:

B d (w, q; β, p) = B Gal d (w, q; β, p) + B S tab (w, q; β, p) = L Gal d (w, q) + L S tab (w, q) B Gal d (w, q; β, p) = (w, Λ -1 β) Ω -(∇ • w, p) Ω + (q, ∇ • β) Ω + [[w]], {p} E i -{q}, [[β]] E i + [[q]], α β [[p]] E i + [[w]], α p [[β]] E i + w • n, p Γ N -q, β • n Γ N + q, α β p Γ D B S tab (w, q; β, p) = 1 2 -Λ -1 w + ∇q, β + Λ∇p Ω + ζ div 2 ∇ • w, ||Λ||h 2 ∇ • β Ω L Gal d (w, q) = (q, φ) Ω -w • n, p D Γ D + q, α β p D Γ D -q, h N Γ N L S tab (w, q) = ζ div 2 (∇ • w, ||Λ||h 2 φ) Ω (13)
where ζ div = O(1) is a stabilization constant. In the stabilization B S tab (w, q; β, p), the first term containing the residual of the momentum equation β + Λ∇p will be referred as "momentum stabilization" while the term involving the divergence of the flux will be referred as "div-div" stabilization. From [START_REF] Douglas | Interior penalty procedures for elliptic and parabolic galerkin methods[END_REF], the continuous stabilized formulation can be obtained by removing terms containing jumps across internal edges/faces, and by replacing averages of the unknowns and test functions by their unique continuous values. The stability of the proposed formulation was proved in [START_REF] Thomas | A stabilized mixed discontinuous galerkin method for darcy flow[END_REF][START_REF] Masud | A stabilized mixed finite element method for darcy flow[END_REF] in the norm:

|||(β, p)||| 2 = 1 2 (||Λ -1/2 β|| 2 Ω + ||Λ 1/2 ∇p|| 2 Ω ) + ζ div 2 ||Λ 1/2 h∇ • β|| 2 Ω + ||( αβ ||Λ||/h ⊥ ) 1/2 p|| 2 Γ D (14)
+ ||( αβ ||Λ||/h ⊥ ) 1/2 [[p]]|| 2 E i + ||( αp h ⊥ /||Λ||) 1/2 [[β]]|| 2 E i (15)
with the following error estimate:

|||β e -β, p e -p||| ≤ C 1 h k+1 |β| k+1 + C 2 h l |p| l+1 (16)
where C 1 and C 2 are two constants. Hence this method converges with order r = min(k + 1, l) for the flux. L 2 -error estimates derived in [START_REF] Badia | Stabilized continuous and discontinuous galerkin techniques for darcy flow[END_REF] show that the pressure converges with order r + 1. For linear approximations, k = l = 1, meaning that the pressure is second-order accurate and the flux is first-order accurate, as will also be verified in Section 6.

Remark 1. The definition of the numerical fluxes [START_REF] Cockburn | Discontinuous Galerkin Methods: Theory, Computation and Applications[END_REF] can be related to the variational multiscale approach proposed by Badia and Codina [START_REF] Badia | Stabilized continuous and discontinuous galerkin techniques for darcy flow[END_REF].

the closest to the true boundary Γ, in some sense. Γ can be constructed, for example, by computing the intersections of the grid and the true boundary Γ and using closest-point projection algorithms to detect the closest face/edge of Γ to Γ. Other choices are of course possible, such as level-set functions, for instance, as long as the overall topology of Γ and Γ are close to each other, that is if Γ has a certain number of holes, the same number of holes needs to be also present in Γ. The surrogate boundary Γ encloses the surrogate domain Ω. In particular, ñ indicates the unit outward normal of Γ, to be distinguished from the outward normal n of Γ. The mapping from the surrogate to true interface is then defined as:

M : Γ → Γ x → x (17)
In particular, the map M can be characterized through a distance vector function:

d M ( x) = x -x = [M -I]( x) (18)
In the following, we will simply define the distance vector function as d, and note the following relation with the normal of the true boundary:

n = d ||d|| (19)
We recall here an assumption of sufficient geometric resolution, discussed in detail in the original development of the SBM [START_REF] Main | The shifted boundary method for embedded domain computations. part i: Poisson and stokes problems[END_REF]:

Assumption 1.
The vector distance d is defined as d = ||d||n, where the normal n to the true boundary and the normal ñ to the surrogate boundary satisfy:

n • ñ > 0 (20)
Remark 2. Condition ( 20) is a resolution condition on the grid used in SB simulations. In practice, a check is performed at the beginning of every simulation to verify that (20) is satisfied. If this is not the case, the grid is refined (for example, splitting every edge of the grid into two), and ( 20) is checked again. This simple strategy proved extremely robust in ensuring the well-posedness of the overall SB approach. In principle, more sophisticated edge swapping strategies can be used to ensure the grid satisfies (20), although we did not perform numerical experiments to confirm these claims. There is also the more general question of the existence of the distance function between the surrogate and true boundary. This is a more fundamental issue common to all immersed/embedded methods and not just the SBM. For a more detailed view on the minimal conditions of existence, the reader is prompted to read [START_REF] Rangarajan | Universal meshes: A method for triangulating planar curved domains immersed in nonconforming meshes[END_REF]. Condition [START_REF] Hansbo | An unfitted finite element method, based on nitsche's method, for elliptic interface problems[END_REF] is somewhat related to the conditions specified in [START_REF] Rangarajan | Universal meshes: A method for triangulating planar curved domains immersed in nonconforming meshes[END_REF], in the sense that it is a sufficient condition (although possibly not necessary) for the well-posedness of the existence of a distance function.

Through the map M, it is possible to define the extension ψ on Γ of a function ψ originally defined on Γ, as

ψ( x) ≡ ψ(M( x)) (21)
For instance, the unit normal n and tangential vectors τ i ( 1 < i < n d -1 ) of the true boundary Γ can be extended to the boundary Γ as follows:

n( x) ≡ n(M( x)) τi ( x) ≡ τ i (M( x)) (22)
In what follows, and for the sake of simplicity, we will omit the bar in the expressions of the extensions of the normal and tangent vectors whenever there is no possibility of confusion, thus writing n( x) and τ i ( x). We can also introduce the derivatives in the directions n and τ i of a function ψ at a point x ∈ Γ: These constructions are particularly useful if we consider the solution u to a partial differential equation and, for example, we desire to extend the Dirichlet boundary condition u D from the true boundary Γ D to a surrogate boundary ΓD . In particular, the following Taylor expansion centered at x ∈ ΓD holds for

ψ ,n ( x) = ∇ψ( x) • n( x) ψ ,τ i ( x) = ∇ψ( x) • τ i ( x) (23)
x = M( x) = x + d( x) ∈ Γ D : u D (x) = u D ( x + d( x)) = u D ( x) + ∇u( x) • d( x) + 1 2 d T ( x)H (u)d( x) + O(||d( x)|| 3 ) (24)
where H (u) denotes the Hessian of u. This last expression can be used to develop a new strategy for the imposition of Dirichlet boundary conditions in the context of embedded methods. Depending on the desired accuracy, additional terms can be considered in the Taylor expansion.

The shifted boundary variational formulation

We now propose to extend the scheme (13) to an embedded formulation, using the SBM [START_REF] Main | The shifted boundary method for embedded domain computations. part i: Poisson and stokes problems[END_REF][START_REF] Main | The shifted boundary method for embedded domain computations. part ii: Linear advection-diffusion and incompressible navier-stokes equations[END_REF]. Looking at equations [START_REF] Cockburn | A hybridizable discontinuous galerkin method for steady-state convectiondiffusion-reaction problems[END_REF], the integrations over internal edges are not impacted by the integration on the surrogate domain. However, the boundary integral on the surrogate domain ∂ Ω now involves the surrogate normal ñ. Thus, equations [START_REF] Cockburn | A hybridizable discontinuous galerkin method for steady-state convectiondiffusion-reaction problems[END_REF] are written as

(w, Λ -1 β) Ω -(∇ • w, p) Ω + [[w]], { p} E i + {w}, [[ p]] E i + w • ñ, p ∂ Ω = 0 (25) (q, ∇ • β) Ω + [[q]], { β -β} E i + {q}, [[ β -β]] E i + q, ( β -β) • ñ ∂ Ω = (q, φ) Ω ( 26
)
The definitions of the numerical fluxes [START_REF] Cockburn | Discontinuous Galerkin Methods: Theory, Computation and Applications[END_REF] need to be modified accordingly on the boundaries. The main differences arise for Neumann boundary conditions, which involve the flux component along the normal to the true boundary. The first step is to decompose the surrogate normal ñ into the normal n and tangents τ i to the true boundary. Then, the boundary term in equation ( 26) yields, on Neumann boundaries,

q, ( β -β) • ñ ΓN = q(n • ñ), ( β -β) • n ΓN + q(τ i • ñ), ( β -β) • τ i ΓN (27)
In addition, we impose the following condition on the tangential flux:

β • τ i = β • τ i , on ΓN (28)
Then, omitting the stabilization terms, the scheme (13) reads:

(w, Λ -1 β) Ω -(∇ • w, p) Ω + (q, ∇ • β) Ω + [[w]], {p} E i -{q}, [[β]] E i + [[q]], α β [[p]] E i + [[w]], α p [[β]] E i + w • ñ, p ΓN -q(n • ñ), β • n -hN ) ΓN + w • ñ, pD ΓD + q, α β (p -pD ) ΓD = 0 (29)
Here we wrote pD and hN to indicate the fact that these values are not equal to p D and h N , which are instead defined on the true boundaries Γ D and Γ N . pD and hN need to be defined appropriately to preserve the accuracy of the schemes.

Shifted Dirichlet boundary conditions

Following the methodology developed in [START_REF] Main | The shifted boundary method for embedded domain computations. part i: Poisson and stokes problems[END_REF], with the theory and notations proposed in Section 4, the Taylor expansion ( 24) is employed to extend the known function p D from the true boundaries Γ D onto the surrogate ones ΓD :

p D (x) = pD ( x) = pD ( x) + ∇p • d + O(||d|| 2 ) (30)
Thus, on Dirichlet boundaries, pD ( x) = p D -∇p • d needs to be substituted in [START_REF] Cockburn | Discontinuous Galerkin Methods: Theory, Computation and Applications[END_REF] in place of p D .

Shifted Neumann boundary condition

Similar to the case of Dirichlet boundary conditions, on Neumann boundaries we have:

h N (x) = hN ( x) = hN ( x) + [(∇β)d] • n + O(||d|| 2 ) (31) so that hN ( x) = h N (x) -[(∇β)d]
• n needs to be substituted in [START_REF] Cockburn | Discontinuous Galerkin Methods: Theory, Computation and Applications[END_REF] in place of h N .

Shifted Boundary Method scheme

Using now the corrected (shifted) values of the boundary conditions ( 30)-( 31), the final discontinuous weak formulation reads:

B d (w, q; β, p) = B Gal d (w, q; β, p) + B S tab (w, q; β, p) = L Gal d (w, q) + L S tab (w, q) B Gal d (w, q; β, p) = (w, Λ -1 β) Ω -(∇ • w, p) Ω + (q, ∇ • β) Ω + [[w]], {p} E i -{q}, [[β]] E i + [[q]], α β [[p]] E i + [[w]], α p [[β]] E i + w • ñ, p ΓN -q(n • ñ), (β + [∇β]d) • n ΓN -w • ñ, ∇p • d ΓD + q + ∇q • d, α β (p + ∇p • d) ΓD L Gal d (w, q) = (q, φ) Ω -w • ñ, pD ΓD + q + ∇q • d, α pD ΓD -q(n • ñ), hN ΓN (32)
Once again, the continuous stabilized variational formulation is obtained from [START_REF] Mazaheri | Efficient high-order discontinuous galerkin schemes with first-order hyperbolic advection-diffusion system approach[END_REF] by neglecting all the internal jump terms:

B c (w, q; β, p) = B Gal c (w, q; β, p) + B S tab (w, q; β, q) = L Gal c (w, q) + L S tab (w, q) B Gal c (w, q; β, p) = (w,

Λ -1 β) Ω -(∇ • w, p) Ω + (q, ∇ • β) Ω + w • ñ, p ΓN -q(n • ñ), (β + [∇β]d) • n ΓN + w • ñ, ∇p • d ΓD + q + ∇q • d, α β (p + ∇p • d) ΓD L Gal c (w, q) = (q, φ) Ω -w • ñ, pD ΓD + q + ∇q • d, α β pD ΓD -q(n • ñ), hN ΓN (33)
From this continuous Galerkin formulation [START_REF] Mazaheri | A first-order hyperbolic system approach for dispersion[END_REF], the Euler-Lagrange equations can be obtained by performing integration by parts. The resulting expression is:

0 = (w, (Λ -1 β + ∇p)) Ω + (q, (∇ • β -φ)) Ω -w • ñ + α β (q + ∇q • d), p + ∇p • d -pD ΓD -q(n • ñ), (β + (∇β)d) • n -hN ΓN (34)
The Euler-Lagrange equations identify the consistency conditions enforced by the SBM: the partial differential equation is enforced on the interior of the surrogate domain Ω; Dirichlet conditions for the pressure are approximately enforced by means of the Taylor expansion, in a Nitsche fashion [START_REF] Nitsche | Uber ein Variationsprinzip zur Losung Dirichlet-Problemen bei Verwendung von Teilraumen, die keinen Randbedingungen unteworfen sind[END_REF]; and Neumann boundary conditions are also approximately enforced using the Taylor expansion. As illustrated later in Section 6, the formulation (32) yields a first-order accurate flux, and a second-order accurate pressure for Dirichlet boundary conditions. The accuracy of the latter reduces to first order when embedded Neumann conditions are considered.

The next section is devoted to improving the formulation just described, so that second-order accuracy can be obtained for both pressure and flux, regardless of the type of embedded boundary condition.

Pressure enrichment and high-order flux finite element scheme

We propose next to modify the variational formulations ( 13) and [START_REF] Mazaheri | Efficient high-order discontinuous galerkin schemes with first-order hyperbolic advection-diffusion system approach[END_REF] with the aim of obtaining at least secondorder accuracy on both pressure and flux for any kind of embedded boundary conditions. Taking inspiration from Caraeni [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for navier-stokes simulations[END_REF], Nishikawa [START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme -i: Basic principle[END_REF] and Mazaheri [START_REF] Mazaheri | Improved second-order hyperbolic residual-distribution scheme and its extension to third-order on arbitrary triangular grids[END_REF] improved the accuracy of advection-diffusion problems by using the mixed form of the equations. In particular, by exploiting the relations between the flux and the pressure gradients in the Darcy equation, we can construct a quadratic pressure approximation, as discussed in the following section.

Enriched pressure approximation

Our aim is to to define a quadratic pressure polynomial in each element. One way to achieve this is to add additional degrees-of-freedom associated to the pressure values at edge midpoints, as illustrated in Figure 2 (for simplicity, in two dimensions). These additional degrees-of-freedom become unknowns of the problem and thus the size of the algebraic system of equations to be solved increases. The resulting approximation u h of the variable u reads:

u h (x) = i,node u i ϕ P 2 i (x) + j,edge u j ϕ P 2 j (x) (35)
where the u i and u j correspond to solution values associated to the nodal degrees-of-freedom (vertices), and to the extra P 2 degrees-of-freedom (edge mid-points). In addition, we denote by ϕ P 1 i and ϕ P 2 i the linear and quadratic local shape functions.

Both in two and three space dimensions one can easily prove the following relations, which will be useful later on (cf. e.g. [START_REF] Ern | Theory and practice of finite elements[END_REF]):

           ϕ P 2 i = ϕ P 1 i (2ϕ P 1 i -1) i is a vertex ϕ P 2 k = 4ϕ P 1 i ϕ P 1 j k is the midpoint of edge [i, j] (36)
The idea is now to replace the mid-point unknowns by values appropriately defined using the known nodal values of the pressure and its gradient. For this purpose, Taylor expansions are used in combination with the mixed form of the equations. Looking at Figure 2, focusing on and edge [i, j], the aim is to define the pressure at the midpoint l using only information stored at nodes i and j. We thus consider a third-order truncated Taylor expansion of the pressure along the edge, as well as a second-order expansion of the flux:

p(x + δx) =p(x) + ∇p(x) • δx + 1 2 δx T H (p)δx + O(||δx|| 3 ) β(x + δx) =β(x) + (∇β(x))δx + O(||δx|| 2 ) ( 37
)
where H denotes the Hessian of the pressure.

Remark 3. In this section, the use of the Taylor expansion is not related to the SBM. In equation [START_REF] Nitsche | Uber ein Variationsprinzip zur Losung Dirichlet-Problemen bei Verwendung von Teilraumen, die keinen Randbedingungen unteworfen sind[END_REF], δx is defined along the edges of the grid, and not along the vector distance d, introduced in Section 4 to map from the true to the surrogate boundary.

If one considers [START_REF] Nitsche | Uber ein Variationsprinzip zur Losung Dirichlet-Problemen bei Verwendung von Teilraumen, die keinen Randbedingungen unteworfen sind[END_REF] in combination with linear approximation, gradients are constant and the Hessian of the pressure is null. However, using the Darcy equation, the gradient of the pressure can be evaluated as ∇p = -Λ -1 β which is actually linear and thus the Hessian can be now evaluated as

H (p) = -∇ Λ -1 β + ∇ Λ -1 β T (38)
In the last expression, the superscript T denotes the transpose of a tensor. Hence,

δx T H (p)δx = -δx T ∇ Λ -1 β δx
and the Taylor expansions [START_REF] Nitsche | Uber ein Variationsprinzip zur Losung Dirichlet-Problemen bei Verwendung von Teilraumen, die keinen Randbedingungen unteworfen sind[END_REF] can thus be rewritten as:

p(x + δx) ≈ p(x) -Λ -1 β • δx - 1 2 δx T ∇(Λ -1 β)δx Λ -1 β (x + δx) ≈ Λ -1 β (x) + ∇ Λ -1 β (x)δx (39)
Denoting e i→ j = x j -x i the edge vector, the Taylor expansion of Λ -1 β on the entire edge gives:

∇ Λ -1 β (x)e i→ j = Λ -1 j β j -Λ -1 i β i ( 40
)
Remark 4. The error of the the Taylor expansion (37) can be estimated using the Lagrange formula of the remainder.

In particular, classical arguments can be used to show that ∃c ∈ e i→ j such that:

(Λ -1 j β j -Λ -1 i β i ) -∇ Λ -1 β (x)e i→ j = (Λ -1 j β j -Λ -1 i β i ) + H (p)(x)e i→ j = 1 2 × 2! e T i→ j ∇ 2 β(c)e i→ j (41)
Applying the expansion at both node i and j, and, for symmetry reasons, taking the average, we end up with the following estimate for the edge mid-point values of the pressure:

p l = pe + 1 8 (Λ -1 j β j -Λ -1 i β i ) • (x j -x i ) ( 42
)
where pe = 1 2 (p j + p i ) is the average pressure on the edge. Feeding this expression back into [START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme -i: Basic principle[END_REF] and using [START_REF] Nishikawa | On hyperbolic method for diffusion with discontinuous coefficients[END_REF], one obtains, after a few manipulations,

p K (x) = i p i ϕ P 1 i (x) + e 1 2 ∆(Λ -1 β) e • ∆x e ϕ e,P 1 1 (x)ϕ e,P 1 2 (x) (43)
where now only P 1 -shape functions are involved, and ∆t denotes the difference of the nodal values of a given vector t along an edge e. Note that the second sum over the edges of the element K contains the remainder of the P 2correction associated to the product of test functions centered at each node of the edge e.

Remark 5. The definition [START_REF] Schott | A face-oriented stabilized nitsche-type extended variational multiscale method for incompressible two-phase flow[END_REF] of the enrichment is also valid for tetrahedra.

We will now define the new pressure-enriched schemes, in which p * will denote the enriched pressure [START_REF] Schott | A face-oriented stabilized nitsche-type extended variational multiscale method for incompressible two-phase flow[END_REF]. We first propose a formulation akin the ones in [START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme -i: Basic principle[END_REF][START_REF] Nishikawa | On hyperbolic method for diffusion with discontinuous coefficients[END_REF], in which [START_REF] Douglas | Interior penalty procedures for elliptic and parabolic galerkin methods[END_REF][START_REF] Mazaheri | Efficient high-order discontinuous galerkin schemes with first-order hyperbolic advection-diffusion system approach[END_REF] are modified by replacing directly p * for the pressure in the bilinear forms. As numerically shown in Section 6, such formulations recover the desired accuracy once proper stabilization procedures are applied. We will then elaborate on this formulation proposing a symmetric enrichment, where both the trial and test function spaces are modified. This approach yields a natural stabilization mechanism without requiring additional ad-hoc terms in the variational form, and can be elegantly proved to satisfy the desired consistency, stability, and error estimate by simple generalizations of classical arguments [START_REF] Masud | A stabilized mixed finite element method for darcy flow[END_REF][START_REF] Thomas | A multiscale discontinuous galerkin method with the computational structure of a continuous galerkin method[END_REF].

Modified schemes: Non symmetric enrichment

Conformal case

The conformal case is straightforward, in that the enrichment [START_REF] Schott | A face-oriented stabilized nitsche-type extended variational multiscale method for incompressible two-phase flow[END_REF] is used in the variational formulation [START_REF] Douglas | Interior penalty procedures for elliptic and parabolic galerkin methods[END_REF]:

B d (w, q; β, p * ) = B Gal d (w, q; β, p * ) + B S tab (w, q; β, p * ) = L Gal d (w, q) + L S tab (w, q) B Gal d (w, q; β, p * ) = (w, Λ -1 β) Ω -(∇ • w, p * ) Ω + (q, ∇ • β) Ω + [[w]], {p * } i -{q}, [[β]] i + [[q]], α β [[p * ]] i + [[w]], α p [[β]] E i + w • n, p * Γ N -q, β • n Γ N + q, α β p * Γ D L Gal d (w, q) = (q, φ) Ω -w • n, p D Γ D + q, α β p D Γ D -q, h N Γ N (44)
While the test functions employed are still linear, that is (w, q) ∈ S 1 (Ω) × V 1 (Ω), the enriched pressure p * is defined to be a quadratic approximation of the pressure, that is p * ∈ S 2 (Ω). Unfortunately, as the test and trial functions spaces are different, namely (w,

q) ∈ V 1 (Ω) × S 1 (Ω) and (β, p * ) ∈ V 1 × S 2 (Ω)
, a standard coercivity argument cannot be used to derive the stability bound [START_REF] Düster | The finite cell method for three-dimensional problems of solid mechanics[END_REF] and error estimate ( 16), as already seen in Section 3. We can only show that the consistency estimate is improved by the enrichment, as detailed in AppendixA. Numerical evidence detailed here and in [START_REF] Mazaheri | Improved second-order hyperbolic residual-distribution scheme and its extension to third-order on arbitrary triangular grids[END_REF][START_REF] Mazaheri | Efficient high-order discontinuous galerkin schemes with first-order hyperbolic advection-diffusion system approach[END_REF] shows that the non-symmetric enrichment yields a third-order pressure and a second-order flux, and allows to recover, within machine precision, quadratic polynomial solutions of the Darcy flow problem (1). Remark 6. To achieve the levels of accuracy mentioned above, additional stabilization mechanisms are however required in some cases. For example, the "div-div" stabilization is necessary to avoid spurious modes, which may originate from the enhancing terms that replace the mid-point values. The authors of [START_REF] Mazaheri | Improved second-order hyperbolic residual-distribution scheme and its extension to third-order on arbitrary triangular grids[END_REF][START_REF] Mazaheri | Efficient high-order discontinuous galerkin schemes with first-order hyperbolic advection-diffusion system approach[END_REF] propose to include a condition of consistency on the tangent variation of the flux on Dirichlet boundaries, adding to the weak Galerkin formulation the additional boundary penalty term

w • τ i , α τ i Λ -1 β • τ i + ∂ τ i p D Γ D ( 45
)
where τ i denotes the tangent of the boundary in the ith direction, ∂ τ i p D = ∇p D • τ i is the gradient of the Dirichlet value along the boundary, and the penalty coefficient α τ is defined as α τ = αh/||Λ||. Such penalty term has a behavior similar to the "div-div" stabilization term ζ 2 ∇ • w, ||Λ||h 2 ∇ • β Ω in B S tab (w, q; β, p). For this reason, the "div-div" stabilization can be removed if (45) is employed.

Embedded case: third-order extrapolation

In the embedded case, the boundary terms of the type q, α β (pp D ) Γ D + w • n, p D Γ D need to be modified to account for a parabolic extrapolation, in order to be consistent with the enrichment of the pressure. For this reason, when extending the value of the pressure from Γ D to ΓD , the Taylor expansion (30) has been enhanced using the Hessian of the pressure:

p D (x) = pD = p D ( x) + ∇p • d + 1 2 d T H (p)d (46)
As in the case of internal edges, to evaluate the Hessian we use the nodal flux, which essentially leads to [START_REF] Parvizian | Finite cell method[END_REF] with δx = d. The resulting shifted boundary condition reads:

p D (x) = pD = p D ( x) -(Λ -1 β) • d - 1 2 d T ∇ Λ -1 β d (47)
Thus, the Dirichlet terms are now modified as follows: [START_REF] Udaykumar | Elafint: A mixed eulerian-lagrangian method for fluid flows with complex and moving boundaries[END_REF] where the expression in the left slot of the first term is designed to guarantee the symmetry and positive semidefiniteness of the associated variational form, in the spirit of the formulation used in [START_REF] Main | The shifted boundary method for embedded domain computations. part i: Poisson and stokes problems[END_REF].

q-Λ -1 w•d- 1 2 d T ∇ Λ -1 w d, α β (p-Λ -1 β•d- 1 2 d T ∇ Λ -1 β d-pD ) ΓD + w• ñ, Λ -1 β•d+ 1 2 d T ∇ Λ -1 β d+ pD ΓD
Since no modification is performed on the flux itself, the Neumann boundary terms remain unchanged, and the final embedded formulation using pressure enrichment reads:

B d (w, q; β, p * ) = B Gal d (w, q; β, p * ) + B S tab (w, q; β, p * ) = L Gal d (w, q) + L S tab (w, q) B Gal d (w, q; β, p * ) = (w, Λ -1 β) Ω -(∇ • w, p * ) Ω + (q, ∇ • β) Ω + [[w]], {p * } i -{q}, [[β]] i + [[q]], α β [[p * ]] i + [[w]], α p [[β]] E i + q -Λ -1 w • d - 1 2 d T ∇ Λ -1 w d, α β (p -Λ -1 β • d - 1 2 d T ∇ Λ -1 β d) ΓD + w • ñ, Λ -1 β • d + 1 2 d T ∇ Λ -1 β d ΓD + w • ñ, p * ΓN -q(n • ñ), (β + [∇β]d) • n ΓN L Gal d (w, q) = (q, φ) Ω -q(n • ñ), hN ΓN -w • ñ, pD ΓD + q -Λ -1 w • d - 1 2 d T ∇ Λ -1 w d, α β pD ΓD (49)
As before, the continuous Galerkin formulation is obtained by removing all internal edge integrals involving jumps and is omitted for the sake brevity. The Euler-Lagrange equations associated to the enriched scheme read:

0 =(w, (Λ -1 β + ∇p * )) Ω + (q, (∇ • β -φ)) Ω -w • ñ + α β (q -Λ -1 w • d - 1 2 d T ∇ Λ -1 w d), p * -Λ -1 β • d - 1 2 d T ∇ Λ -1 β d -pD ΓD -q(n • ñ), (β + (∇β)d) • n -hN ΓN (50)
This Euler-Lagrange formulation differs from [START_REF] Mittal | Immersed boundary methods[END_REF], due to the use of the enriched pressure and the mixed formulation to evaluate pressure gradients in the boundary terms. On ΓD , to enforce with higher accuracy the boundary condition, the Taylor expansion is enhanced to third order. As illustrated by the numerical examples of Section 6, such a scheme allows to recover overall second-order accuracy.

Remark 7. The extrapolation of the pressure could be already enhanced to third order without the pressure enrichment. However, such a procedure is unnecessary as the expected accuracy is second order and is already satisfied by a second-order Taylor expansion.

Remark 8. When employed, the additional high-order correction on the tangential gradient ( 45) is also extended using the Taylor expansion:

∂ τ i pD = ∇ pD • τ i = Λ -1 β( x) • τ i + ∇ Λ -1 β d • τ i + O(||d|| 2 ) (51)

Modification of the schemes -Symmetric enrichment

We propose next an improved variant of the ideas discussed so far. For this purpose, let us introduce the mapping R l+1 l : (V k (Ω), S l (Ω)) → S l+1 (Ω), which defines the enriched pressure in the case of linear approximations, (β, p) ∈

V 1 (Ω) × S 1 (Ω): R 2 1 (Ω) : (V 1 (Ω), S 1 (Ω)) → S 2 (Ω) (β, p) → p * (52)

Conformal case

We propose here a symmetrically enriched method in which the operator (52) is used on both the primary variable (the pressure) and on the corresponding test function q in the bilinear form:

B d (w, R 2 1 (w, q); β, R 2 1 (β, q)) = L d (w, R 2 1 (w, q)) (53)
12

With respect to the previous non-symmetric enriched formulation, now the test and trial spaces are identical. Denoting q * = R 2 1 (w, q), we consider the following stabilized variational form:

B d (w, q * ; β, p * ) = B Gal d (w, q * ; β, p * ) + B S tab (w, q * ; β, p * ) = L Gal d (w, q * ) B Gal d (w, q * ; β, p * ) = (w, Λ -1 β) Ω -(∇ • w, p * ) Ω + (q * , ∇ • β) Ω + [[w]], {p * } i -{q * }, [[β]] i + [[q * ]], α β [[p * ]] i + [[w]], α p [[β]] E i + w • n, p * Γ N -q * , β • n Γ N + q * , α β p * Γ D B S tab (w, q * , β, p * ) = 1 2 -Λ -1 w + ∇q * , β + Λ∇p * Ω L Gal d (w, q * ) = (q * , φ) Ω -w • n, p D Γ D + q * , α β p D Γ D -q * , h N Γ N (54)
The proposed symmetric enriched formulation allows for a more straightforward analysis. Setting W = (w, q), W * = (w, R 2 1 (w, q)) and indicating by V e = (β e , p e ) a sufficiently smooth solution, the following theorem provides stability, consistency and accuracy of the formulation (for the proof, see AppendixA):

Theorem 1 (Symmetric enrichment -stability, consistency and error estimate). For all (w, q * ) and (β,

p * ) ∈ V 1 × R 2 1 (V 1 , S 1 (Ω))
, the enriched formulation ( 54) is: a) stable (coercive) with respect to the norm

|||(β, p * )||| 2 = 1 2 (||Λ -1/2 β|| 2 Ω + ||Λ 1/2 ∇p * || 2 Ω ) + ||( αβ ||Λ||/h ⊥ ) 1/2 [[p * ]]|| 2 E i + ||( αβ ||Λ||/h ⊥ ) 1/2 p * || 2 Γ D + ||( αp h ⊥ /||Λ||) 1/2 [[β]]|| 2 (55)
b) consistent (in the sense of the Galerkin orthogonality property), that is

B d (W * , E) = 0 (56)
c) convergent, with the error satisfying the estimate

|||E||| ≤ C 1 h 2 |β| 2 + C 2 h 2 |p| 3 ( 57
)
where E = (β -β e ; p *p e ).

This last error estimates differs from [START_REF] Ern | Theory and practice of finite elements[END_REF], since an extra order of accuracy is gained for the flux, namely r = 2. One additional order of convergence is also expected for the L 2 -error of the pressure. In conclusion the expected rates of convergence for pressure and flux are third and second order, respectively. Remark 9. In the discontinuous case, the B S tab (w, q * , β, p * ) term is not necessary. For a continuous approximation, the "div-div" term was found necessary for the non-symmetric scheme, but unnecessary for the symmetric scheme, for which the momentum stabilization is sufficient to guarantee stability and convergence, as also confirmed in a wide variety of computations. In fact, the symmetric enrichment naturally introduces a stabilization mechanism akin to the "div-div" term:

1 2 -Λ -1 w + ∇q * , β + Λ∇p * Ω = 1 2 -Λ -1 w + ∇q, β + Λ∇p * Ω + 1 2 (∇(q * -q), β + Λ∇p) Ω + 1 2 (∇(q * -q), β + Λ∇(p * -p)) Ω (58)
Focusing on the last term and using (43),

1 2 (∇(q * -q), β + Λ∇(p * -p)) Ω = K∈Ω 1 2 (∇(q * -q), β + Λ∇(p * -p)) K = K∈Ω 1 8        e ∇(∆(Λ -1 w) e • ∆x e ϕ e 1 ϕ e 2 ), Λ e ∇(∆(Λ -1 β) e • ∆x e ϕ e 1 ϕ e 2 )        K ( 59
)
This last term in equality (59) contains the inner product of the gradient of the flux and the gradient of the shape function, similar in nature to the "div-div" stabilization. Note also that the scaling with the square of the grid size (h 2 ) is accounted for through the dot products with ∆x e . Thus, when mapping from the linear unknowns to the quadratic pressure, the "div-div" stabilization is directly included in the momentum stabilization and can be removed from the variational form. This provides a natural and order-preserving stabilization, independent from any stabilization constant.

Embedded case

Proceeding as above, an enhanced extrapolation is applied to the SBM in order to preserve the required order of accuracy. This leads to the following stabilized variational form:

B d (w, q * ; β, p * ) = B Gal d (w, q * ; β, p * ) + B S tab (w, q * ; β, p * ) = L Gal d (w, q * ) + L S tab (w, q * ) B Gal d (w, q * ; β, p * ) = (w, Λ -1 β) Ω -(∇ • w, p * ) Ω + (q * , ∇ • β) Ω + [[w]], {p * } i -{q * }, [[β]] i + [[q * ]], α β [[p * ]] i + [[w]], α p [[β]] E i + q * -Λ -1 w • d - 1 2 d T ∇ Λ -1 w d, α β (p -Λ -1 β • d - 1 2 d T ∇ Λ -1 β d) ΓD + w • ñ, Λ -1 β • d + 1 2 d T ∇ Λ -1 β d ΓD + w • ñ, p * ΓN -q * (n • ñ), (β + [∇β]d) • n ΓN L Gal d (w, q * ) = (q * , φ) Ω -q * (n • ñ), hN ΓN -w • ñ, pD ΓD + q * -Λ -1 w • d - 1 2 d T ∇ Λ -1 w d, α β pD ΓD (60)
Here, as for the original SBM formulation [START_REF] Mazaheri | Efficient high-order discontinuous galerkin schemes with first-order hyperbolic advection-diffusion system approach[END_REF] and for the non-symmetric formulation [START_REF] Wang | Dual-scale galerkin methods for darcy flow[END_REF], we do expect a reduction in accuracy for the pressure when Neumann conditions are embedded. However, the overall accuracy is still secondorder for both flux and pressure. Hence, this enriched scheme and its non-symmetric counterpart deliver second-order accuracy for the flux and pressure, for any kind of embedded boundary conditions. In addition, if only Dirichlet conditions are embedded, the approximated pressure is third-order accurate.

Numerical results

In this section, several numerical examples are presented, to confirm the the theoretical expectations, as well as to provide further insights on the previous discussion on the accuracy and stability of the schemes.

For all two-and three-dimensional simulations, the results are obtained with the symmetric enriched method, using, where necessary, penalty coefficients set to α p = 2 and α β = 1. The non-symmetric formulation is only used in the convergence tests and in comparative studies among the different options for stabilization.

Convergence tests

The convergence tests presented next aim at demonstrating the accuracy of the proposed schemes, and, in the case of the pressure enrichment, tradeoffs between accuracy and overall computational cost. In addition, for Dirichlet boundary conditions, we also illustrate the differences between the formulations and, for the non-symmetric enrichment, the importance of the "div-div" stabilization or the tangential penalty term.

Setup of the tests cases

We use the method of manufactured solutions with matching boundary conditions on a geometry given by the circular annulus of outer radius r out = 0.35 and inner radius r in = 0.1 (see Fig. 3(a)). Note that in one of the tests for the embedded case, the inner circular boundary shape is replaced by a square shape of size l = 0.15 (see Fig. 3(c)), with the goal of demonstrating the robustness of the approach on geometries with singularities. The manufactured solutions employed for this study are defined as

                             β(x, y) =       -e x+y (x 2 + y) -e x+y (x + y 2 )       p(x, y) = 1 3 (x 3 + y 3 ) + xy φ(x, y) = -(y 2 + x 2 + 3x + 3y)e x+y Λ = e x+y I (61)
for the variable permeability case, and as

                                     β(x, y) = -π      
λ xy sin(4πx) cos(4πy) + λ xx cos(4πx) sin(4πy) λ yy sin(4πx) cos(4πy) + λ yx cos(4πx) sin(4πy)

      p(x, y) = sin(4πx) sin(4πy) 4 φ(x, y) = 4π 2 (λ yy + λ xx ) sin(4πx) sin(4πy) -(λ yx + λ xy ) cos(4πx) cos(4πy) Λ =       λ xx λ xy λ yx λ yy       =       3 1 1 3       (62)
for the anisotropic permeability case.

Conformal Results

Dirichlet and Neumann boundary conditions are applied in the outer and inner boundaries, respectively (see Fig. 3(a)). For the variable permeability case (61), convergence rates for both continuous and discontinuous discretizations are given in Figure 4(a) for the flux and Figure 4(b) for the pressure. For the anisotropic permeability case (62), convergence rates for both continuous and discontinuous discretizations are given in Figure 4(c), for the flux, and Figure 4(d), for the pressure. The numerical results confirm the expected convergence rates, in the sense that the original formulation has a second-order accurate pressure and a first-order accurate flux, while the enriched pressure formulation yields a thirdorder pressure and a second-order flux.

These convergence rates are similar to the ones obtained in [START_REF] Thomas | A stabilized mixed discontinuous galerkin method for darcy flow[END_REF] when pressure and velocity are discretized with quadratic and linear interpolation, respectively. However, in the present study, both variables are discretized using piecewise linear interpolation. Therefore, a reduced number of degrees-of-freedom still allows to recover the same accuracy.

Embedded results

We consider two boundary-condition configurations: in a first configuration, termed Dirichlet/Dirichlet, Dirichlet boundary conditions are applied on both the inner and outer boundaries (see Fig. The continuous blue curve corresponds to the symmetric formulation; the continuous red curve corresponds non-symmetric formulation without the "div-div" stabilization or tangential boundary penalty; the dashed red curve corresponds to the non-symmetric formulation with "div-div" stabilization and the dashed magenta curve to the same non-symmetric formulation but with the tangential penalty term instead.

(a) No "div-div" stabilization, no tangential penalty (b) "div-div" stabilization active, no tangential penalty (c) No "div-div" stab, tangential penalty active The reported numerical results match the expected convergence rates. When the pressure enrichment is activated, the flux is always second-order accurate, while the pressure is third-order accurate for pure Dirichlet boundary conditions and second-order accurate for Neumann boundary conditions. We thus validated the gain of one order of accuracy for both the flux and pressure variables using the pressure enrichment.

For reference, the plots of the solutions (with variable permeability distribution on the surrogate domain) are displayed in Figure 7.

The condition number of the algebraic system associated with the proposed methods is plotted, as a function of the mesh size, in Figure 8. One can see that the use of the pressure enrichment has a minimum influence on the condition number.

Comparative study among various weak formulations

This section presents a comparison between the non-symmetric formulation [START_REF] Schott | A new face-oriented stabilized xfem approach for 2d and 3d incompressible navier-stokes equations[END_REF] and the symmetric formulation (54). We focus on the CG case, in which the edge penalty term has no beneficial effect.

In section (5.2), we argued that the non-symmetric modification of the schemes requires the use of "div-div" stabilization or some additional tangential penalty. Here, we investigate the different among the following approaches:

• Non-symmetric enrichment with momentum stabilization and no additional stabilization or tangential penalty term; For r = 1 and the exact solution (64), plots of the pressure and one flux component, with their respective errors, are shown in Figure 12, for both the CG and DG schemes. For r = 1 and the exact solution (65), similar plots are shown in Figure 13, only for the DG scheme. As expected, for the continuous variables, both schemes perfectly match the exact solution. For all these settings, the CG and DG schemes can represent exactly the solutions. However, note that a continuous flux approximation would not work for the exact solution (65), as evident in Figure 14, which shows oscillations near the interface x = 0.5.

Considering now the case r = 2 and the continuous exact solution (64), since CG and DG schemes behave identically, we plot in Figure 15 only the results obtained with the continuous formulation along with the errors obtained with both original and enriched schemes. For the discontinuous solution (65), the simulations are only performed with the discontinuous scheme and the results are displayed in Figure 16. It is clearly seen that while the original linear CG and DG schemes do not allow to recover the exact solutions, the enriched formulations perfectly match quadratic solutions (i.e., within machine precision).

Domain with a circular impermeable obstruction

The domain is Ω = [0, 1] × [0, 1] \ C where C is a circle of radius 0.2 centered in (0.5, 0.5). The permeability in the domain is Λ = I. On the left and right side we apply Dirichlet conditions with respectively p = 1 and p = 0 while all other boundaries are set to slip wall: β • n = 0 (see Fig. 17(a)). As there is no discontinuity in the permeability, the continuous scheme is employed. We compare the solutions for a conformal simulation and one where the inner impermeable obstruction is embedded. The meshes employed are shown in Figures 17(b 19(a)). For the embedded simulation, the entire domain is embedded in a larger square domain. Since the permeability is discontinuous, the DG scheme is used. Meshes are presented in Figures 19(b) and 19(c). We propose a comparison of the solution obtained with the conformal and embedded simulation in Figure 20. Once again, one can appreciate the similarities between the two results, which validate the discontinuous embedded discretization.

A three-dimensional convergence test

We embed a sphere of radius r in = 0.2 in a bigger sphere of radius r out = 0.5. The simulations are performed with 6 levels of refinement, from a mesh of 330 nodes and 1, 420 tetrahedrons to a mesh of 7, 838, 897 nodes and

Conclusion

We have developed enhanced embedded mixed finite elements methods for Darcy-type equations that, for linear approximation spaces, deliver second-order error convergence for the primary and gradient variables. In the case when only pure Dirichlet boundary conditions are embedded, the primary variable becomes third-order accurate.

The constructions performed exploit the mixed form of the problem to develop an enrichment of the primary variable. When both test and trial spaces are enriched, all the classical estimates can be generalized in the body fitted case. Extensive numerical simulations has confirmed all the theoretical expectations concerning both accuracy and stability, also when strong discontinuities in the permeability exist.

Acknowledgements

The first and the last author would like to gratefully acknowledge the support of the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research under the Early Career Research Program Grant SC0012169, U.S. Army Research Office, under Grant W911NF1810308, and ExxonMobil Upstream Research Company (Houston, TX). This work has also greatly benefitted from the support of the Associated Team HAMSTER between Duke University and the CARDAMOM INRIA team, in Bordeaux (France).

AppendixA. Conformal enriched scheme analysis

In this appendix, we provide some analysis of the enriched schemes [START_REF] Schott | A new face-oriented stabilized xfem approach for 2d and 3d incompressible navier-stokes equations[END_REF] and (54). We start by providing the detailed version of Theorem 1 for the symmetric enrichment. We then provide some additional comments on the non symmetric formulation. Also, we only focus on the conformal case, and refer to [START_REF] Main | The shifted boundary method for embedded domain computations. part i: Poisson and stokes problems[END_REF][START_REF] Main | The shifted boundary method for embedded domain computations. part ii: Linear advection-diffusion and incompressible navier-stokes equations[END_REF] for a discussion on the issue of the stability of the SBM approximation.

AppendixA.1. Symmetric enrichment

The main idea in the proof of Theorem 1 resides in the fact that the test and trial function spaces for the pressure are identical. Therefore, the analysis of the scheme can be performed in a similar manner to the original continuous and discontinuous schemes [START_REF] Douglas | Interior penalty procedures for elliptic and parabolic galerkin methods[END_REF] proposed in [START_REF] Thomas | A stabilized mixed discontinuous galerkin method for darcy flow[END_REF][START_REF] Masud | A stabilized mixed finite element method for darcy flow[END_REF].

Proposition 1 (Stability). For all β, p * ∈ V 1 (Ω)×R 2 1 (V 1 (Ω), S 1 (Ω), the enriched formulation (54) is stable (coercive) in the norm:

|||(β, p * )||| 2 = 1 2 ||Λ -1/2 β|| 2 Ω + ||Λ 1/2 ∇p * || 2 Ω + ||( αβ ||Λ||/h ⊥ ) 1/2 [[p * ]]|| 2 E i + ||( αβ ||Λ||/h ⊥ ) 1/2 p * || 2 Γ D + ||( αp h ⊥ /||Λ||) 1/2 [[β]]|| 2 (A.1)
Proof. The proof is immediate, by substitution. For all (w, q

* ) ∈ V 1 (Ω) × R 2 1 (V 1 (Ω) × S 1 (Ω)): B(w, q * ; w, q * ) = |||(w, q * ))||| 2 (A.2)
Setting W = (w, q), W * = (w, R 2 1 (w, q)), we indicate with V e = (β e , p e ) a sufficiently smooth solution. We denote with the symbol • an interpolation operator. Hence, Ṽ * is the interpolate of the reconstruction (52), and ṼP 1 and ṼP 2 are, respectively, the P 1 and P 2 Lagrange interpolants. We then decompose the error as

E = V * -V e = E * V * -Ṽ * e + H * Ṽ * e -V e (A.
(W * , V e ) = L d (W * ): Proposition 2 (Consistency). For all W * ∈ V 1 (Ω) × R 2 1 (V 1 (Ω), S 1 (Ω)): B d (W * , E) = 0 (A.4)
Finally, we can focus on the accuracy of the scheme. Using the decomposition of the error between the approximation and truncation errors, the following theorem holds, for the "classical" non enriched formulation (see [START_REF] Masud | A stabilized mixed finite element method for darcy flow[END_REF] for detailed proof): In our case, the same consistency/stability argument that yields the proof of Theorem 2 can be followed, step-bystep, to obtain |||E||| ≤ C|||H * |||. We then seek an upper bound of H * . To this purpose, we introduce the vector of the P 1 and P 2 interpolates of the flux and pressure Ṽe = ( βP 1 e , pP 2 e ) and we decompose H * as:

H * = Ṽ * e -V e = H * ,2
V * e -Ṽe + , p * -pP 2 ) = (0, p * -p2 ) represents the error between the quadratic reconstruction of the pressure with (52) and the P 2 one, as the mapping (52) approximate linearly the flux. We then have the following Lemma 1. The norm of the error between the two quadratic interpolations H * ,2 is bounded by the interpolation error on the flux β:

|||H * ,2 ||| ≤ C 3 h 2 |β| 2 (A.9)

Proof. Let first recall some useful inequalities. The trace inequality for a simplex K of the mesh writes:

||q|| 2 ∂K ≤ C T (h -1 ||q|| 2 K + h||∇q|| 2 K ) (A.10)
where C T is a constant. In addition, assuming the mesh regular enough, we have .11) and there exists constants C1 and C2 such that, for all edges on the mesh e i→ j , C1 h ≤ ||e i→ j || ≤ C2 h (A.12)

||∇q|| K ≤ C inv h -1 ||q|| K (A
Now, from (A.11), we can provide a stronger form of the trace inequality (A.10) reading:

||q|| 2 ∂K ≤ C T h -1 ||q|| 2 K (A.13)
Then, by definition, H * ,2 = (0, p * -pP 2 ) and thus we have:

|||H * ,2 ||| 2 Ω = 1 2 ||Λ 1/2

AppendixA.2. Non-symmetric enrichment

As explained in Section 5.2, for the non symmetric scheme, the solution (trial) space V 1 (Ω) × R 2 1 (S 1 (Ω), V 1 (Ω)) is not the same as the test space V 1 (Ω) × S 1 (Ω). Note that both S 1 (Ω) and R 2 1 (S l (Ω), V 1 (Ω)) are however sub-spaces of H 1 (Ω), which is the infinite dimensional pressure space for the problem under consideration.

Figure 1 :

 1 Figure 1: Surrogate and true domain/boundary notations.

Figure 2 :

 2 Figure 2: A triangle and the additional degrees-of-freedom for pressure enrichment (indicated with crosses).

  (a) Conformal simulation domain and mesh (b) Embedded simulation with inner embbeded circle, domain and mesh (c) Embedded simulation with inner embbeded square, domain and mesh

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Typical grids and geometries for the convergence tests.

Figure 6 :

 6 Figure 6: Convergence test for the embedded simulations: Neumann/Dirichlet.

Figure 7 :Figure 8 :

 78 Figure 7: Convergence test: Plots of the solution for the isotropic permeability distribution.

  3(b),3(c)); in a second configuration, termed Neumann/Dirichlet, Neumann conditions are applied on the inner boundary while Dirichlet conditions are applied on the outer boundary.For the Dirichlet/Dirichlet simulations, convergence rates obtained with isotropic permeability for both the CG and DG schemes are shown in Figure5(a) and 5(d) for the flux and the pressure variables, respectively. The anisotropic

Figure 9 :

 9 Figure9: Convergence rates: Comparison among various symmetric and non-symmetric formulations. The continuous blue curve corresponds to the symmetric formulation; the continuous red curve corresponds non-symmetric formulation without the "div-div" stabilization or tangential boundary penalty; the dashed red curve corresponds to the non-symmetric formulation with "div-div" stabilization and the dashed magenta curve to the same non-symmetric formulation but with the tangential penalty term instead.

Figure 10 :

 10 Figure 10: Comparison of various non-symmetric formulations. Analysis of spurious modes in the flux.

Figure 11 :Figure 12 :

 1112 Figure 11: Patch tests mesh.

Figure 13 :

 13 Figure 13: Patch test with DG scheme: Elevation plots of the pressure and the flux.

Figure 14 :Figure 15 :

 1415 Figure 14: Flux obtained with CG scheme for the discontinuous exact solution.

Figure 16 :Figure 17 :Figure 18 :

 161718 Figure 16: Patch test with DG scheme. Elevation plots of the quadratic pressure and linear flux.

Figure 21 :

 21 Figure 21: Embedded Neumann 3D simulation: Convergence rates.

Figure 22 :

 22 Figure 22: Monkey trefoil geometry.

Figure 23 :

 23 Figure 23: Monkey trefoil test.

Theorem 2 .

 2 For (β, p) ∈ V l (Ω) × S k (Ω): |||E||| ≤ C|||H|||, (A.5) with |||H||| ≤ C 1 h k+1 |β| k+1 + C 2 h l |p| l+1 (A.6)and |.| s the sth Sobolev seminorm.

H

  Ṽe -V e (A.[START_REF] Burman | Ghost penalty[END_REF] where H = (β P 1 -β e , pP 2p e ) the P 1 and P 2 Lagrange interpolation error of respectively the flux and the pressure, satisfying (A.6)||| H||| ≤ C 1 h 2 |β| 2 + C 2 h 2 |p| 3 (A.8)and H * ,2 = ( β * -βP 1

  ∇( p * -pP 2 )|| 2 Ω + ||( αβ ||Λ/h ⊥ ) 1/2 [[ p * -pP 2 ]]||2 approximate the mid-edge values pP 2 m using the same Taylor expansion as in Section 5.1 and define:β i -Λ i β j ) + H i e i→ j + (Λ j β i -Λ i β j ) + H j e i→ j • e i→ j (A.18)which can be written using the Taylor expansion with Lagrange's reminder (see remark 4):∇ 2 β(c)e i→ j • e i→ j ϕ m (x) (A.20)We can now focus on each term of (A.14).∇ 2 β(c)e i→ j • e i→ j ϕ m (x) ∇( p -pP 2 )|| Ω ≤ C 3 h 4 |β| 2 (A.22)Denoting e p = p * -pP 2 , an identical process can bound the error on edges:||( αβ ||Λ||/h ⊥ ) 1/2 [[e p ]]|| 2 E i +||( αβ ||Λ||/h ⊥ ) 1/2 [[e p ]]|| 2 Γ D ≤ ||( αβ ||Λ||/h ⊥ ) 1/2 [[e p ]]|| E ≤ C 4 h 2 |β| 2 (A.24)Combining (A.22) and (A.24), we have the following final upper bound:|||H * ,2 ||| ≤ (C 3 + C 4 )h 2 |β| 2 (A.25)Thus, the estimate for the interpolation error reads:Proposition 3 (Error estimate). The scheme (54) is accurate with the error estimate:|||E||| ≤ C1 h 2 |β| 2 + C2 h 2 |p| 3 (A.26)Proof. By a consistency/stability argument, we obtain a similar inequality to the case of the scheme without enrichment, namely |||E||| ≤ C|||H * |||. Then, using the triangular inequality on (A.7), as a direct consequence of Theorem 2 and Lemma 1, we finally obtain:|||E||| ≤ C((C 1 + C 4 )h 2 |β| 2 + C 2 h 2 |p| 3 ) (A.27)

	i,node p i + p j p i ϕ P 2 i (Λ pP 2 = + m,edge p i + p j 2 + 1 8 i,node p i ϕ P 2 i + m,edge 2 -1 4 (Λ 3 ) pP 2 m ϕ P 2 m This allows to define the difference in pressure interpolations as: where pP 2 m = p * m -pP 2 m = 1 8 (Λ p * m -pP 2 m = 1 8 e T i→ j ∇ 2 β(c)e i→ j • e i→ j and thus: p * (x) -pP 2 (x) = m,edge 1 8 e T 2 ||Λ 1/2 ∇( p -pP 2 )|| 2 Ω ≤ C inv 2h 2 sup x ||Λ|| K || p -pP 2 || 2 K ≤ C inv 2h 2 sup x ||Λ|| K K p * -pP 2 ≤ C inv 2h 2 sup x ||Λ|| K K m,edge 1 8 e T i→ j 2 ≤ C inv C2 h 4 2 sup x ||Λ|| K K m,edge 1 8        e T i→ j ||e i→ j || ∇ 2 β(c) e i→ j ||e i→ j ||        • e i→ j ||e i→ j || and we can get: 1 2 ||Λ 1/2 αβ h sup x ||Λ|| E i |e p | 2 + + |e p | 2 -+ Γ D |e p | 2 ≤ αβ h sup ∂K ≤ αβ C T h 2 sup x ||Λ|| K K m,edge 1 8        e T i→ j ||e i→ j || ∇ 2 β(c) e i→ j ||e i→ j ||        • e i→ j ϕ m (x) ||e i→ j || ϕ m (x) 2 i→ j 1 and thus:	2	(A.16) (A.17) (A.19) (A.21) (A.23)

E

(A.

[START_REF] Douglas | Self-adaptive finite element simulation of miscible displacement in porous media[END_REF]

where E = E i Γ D . For every simplex K we have, by definition (

42

)

p * = i β i -Λ j β j ) • e i→ j ϕ P 2 m (A.15) and we can i β i -Λ j β j ) + 1 16 e i→ j (H i + H j )e i→ j + O(h j x ||Λ|| K ||pp * || 2 x ||Λ|| K ||pp * || 2 K ≤ αβ C T C2 h 4 sup

Shifted boundary methods in mixed form4.1. General notations, the true and surrogate domainsLet us consider now an embedded discretization, in which the computational grid does not conform to the boundary (see Figure1). We introduce the notion of a surrogate boundary Γ, composed of the faces/edges of the mesh that are

• Non-symmetric enrichment with momentum and "div-div" stabilization; no tangential penalty term;

• Non-symmetric enrichment with momentum stabilization and tangential penalty term (with α τ = 2.0);

• Symmetric enrichment with momentum stabilization.

The convergence plots obtained with Dirichlet/Dirichlet combination are displayed in Figure 9 show that the nonsymmetric enrichment with only the momentum stabilization is unstable (see also the spurious modes appearing near the embedded boundary in Figure 10). The addition of either the "div-div" stabilization or the tangential penalty term reestablishes stability in the solution with the expected convergence rate, similar to the one obtained with the symmetric enriched formulation.

Patch test cases

We now tackle the exact representation of polynomial solutions of the same degree of the finite element approximation spaces utilized. For this purpose, we perform tests similar to the ones presented in [START_REF] Nishikawa | On hyperbolic method for diffusion with discontinuous coefficients[END_REF] with discontinuous permeabilities. We consider again an annulus geometry with inner radius r in = 0.1 and outer radius r out = 0.5, embedded in a disk of radius R = 0.6. Neumann and Dirichlet conditions are imposed on the inner and outer boundaries, respectively. The computational grid is composed of 2, 531 nodes and 4, 908 triangles, and the interface x = 0.5 is meshed explicitly (i.e., as body fitted, see Fig. 11). The permeability distribution is chosen to be Λ = κ(x, y)I, where

We consider a generalization of the exact solutions considered in [START_REF] Nishikawa | On hyperbolic method for diffusion with discontinuous coefficients[END_REF] to any polynomial degree r:

with κ = 0.5(κ 1 + κ 2 ) + 4κ 1 κ 2 , and

Note that in the last two equations of (64)-(65) we define x α = 0 if α < 0. Observe also that (64) implies that both the pressure and the flux (gradient) are continuous, while in (65) the flux is discontinuous. Hence, a patch test with (64) can be in principle passed by both DG and CG schemes with pressure enrichment, while (65) can be solved exactly only by the DG scheme. We perform the simulations for a linear pressure / constant flux (r = 1) and quadratic pressure / linear flux (r = 2) with and without the enrichment. 46, 530, 560 tetrahedrons. The exact solution is defined as:

cos(y) sinh(z) + sin(x) cosh(y) + 2y cosh(z)(sin(y) + sin(x)) -2z

The convergence rates obtained with and without the pressure enrichment are plotted in Figure 21. There we can appreciate an analogous behavior with respect to the two-dimensional cases, that is, second-order convergence for the flux and pressure in the case of pressure enrichment, and first-order convergence for both fields if no enrichment is provided.

Flow past a complicated three-dimensional object

In order to demonstrate the applicability of the proposed algorithms, we simulate the flow around a complicated object, named "Monkey Trefoil" and already employed in [START_REF] Main | The shifted boundary method for embedded domain computations. part i: Poisson and stokes problems[END_REF] (see Fig. 22). The surface is embedded in a domain of Thus, the stability cannot be proven with same coercivity argument as for the symmetric formulation, but an argument can be provided regarding the consistency and error of consistency.

Consistency can be proved in a very similar manner to case of the symmetric formulation. If (β e , p e ) are smooth fields verifying exactly the PDE and the boundary conditions, one can readily show that 0 = B d (w, q; ββ e , R 2 1 (β, p)p e ) = B d (w, q; ββ e , p *p e ) Error analysis. As before, the global error can be split into the approximation and the truncation error. Note that the consistency/orthogonality relation still implies

And we focus on the consistency error |B d (W, E * )|, as shown in [START_REF] Abgrall | High order methods for CFD[END_REF][START_REF] Abgrall | High order fluctuation schemes on triangular meshes[END_REF]. Using the continuity of B d provides an estimate of the consistency error, namely,

where the norms on the right hand side are the natural norms in H div × H 1 , and we have:

Here p represents the interpolation with or without the enrichment. We then have, using inverse inequality (A.11),

Here the interpolation error can be used with p = pP 1 for the original formulation We can note that although the final estimate is O(h), due to the divergence part in the classical H div norm, the consistency error estimates coming from the pressure are improved when switching from the non enriched to the enriched scheme.