domain: math spi

text:
Pablo G S Carvalho
email: pablogscarvalho@gmail.com

Philippe R B Devloo

Sônia M Gomes
email: soniag@ime.unicamp.br

A hybrid-mixed method for Stokes-Brinkman-Darcy flows with H(div)-velocity fields

Keywords: Hybrid-mixed method, Stokes, Brinkman and Darcy's flows, Divergence-free simulation, Static condensation

We consider hybrid-mixed finite element formulations for Stokes-Brinkman problems. Using H(div)-conforming approximate velocity fields, the continuity of normal components over element interfaces is taken for granted, and pressure is searched in discontinuous spaces preserving the divergence compatibility property. Tangential continuity is weakly imposed by a traction Lagrange multiplier. The method is strongly mass-conservative, leading to exact divergence-free simulations of incompressible flows. The multiplier space requires specific choices according to the velocity approximations implemented in each element geometry. In certain cases, classic divergence-compatible pairs adopted for Darcy's flows may require divergence-free bubble enrichment to enforce tangential continuity in some extent, avoiding any extra stabilization technique. An error analysis typically used for non-conforming methods reveals estimates in terms of optimal errors and consistency errors. Considerable improvement in computational performance is achieved by the application of static condensation: the global system is solved only for a piecewise constant pressure variable, velocity normal trace and tangential traction over interfaces. The remaining solution components are recovered by solving independent local Neumann problems in each element. Numerical results are presented for verification of the main convergence properties of the method in the whole range of parameters, from Stokes to Darcy limits, as well as for the combined Stokes-Darcy scenario.

Introduction

Fluid mechanics modeling has been the subject of different computational and experimental investigations, with a variety of applications in science and engineering. Our focus is on a finite element (FE) numerical approach for the fundamental Stokes, Brinkman and Darcy's models for incompressible fluid flows. For instance, they are usually applied in modeling flow through rough-walled rocks, fractures and small cavities in porous media, called "vugs", and also appear in surface-water hydrology, filtration and lubrication problems [START_REF] Slattery | Advanced transport phenomena[END_REF].

Stokes and Brinkman problems consider fluid flows occupying a region Ω, with constant viscosity 𝜇 > 0, with Dirichlet boundary condition 𝑢 𝐷 ∈ 𝐻 1/2 (𝜕Ω, R 𝑑 ) enforced in 𝜕Ω. We assume that ∫ 𝜕Ω 𝑢 𝐷 • 𝑛 𝑑𝑠 = 0, and the constitutive equation 𝜎 = 𝜎(𝑢, 𝑝) = 2𝜇 𝐷(𝑢) -𝑝𝐼 for the flow stress tensor, where 𝐷(𝑢) = 1 2 (∇𝑢 + ∇𝑢 𝑇 ) is the symmetric gradient tensor. Both cases can be formulated in a unified form, in terms of the velocity 𝑢 and pressure 𝑝 fields, by the system of momentum balance and conservation of mass equations

-∇ • 𝜎(𝑢, 𝑝) + 𝛼𝑢 = 𝑓 in Ω, (1)
∇ • 𝑢 = 𝑔 in Ω, ( 2
)
𝑢 = 𝑢 𝐷 in 𝜕Ω, (3)
where 𝑓 ∈ 𝐿 2 (Ω, R 𝑑 ) and 𝑔 ∈ 𝐿 2 (Ω) are force terms. The Brinkman problem occurs for 𝛼 > 0 defining a dynamic viscosity divided by the permeability. The Stokes problem is the limit case when 𝛼 = 0 together with 𝑔 = 0. Observe that the model can also be interpreted as a Darcy's problem with Neumann boundary condition in the limit case 𝜇 = 0 with 𝑓 ≡ 0. Our goal is to propose, analyze and implement a new hybrid-mixed nonconforming formulation for this problem using divergence-compatible FE pairs in H(div, Ω) × L 2 (Ω) for velocity and pressure discretizations.

There exist in the literature several FE discretizations for this mathematical model [START_REF] Girault | Finite element methods for Navier-Stokes equations: Theory and algorithms[END_REF]. For the Stokes problem, a usual approach is based on a mixed variational formulation, having the velocity of the fluid displacement as primary variable, with continuous traces over element interfaces. It can be characterized as the solution of a constrained minimization problem, where discontinuous pressure fields play the role of a Lagrange multiplier to impose the incompressibility constraint. The main difficulty faced by Stokes 𝐻 1 -conforming FE models, using continuous velocity fields, occurs in the treatment of this divergence-free condition. For well-posedness, velocity and pressure approximation spaces should verify a rather strict compatibility property, the well known BBL or inf-sup condition [START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from lagragian multipliers[END_REF][START_REF] Boffi | Mixed and hybrid finite element methods[END_REF]. For example, using equal order Lagrangian elements for both velocity and pressure results in unstable approximations.

To relax this constraint, a variety of nonconforming velocity-pressure mixed methods have been proposed, violating the inter-element continuity condition of the velocity fields to some extent. Let us recall some lines of investigations on this topic, remarking that giving a detailed account of related literature lies out of the scope of the present work.

• In the classic seminal work by Crouzeix-Raviart [START_REF] Crouzeix | Conforming and non-conforming finite element methods for solving the stationary Stokes equations[END_REF], the principle is to take velocities that are only continuous at appropriate Gauss points on triangle edges. A general theory for analyzing nonconforming methods applied to Stokes problems can also be found in this publication. We refer to [START_REF] Brenner | Forty years of the Crouzeix-Raviart element[END_REF] for a survey review on this matter.

• Discontinuous Galerkin (DG) methodology discretizes the variables in full broken spaces, without assuming any continuity constraint in the approximation spaces [START_REF] Toselli | hp discontinuous Galerkin approximations for the Stokes problem[END_REF][START_REF] Girault | A discontinuous galerkin method with nonoverlapping domain decomposition for the Stokes and Navier-Stokes problems[END_REF][START_REF] Lazarov | Stabilized discontinuous finite element approximations for Stokes equations[END_REF]. As requested in any DG FE method, extra stabilization terms are incorporated to the variational formulation in order to have stable and accurate simulations.

• Nonconforming approaches in between 𝐻 1 -conforming FE methods and full DG methods use divergence-compatible FE pairs in H(div, Ω) × L 2 (Ω) for velocity and pressure discretizations. Namely, the approximate velocity fields are assumed from start to have continuous normal traces across element interfaces. In this direction, there is the semi-DG proposal in [START_REF] Wang | A robust numerical method for Stokes equations based on divergence-free H(div) finite element methods[END_REF], where the penalization is introduced just on the velocity tangential components over element interfaces. We refer to [START_REF] Carvalho | On the use of divergence balanced H(div)-L2 pair of approximation spaces for divergence-free and robust simulations of Stokes, coupled Stokes-Darcy and Brinkman problems[END_REF] for a unified representation and a comparative study of some known full DG and semi-DG schemes for Stokes, as well for Brinkmann and coupled Stokes-Darcy, problems.

• Penalization can be avoided when divergence-compatible FE pairs in H(div, Ω)×L 2 (Ω) are adopted, but the enforcement of some kind of tangential continuity is needed.

For that, the velocity FE space is enriched with some properly chosen divergencefree bubble velocity FE of higher degree. The original pressure space is maintained, without spoiling the divergence consistency property. In this direction, we refer to [START_REF] Mardal | A robust finite element method for Darcy-Stokes flow[END_REF][START_REF] Tai | A discrete de rham complex with enhanced smoothness[END_REF][START_REF] Xie | Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models[END_REF][START_REF] Guzmán | A family of nonconforming elements for the Brinkman problem[END_REF] for methods based on simplexes and to [START_REF] Chen | Uniformly convergent H(div)-conforming rectangular elements for Darcy-Stokes problem[END_REF][START_REF] Chen | Uniformly convergent cubic nonconforming element for Darcy-Stokes problem[END_REF] for quadrilateral and hexahedral meshes. Excepting [START_REF] Guzmán | A family of nonconforming elements for the Brinkman problem[END_REF], all these methods are of lowest order.

Another class of nonconforming methods occurs by hybridization: by relaxing the velocity continuity between elements of a domain partition, an appropriate Lagrange multiplier is introduced to weakly impose it. The term hybrid comes from this additional variable used for the approximation of fields over the skeleton of the partition. In the context of structural mechanics, the importance of hybridization was recognized since the beginning of FE history, and originally presented in [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF]. Since then, these principles have been widely explored in the design of new formulations in many fields, motivated by the growing interest in solving high precision problems with optimized computational processes. For them, model resolution is divided into local independent problems, favoring parallel computing implementations, and profiting from static condensation techniques (e.g., see [START_REF] Arnold | Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates[END_REF]). Moreover, when the multiplier field admits a physical interpretation, giving further information about the solution, it can be exploited in a postprocess to improve accuracy and/or for posteriori error analysis .

For the hybrid-mixed formulation we have in mind, using divergence-compatible FE pairs in H(div, Ω) × L 2 (Ω) for velocity and pressure discretizations, the velocity normal trace continuity is taken for granted, and requirement for continuous tangential velocity components is weakened but not removed completely from the approximation spaces (as it is done in full DG methods). Given a geometric partition 𝒯 = {𝐾} of the flow domain, a new traction variable (multiplier) is introduced over the facets of the mesh skeleton ℰ (edges or faces of element boundaries). Well-posedness of the method holds with proper choices of the FE space for the traction variable, avoiding the modification of the formulation with extra stabilization treatment. An error analysis reveals typical estimates held by nonconforming methods for these problems.

The proposed hybrid-mixed formulation offers some additional attractive characteristics. As listed in the sequel, some of them correspond to the characteristics of other hybrid meth-ods, specially with respect to the MHM multiscale hybrid mixed method (MHM) recently proposed in [START_REF] Araya | Multiscale hybrid-mixed method for the Stokes and Brinkman equations -the method[END_REF] for the Stokes-Brinkman problem. However, the use divergence-compatible FE pairs in H(div, Ω) × L 2 (Ω) for velocity and pressure discretizations offers some additional positive properties.

1. The method is strongly locally conservative, and naturally gives exact divergence-free velocity fields and pressure robustness with respect to the limiting case of Darcy flows, properties that few schemes can achieve, as discussed in the survey [START_REF] John | On the divergence constraint in mixed finite element methods for incompressible flows[END_REF]. 2. A static-condensation scheme can be applied, by taking the traction multiplier, the velocity normal trace over element interfaces, and the piecewise constant element average pressure as primary variables, given by the unique solution of a global system of equations. 3. The remaining velocity and pressure degrees of freedom are secondary variables. They are recovered by solving local Neumann problems restricted to each element 𝐾. 4. The method allows much more efficient parallel implementations since the degrees of freedom on the element boundary are located at the facets (edges or faces) such that the local communication of common data affects at most two elements. 5. The coupling of Stokes, Brinkman and porous media flows can be implemented naturally.

Other hybrid formulations for Stokes and/or Brinkman problems have been proposed in the literature from different approaches, but also following the "divide-and-conquer" principle. The hybrid DG method introduced and analyzed in [START_REF] Egger | hp analysis of a hybrid DG method for Stokes flow[END_REF] adopts Lagrange multipliers associated with the trace of the velocity field. Another variant of this method is considered in [START_REF] Igreja | Stabilized velocity and pressure mixed hybrid DGFEM for the Stokes problem[END_REF]. There are also the so called Hybrid Discontinuous Galerkin methods (HDG). In addition to velocity and pressure, the technique has been developed by the incorporation of vorticity, velocity gradient or stress variables. The local variables are piecewisely defined by polynomials of the same degree for all the components of the approximate solution, and solved globally for the coupled condensed systems for the approximate trace of the velocity over interfaces and the mean of the pressure on the elements. The local solvers approximate the remaining field components in terms of the numerical trace of the velocity and the mean of pressure. Divergence-free and 𝐻(div)-conformity properties require a post-processed velocity. We refer the reader to the survey in [START_REF] Cockburn | A comparison of HDG methods for Stokes flow[END_REF] for the specific case of Stokes problems. The method considered in a more recent context, called Hybrid High-Order (HHO) [START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF], is formulated in terms of face-based and internal degrees of freedom in each element. The variational formulation takes into account a velocity reconstruction operator on a polynomial space of higher degree, and requires the adoption of reconstruction operators. For the sake of well posedness, all the aforementioned methods require the incorporation in the equations of extra stabilization terms.

Outline of the paper

After some preliminary comments concerning notation and the adopted functional framework, the model problem and some of its classic mixed and hybrid-mixed variational formulations are presented in Section 2. We introduce a new hybrid-mixed variant using H(div)conforming velocity field in Section 2.3 and discuss FE discretizations for it in Section 3. They are based on divergence-compatible FE pairs for approximate velocity and pressure fields. We describe the additional trace compatibility conditions on the FE spaces required for solution existence and uniqueness. We also discuss convergence issues for the method. In the sequel, the focus in Section 4 is on some specific examples of FE space configurations using classic divergence-compatible FE pairs based on triangular, quadrilateral, tetrahedral or hexahedral meshes, for which the trace compatibility property is documented. We propose a static-condensation procedure for the hybrid-mixed finite element method in Section 5. Simulation results for some test problems integrate Section 6 for the verification of predicted convergence rates, velocity divergence-free and pressure robustness properties.

Weak formulations for the model problem

We introduce some basic notation and recall some variational formulations for the model problem: the classic mixed method and the mixed method of current interest, combined with hybridization.

Comments about notation

Let Ω ⊂ R 𝑑 , 𝑑 ∈ {2, 3}, be a connected polyhedral domain. Variational formulations shall be presented here as subordinated to a partition 𝒯 = {𝐾} of Ω, 𝐾 having one of the usual element geometry. The set Γ formed by all element facets 𝐸 (edges in 2D or faces in 3D) is called the mesh skeleton, and Γ o = {𝐸 ∈ Γ : 𝐸 ⊂ Ω} denotes the set of internal facets. Over Γ, once and for all, a vector field 𝑛 is fixed such that 𝑛| 𝐸 is a unit normal associated to the facet 𝐸 ∈ Γ.

Throughout the text, for any subregion 𝐷 ⊆ Ω, denote by 𝑛 𝐷 the external unitary normal to 𝜕𝐷. The scalar Hilbert spaces 𝐿 2 (𝐷) and 𝐻 𝑠 (𝐷) and norms || • || 𝐿 2 (𝐷) and || • || 𝐻 𝑠 (𝐷) have the usual meaning. The 𝐿 2 -inner product is denoted by (., .) 𝐷 , and 𝐿 2 0 (𝐷) stands for the space of functions 𝜓 ∈ 𝐿 2 (𝐷) with zero mean value, i.e., (𝜓, 1) 𝐷 = 0. Associated spaces 𝐿 2 (𝐷, Z) and 𝐻 𝑠 (𝐷, Z) shall also be considered, taking values in either Z = R 𝑑 or in matrix spaces Z = M = R 𝑑×𝑑 . These spaces inherit the corresponding norms associated to the inner products in 𝐿 2 (𝐷) and 𝐻 𝑠 (𝐷). The space 𝐻(div, 𝐷) denotes the square-integrable vector functions, taking values in R 𝑑 , for which the divergence is also square integrable, and 𝐻 0 (div, 𝐷) = {𝑣 ∈ 𝐻(div, 𝐷); 𝑣 • 𝑛 𝐷 = 0}. Similarly, we shall consider tensor functions in 𝐻(div, 𝐷, M), recalling that for 𝜎 ∈ 𝐻(div, 𝐷, M) the divergence ∇ • 𝜎 is the vector field obtained by taking the divergence of each row. These spaces inherit the corresponding norms associated to the inner products (., .) 𝐷 in 𝐿 2 (𝐷, R 𝑑 ), and 𝐿 2 (𝐷, M). The notation ⟨ , ⟩ 𝜕𝐷 defines the duality pairing between the two function spaces

𝐻 1/2 (𝜕𝐷, R 𝑑 ) = { 𝜈 = 𝑢| 𝜕𝐷 , 𝑢 ∈ 𝐻 1 (𝐷, R 𝑑 ) } , 𝐻 -1/2 (𝜕𝐷, R 𝑑 ) = { 𝜂 = 𝜂 𝑛 𝐷 , 𝜂 ∈ 𝐻(div, 𝐷, M) } .
The subscript 𝐷 will be dropped if 𝐷 = Ω. We also denote ⟨𝜂, 𝜈 ⟩ 𝜕𝒯 = ∑ 𝐾∈𝒯 ⟨𝜂, 𝜈⟩ 𝜕𝐾 . We set Ψ = 𝐿 2 (Ω) = Ψ ⊕ Ψ ⊥ , where Ψ ⊂ Ψ denotes the piecewise constant functions based on 𝒯 , and Ψ ⊥ is its 𝐿 2 -orthogonal complement in Ψ (i.e., of piecewise zero-meam functions over 𝐾 ∈ 𝒯 ). Moreover, consider Ψ 0 = 𝐿 2 0 (Ω). To relax the constraint of interelement continuity, we introduce the spaces

V = 𝐻 1 (𝒯 , R 𝑑 ) = { 𝑣 ∈ 𝐿 2 (Ω, R 𝑑 ) : 𝑣| 𝐾 ∈ 𝐻 1 (𝐾, R 𝑑 ), ∀𝐾 ∈ 𝒯 } 𝑋 = V ∩ 𝐻(div, Ω).
A direct sum decomposition 𝑋 = 𝑋 𝜕 + X holds in terms of bubble velocity functions X = {v ∈ 𝑋; v • 𝑛| 𝜕𝐾 = 0, ∀𝐾 ∈ 𝒯 }, whilst 𝑋 𝜕 denotes a complementary space of velocity fields with some sort of non-vanishing normal trace on element boundaries.

Hybridization requires stress trace spaces defined over the mesh skeleton

Λ = { 𝜎𝑛| 𝜕𝐾 ∈ 𝐻 -1/2 (𝜕𝐾, R 𝑑 ), ∀𝐾 ∈ 𝒯 : 𝜎 ∈ 𝐻(div, Ω, M) } .
A direct decomposition Λ = Λ 𝑛 + Λ 𝑡 holds in terms of normal and tangential components

Λ 𝑛 = {𝜆 𝑛 = (𝜆 • 𝑛)𝑛, 𝜆 ∈ Λ} , Λ 𝑡 = { 𝜆 𝑡 = 𝜆 -(𝜆 • 𝑛)𝑛, 𝜆 ∈ Λ } .
The spaces V (and and also 𝑋) are provided with the broken 𝐻

(Ω, R 𝑑 ) = {𝑣 ∈ 𝑋; ⟨𝜂 𝑡 , 𝑣⟩ 𝜕𝒯 = 0, ∀𝜂 𝑡 ∈ Λ 𝑡 }.

Mixed and hybrid-mixed formulations for Stokes and Brinkman problems

In this section we summarize some mixed and hybrid formulations for the Stokes-Brinkman problem (1)-(3). The classic weak formulation is expressed in the mixed form: find 𝑢 ∈ 𝐻 1 (Ω, R 𝑑 ), 𝑢 = 𝑢 𝐷 in 𝜕Ω, and 𝑝 ∈ Ψ 0 , such that

𝑎(𝑢, 𝑣) -𝑏(𝑣, 𝑝) = (𝑓 , 𝑣), ∀𝑣 ∈ 𝐻 1 0 (Ω, R 𝑑 ), (4)
𝑏(𝑢, 𝑞) = (𝑔, 𝑞), ∀𝑞 ∈ Ψ 0 , ( 5
)
where the bilinear forms are 𝑎(𝑤, 𝑣) = (2𝜇𝐷(𝑤), 𝐷(𝑣)) + 𝛼(𝑤, 𝑣) and 𝑏(𝑣, 𝑝) = (∇ • 𝑣, 𝑝). We observe from the assumption 𝑓 ∈ 𝐿 2 (Ω, R 𝑑 ) that 𝜎 ∈ 𝐻(div, M). Notice that in Darcy's limit (𝜇 = 0 and 𝑓 ≡ 0) 𝐻 1 (Ω, R 𝑑 ) is not the appropriate space for 𝑢, and it should be replaced by 𝐻(div, Ω).

There is another mixed formulation counterpart, which is combined with hybridization. It is based on the broken space setting V × Ψ × Λ, and use the bilinear forms

𝑎 𝒯 (𝑢, 𝑣) = ∑ 𝐾∈𝒯 (2𝜇𝐷(𝑢), 𝐷(𝑣)) 𝐾 + (𝛼𝑤, 𝑣), 𝑏 𝒯 (𝑣, 𝑝) = ∑ 𝐾∈𝒯 (∇ • 𝑣, 𝑝) 𝐾 .
This hybrid-mixed formulation searches for (𝑢, 𝑝, 𝜆, 𝜌)

∈ V × Ψ × Λ × R such that 𝑎 𝒯 (𝑢, 𝑣) -𝑏 𝒯 (𝑣, 𝑝) + ⟨𝜆, 𝑣⟩ 𝜕𝒯 = (𝑓 , 𝑣), ∀𝑣 ∈ V, (6)
𝑏 𝒯 (𝑢, 𝑞) + (𝜌, 𝑞) = (𝑔, 𝑞), ∀𝑞 ∈ Ψ, (7)
⟨𝜂, 𝑢⟩ 𝜕𝒯 = ⟨𝜂, 𝑢 𝐷 ⟩, ∀𝜂 ∈ Λ, (8)
(𝜉, 𝑝) = 0, ∀𝜉 ∈ R. (9)
The new variable 𝜆 plays the role of Lagrange multiplier to weakly enforce the velocity continuity over Γ o and the boundary condition, whilst 𝜌 is introduced to enforce the pressure zero-mean constraint [START_REF] Lazarov | Stabilized discontinuous finite element approximations for Stokes equations[END_REF]. As proved in [START_REF] Araya | Multiscale hybrid-mixed method for the Stokes and Brinkman equations -the method[END_REF], these two mixed formulations are equivalent in the following sense.

Theorem 1. The function (𝑢, 𝑝) ∈ 𝐻 1 (Ω, R 𝑑 )×𝐿 2 0 (Ω) is the unique solution of the standard weak formulation (4)-( 5) if and only if (𝑢, 𝑝, 𝜆, 0) ∈ V × Ψ × Λ × R is the unique solution of the hybrid-mixed formulation (6)- [START_REF] Lazarov | Stabilized discontinuous finite element approximations for Stokes equations[END_REF], for

𝜆| 𝐸 = -𝜎𝑛| 𝐸 , ∀𝐸 ∈ Γ,
where 𝜎 = 𝜎(𝑢, 𝑝) = 2𝜇 𝐷(𝑢) -𝑝.

Hybrid-mixed formulation for velocity fields in H(div, Ω)

Our focus is on a variant of the hybrid-mixed formulation (6)- [START_REF] Lazarov | Stabilized discontinuous finite element approximations for Stokes equations[END_REF]. The principle is to assume, from start, velocity fields with continuous normal traces across element interfaces, i.e. 𝑢 ∈ 𝑋. In this sense, the continuity requirement of the classic 𝐻 1 -conforming mixed formulation is weakened in this nonconforming method, but it is stronger than in fully discontinuous methods, as for the hybrid-mixed formulation (u ∈ 𝐻 1 (𝒯 , R 𝑑 )). As shall be verified, the three approaches give the same solution. Differences appear in their FE discretizations.

For that, firstly observe that for 𝑣 ∈ 𝑋, 𝑏 𝒯 (𝑣, 𝑝) = 𝑏(𝑣, 𝑝). Moreover, let the traction variable 𝜆 ∈ Λ be expressed as 𝜆 = 𝜆 𝑛 +𝜆 𝑡 in terms of its normal and tangential components:

𝜆 𝑛 = (𝜆 • 𝑛)𝑛, 𝜆 𝑡 = 𝜆 -𝜆 𝑛 ∈ Λ 𝑡 , Thus, ⟨𝜆, 𝑣⟩ 𝜕𝒯 = ⟨𝜆 𝑛 , 𝑣⟩ 𝜕𝒯 + ⟨𝜆 𝑡 , 𝑣⟩ 𝜕𝒯 = ⟨𝜆 𝑡 , 𝑣⟩ 𝜕𝒯 holds for 𝑣 ∈ 𝑋.
Consequently, by setting 𝒮 = 𝑋 × Ψ × Λ 𝑡 , we may consider the alternative hybrid-mixed formulation (denoted by HM-H(div)): find (𝑢, 𝑝, 𝜆 𝑡 , 𝜌) ∈ 𝒮 × R such that 𝑎 𝒯 (𝑢, 𝑣) -𝑏(𝑣, 𝑝) + ⟨𝜆 𝑡 , 𝑣⟩ 𝜕𝒯 = (𝑓 , 𝑣), ∀𝑣 ∈ 𝑋, [START_REF] Wang | A robust numerical method for Stokes equations based on divergence-free H(div) finite element methods[END_REF])

𝑏(𝑢, 𝑞) + (𝜌, 𝑞) = (𝑔, 𝑞), ∀𝑞 ∈ Ψ, ( 11
)
⟨𝜂 𝑡 , 𝑢⟩ 𝜕𝒯 = ⟨𝜂 𝑡 , 𝑢 𝐷 ⟩, ∀𝜂 𝑡 ∈ Λ 𝑡 , ( 12
) (𝜉, 𝑝) = 0, ∀𝜉 ∈ R. ( 13
)
Theorem 2. The function (u, 𝑝, 𝜆, 0) ∈ V × Ψ × Λ × R is the unique solution of the hybrid-mixed formulation (6)-( 9) if and only if (𝑢, 𝑝, 𝜆 𝑡 , 0) ∈ 𝒮 × R solves the HM -H(div) formulation (10)-( 13), 𝜆 𝑡 being the tangential component of 𝜆.

Proof. Suppose (𝑢, 𝑝, 𝜆, 𝜌) ∈ 𝒮 × R is the unique solution of ( 6)-( 9). For comparison, we adopt the momentaneous notation ( û, p, λ𝑡 , ρ) ∈ 𝒮 div × R for a solution of the hybridmixed-H(div) formulation ( 10)- [START_REF] Tai | A discrete de rham complex with enhanced smoothness[END_REF]. Expressing 𝜂 = 𝜂 𝑡 + 𝜂 𝑛 ∈ Λ, equations ( 8) and ( 12) imply that ⟨𝜂, u -û⟩ 𝜕𝒯 = ⟨𝜂, u⟩ 𝜕𝒯 -⟨𝜂 𝑡 , û⟩ 𝜕𝒯 = 0, ∀𝜆 ∈ Λ. Thus, by Lemma 4 in [START_REF] Araya | Multiscale hybrid-mixed method for the Stokes and Brinkman equations -the method[END_REF],

u -û ∈ 𝐻 1 0 (Ω, R 𝑑 ). That is, û ∈ 𝐻 1 (Ω, R 𝑑
) and û| 𝜕Ω = 𝑢 𝐷 . Taking any 𝑞 ∈ Ψ ⊥ in equations ( 7) and ( 11), we obtain that 𝑏 𝒯 (𝑢 -û, 𝑞) + (𝜌 -ρ, 𝑞) = 𝑏(𝑢 -û, 𝑞) = 0. Moreover, equations ( 9) and ( 13) mean that 𝑝 -p ∈ Ψ ⊥ . Consequently, choosing any 𝑣 ∈ 𝐻 1 0 (Ω, R 𝑑 ) in equations ( 6) and ( 10) we obtain 𝑎(𝑢 -û, 𝑣) + 𝑏(𝑣, 𝑝 -p) = 0. Consequently 𝑢 -û and 𝑝 -p give the unique classic solution of the weak formulation ( 4)-( 5) with homogeneous data 𝑓 = 0 and u 𝐷 = 0. That is, 𝑢 = û and 𝑝 = p. Moreover, for constant test function 𝑞 = ρ, equation [START_REF] Carvalho | On the use of divergence balanced H(div)-L2 pair of approximation spaces for divergence-free and robust simulations of Stokes, coupled Stokes-Darcy and Brinkman problems[END_REF] becomes ||ρ|| 2 = -𝑏( û, ρ) = (∇ • û, ρ) = ⟨𝑢 𝐷 • 𝑛 Ω , ρ⟩ = 0, implying ρ = 0. Finally, the difference of equations ( 6) and ( 10) reduces to ⟨𝜆 -λ𝑡 , 𝑣⟩ 𝜕𝒯 = ⟨𝜆 𝑡 -λ𝑡 , 𝑣⟩ 𝜕𝒯 = 0, for arbitrary test functions 𝑣 ∈ 𝑋, and we conclude that λ𝑡 = 𝜆 𝑡 . Remark 1. In the above proof we used the property stating that if ⟨𝜆 𝑡 , 𝑣⟩ 𝜕𝒯 = 0, for arbitrary 𝑣 ∈ 𝑋, then 𝜆 𝑡 = 0 must hold. That is, ⟨𝜆 𝑡 , 𝑣⟩ 𝜕𝒯 defines an injective linear operator in 𝑋 for any 𝜆 𝑡 ∈ Λ 𝑡 .

Hybrid-mixed FE method for velocity fields in H(div, Ω)

The purpose is to construct discrete versions for the hybrid-mixed formulation ( 10)-( 13) using FE space settings 𝒮 ℎ = 𝑋 ℎ × Ψ ℎ × Λ 𝑡 ℎ ⊂ 𝒮 div , indexed by the characteristic mesh width ℎ of the associated partitions 𝒯 ℎ (and of the mesh skeleton Γ ℎ ), which are assumed to be shape-regular, with shape-regularity factors independent of ℎ.

HM-H(div)(𝒮 ℎ ): find (𝑢 ℎ , 𝑝 ℎ , 𝜆 𝑡 ℎ , 𝜌 ℎ ) ∈ 𝒮 ℎ × R such that 𝑎 𝒯 ℎ (𝑢 ℎ , 𝑣 ℎ ) -𝑏(𝑣 ℎ , 𝑝 ℎ ) + ⟨𝜆 𝑡 ℎ , 𝑣 ℎ ⟩ 𝜕𝒯 ℎ = (𝑓 , 𝑣 ℎ ), ∀𝑣 ℎ ∈ 𝑋 ℎ , ( 14
)
𝑏(𝑢 ℎ , 𝑞 ℎ ) + (𝜌 ℎ , 𝑞 ℎ ) = (𝑔, 𝑞 ℎ ), ∀𝑞 ℎ ∈ Ψ ℎ , ( 15
)
⟨𝜂 𝑡 ℎ , 𝑢 ℎ ⟩ 𝜕𝒯 ℎ = ⟨𝜂 𝑡 ℎ , 𝑢 𝐷 ⟩, ∀𝜂 𝑡 ℎ ∈ Λ 𝑡 ℎ , ( 16
) (𝜉, 𝑝 ℎ ) = 0, ∀𝜉 ∈ R. ( 17
)
Remark 2. As for the original weak formulation (6)-( 9), equation ( 17) is a solvability constraint, meaning that 𝑝 ℎ has zero mean over Ω. In principle, we could formulate HM-

H(div)(𝒮 ℎ ) without the multiplier 𝜌 ℎ , searching 𝑝 ℎ directly in Ψ 0ℎ = Ψ ℎ ∩ 𝐿 2 0 (Ω). In fact, equation (15) implies 𝜌 ℎ = 0, for taking constant 𝑞 ℎ = 𝜌 ℎ we obtain ||𝜌 ℎ || 2 𝐿 2 = -𝑏(𝑢 ℎ , 𝜌 ℎ ) = -⟨𝑢 𝐷 • 𝑛, 𝜌 ℎ ⟩ = 0.
In practice, the introduction of 𝜌 ℎ is for sake of simplifying coding, to avoid the construction of pressure shape functions with vanishing mean. Consequently, equation [START_REF] Guzmán | A family of nonconforming elements for the Brinkman problem[END_REF] ensures that for the divergence-free case ∇ • 𝑢 = 0 (i.e., 𝑔 = 0) this property also holds for the velocity approximation 𝑢 ℎ .

Even though the original formulation ( 10)-( 13) is well posed, this may not be true for the HM-H(div)(𝒮 ℎ ) formulation. In fact, the FE spaces of each field cannot be chosen independently one from the other: there must be compatibility relations between them. Following Brezzi's theory, as applied in [START_REF] Raviart | Primal hybrid finite element methods for 2nd order elliptic equations[END_REF][START_REF] Roberts | Mixed and hybrid methods[END_REF][START_REF] Girault | Finite element methods for Navier-Stokes equations: theory and algorithms[END_REF] for similar contexts, let us introduce the following subspaces associated to the pairs 15) and ( 16) that

{𝑋 ℎ , Ψ ℎ } and {𝑋 ℎ , Λ 𝑡 ℎ }: 𝑍 ℎ (𝑔) = {𝑣 ∈ 𝑋 ℎ ; 𝑏(𝑣, 𝑞) = (𝑔, 𝑞), ∀𝑞 ∈ Ψ ℎ }, 𝑋 ℎ (𝑢 𝐷 ) = {𝑣 ℎ ∈ 𝑋 ℎ ; ⟨𝜂 𝑡 ℎ , 𝑣 ℎ ⟩ 𝜕𝒯 ℎ = ⟨𝜂 𝑡 ℎ , 𝑢 𝐷 ⟩, ∀𝜂 𝑡 ℎ ∈ Λ 𝑡 ℎ }. It is clear from (
𝑢 ℎ ∈ 𝑋 ℎ (𝑢 𝐷 ) ∩ 𝑍 ℎ (𝑔).
According to Brezzis's theory, the HM-H(div)(𝒮 ℎ ) formulation is well posed provided the following conditions are satisfied:

(S1) inf 𝑎 𝒯 ℎ (𝑣, 𝑣) 1/2 > 0 whenever 𝑣 ∈ 𝑍 ℎ (0) ∩ 𝑋 ℎ (0), ‖𝑣‖ 𝐻(div,Ω) = 1.
(S2) There exists a positive constant 𝑐 2 such that for each 𝑞 ∈ Ψ ℎ and 𝜂 𝑡 ∈ Λ 𝑡 ℎ there is a nonzero 𝑣 ∈ 𝑋 ℎ such that

𝑏(𝑣, 𝑞) + ⟨𝜂 𝑡 , 𝑣⟩ 𝜕𝒯 ℎ ≥ 𝑐 2 ‖𝑣‖ 𝐻(div,Ω) (‖𝑞‖ 𝐿 2 (Ω) + ‖𝜂 𝑡 ‖ Λ ). ( 18
)
In fact, (S1) holds for all 𝑣 ∈ 𝑋 ℎ . It is clearly true for Brinkman problems. In the Stokes limit, suppose that, for some 𝑣 ∈ 𝑋 ℎ , 𝑎 𝒯 ℎ (𝑣, 𝑣) = ∑ 𝐾∈𝒯 ℎ (2𝜇𝐷(𝑣), 𝐷(𝑣)) 𝐾 = 0. Then 𝑣 is piecewise defined by rigid body motions, i.e. 𝐷(𝑣)| 𝐾 = 0 over 𝐾 ∈ 𝒯 ℎ , which is not possible for fields 𝑣 ̸ = 0 in 𝐻(div, Ω).

Concerning (S2), the condition [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF] in the absence of the trace terms (e.g., taking 𝜂 𝑡 = 0) is a property verified for classic stable FE pairs 𝑋 ℎ × Ψ ℎ used for stable Darcy's mixed formulations, verifying the inf-sup condition sup

𝑣 ℎ ∈𝑋 ℎ 𝑏(𝑣 ℎ , 𝑞 ℎ ) ≥ 𝐶 ‖𝑣‖ 𝐻(div,Ω) ||𝑞 ℎ || 𝐿 2 (Ω) , ∀𝑞 ℎ ∈ Ψ ℎ , ( 19
)
which is as a consequence of the divergence-compatibility for velocity and pressure spaces

∇ • 𝑋 ℎ Ψ ℎ . ( 20
)
However, as it shall become clear in the forthcoming sections, to consider the trace terms in the condition [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF], the choice of the tangential traction space Λ 𝑡 ℎ ⊂ Λ 𝑡 is a more subtle task. For instance, the following trace-compatibility property becomes crucial:

Λ 𝑡 ℎ (0) = { 𝜂 𝑡 ∈ Λ 𝑡 ℎ ; ⟨𝜂 𝑡 ℎ , 𝑣⟩ 𝜕𝒯 ℎ = 0, ∀𝑣 ∈ 𝑍 ℎ (0) } = {0}. (21)
3.1. Solvability Theorem 3. In addition to the divergence-compatibility property [START_REF] Araya | Multiscale hybrid-mixed method for the Stokes and Brinkman equations -the method[END_REF], suppose the trace compatibility condition [START_REF] John | On the divergence constraint in mixed finite element methods for incompressible flows[END_REF] holds. Then the HM-H(div)(𝒮 ℎ ) formulation has a unique solution.

Proof. Thus, taking 𝑓 = 0, 𝑔 = 0, and 𝑢 𝐷 = 0, we have to evaluate the functions (𝑢 ℎ , 𝜆 𝑡 ℎ ) ∈ Z ℎ (0) × Λ 𝑡 ℎ solving the homogeneous system

𝑎 𝒯 ℎ (𝑢 ℎ , 𝑣 ℎ ) + ⟨𝜆 𝑡 ℎ , 𝑣 ℎ ⟩ 𝜕𝒯 ℎ = 0, ∀𝑣 ℎ ∈ Z ℎ (0) ⟨𝜂 𝑡 ℎ , 𝑢 ℎ ⟩ 𝜕𝒯 ℎ = 0, ∀𝜂 𝑡 ℎ ∈ Λ 𝑡 ℎ ,
Following the algebraic aspects of sadle point problems in Chapter 3 of [START_REF] Boffi | Mixed and hybrid finite element methods[END_REF], let us expand the fields in terms of shape functions {𝑣 𝑗 } of 𝑍 ℎ (0) and {𝜂 𝑡 ℓ } of Λ 𝑡 ℎ and approach this problem in the matrix form

[ A C 𝑇 C 0 ] [ û λ𝑡 ] = [ 0 0 ] . ( 22
)
The vector unknowns û and λ𝑡 represent the degrees of freedom of 𝑢 ℎ and 𝜆 𝑡 ℎ with respect to the corresponding shape functions. The matrices are A = [(𝑎 𝒯 ℎ (𝑣 𝑖 , 𝑣 𝑗 )] and

C 𝑇 = [⟨𝜂 𝑡 ℓ , 𝑣 𝑗 ⟩ 𝜕𝒯 ℎ ],
and the vectors on the right hand side are f = [(f , 𝑣 𝑗 )] and û𝐷 = [⟨𝜂 𝑡 ℓ , 𝑢 𝐷 ⟩] 𝑇 . The second equation in [START_REF] Egger | hp analysis of a hybrid DG method for Stokes flow[END_REF] says that û ∈ K𝑒𝑟 C. In such case, multiplying the first equation by û𝑇 and observing that û𝑇 C 𝑇 λ𝑡 = ( λ𝑡 ) 𝑇 Cû = 0, we obtain û𝑇 Aû = 0. Since A is positive definite, we oconclude that û = 0. Consequently, the first equation reduces to C 𝑇 λ𝑡 = 0. The question is: C 𝑇 λ𝑡 = 0 ⇒ λ𝑡 = 0? This is true if C 𝑇 defines an injective mapping, which is the algebraic expression of the trace compatibility condition [START_REF] John | On the divergence constraint in mixed finite element methods for incompressible flows[END_REF].

Remark 3. The fact that 𝐶 𝑇 defines an injective mapping can also be expresses by the property rank 𝐶 𝑇 = dim Λ 𝑡 ℎ , from where the term rank condition comes from, according to [START_REF] Babuska | Mixed-hybrid finite element approximations of second-order elliptic boundary-value problem[END_REF]. As discussed in the Chapter 3 of [START_REF] Boffi | Mixed and hybrid finite element methods[END_REF], this solvability condition is not sufficient to provide a good method and stability is actually needed. There is where the inf-sup condition (S2) comes into play.

Error estimates

There is a classic methodology for general nonconforming mixed methods for Stokes problems, e.g., see [START_REF] Crouzeix | Conforming and non-conforming finite element methods for solving the stationary Stokes equations[END_REF][START_REF] Han | Nonconforming elements in the mixed finite element method[END_REF]. We follow a similar approach, which is also closely related to the nonconforming formulations adopted in [START_REF] Tai | A discrete de rham complex with enhanced smoothness[END_REF][START_REF] Xie | Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models[END_REF][START_REF] Guzmán | A family of nonconforming elements for the Brinkman problem[END_REF] for Stokes-Brinkman problems using 𝐻(div)-conforming velocity approximations (see Section 4.2). Theorem 4. Let (𝑢, 𝑝, 𝜆 𝑡 , 0) ∈ 𝒮 × R be the exact solution of formulation (10)-( 13), and suppose (𝑢 ℎ , 𝑝 ℎ , 𝜆 𝑡 ℎ , 0) ∈ 𝒮 ℎ × R solves the HM-H(div)(𝒮 ℎ ) system. Then,

||𝑢 -𝑢 ℎ || 𝑎 𝒯 ℎ [ inf 𝑧 ℎ 𝑍 ℎ (𝑔) ||𝑢 -𝑧 ℎ || 𝑎 𝒯 ℎ + sup 𝑣 ℎ ∈𝑍 ℎ (0)∖{0} ⟨𝜆 𝑡 ℎ -𝜆 𝑡 , 𝑣 ℎ ⟩ 𝜕𝒯 ℎ ||𝑣 ℎ || 𝑎 𝒯 ℎ ] , ( 23
)
||𝑝 -𝑝 ℎ || 𝐿 2 [ inf 𝑞 ℎ ∈Ψ 0ℎ ||𝑝 -𝑞 ℎ || 𝐿 2 + inf 𝑧 ℎ 𝑍 ℎ (𝑔) ||𝑢 -𝑧 ℎ || 𝑎 𝒯 ℎ + sup 𝑣 ℎ ∈𝑋 ℎ ∖{0} ⟨𝜆 𝑡 ℎ -𝜆 𝑡 , 𝑣 ℎ ⟩ 𝜕𝒯 ℎ ||𝑣 ℎ || 𝑎 𝒯 ℎ ] . ( 24
)
Proof. We already know that 𝑢 ∈ 𝐻 1 (Ω), ∇ • 𝑢 = 𝑔, 𝑢| 𝜕Ω = 𝑢 𝐷 , and 𝑝 ∈ Ψ 0 . In the discrete context, 𝑢 ℎ ∈ 𝑋 ℎ (𝑢 𝐷 ) ∩ 𝑍 ℎ (𝑔) and 𝑝 ℎ ∈ Ψ 0ℎ . Consider the errors 𝑢 -𝑢 ℎ , 𝑝 -𝑝 ℎ and 𝜆 𝑡 -𝜆 𝑡 ℎ . Then

𝑎 𝒯 ℎ (𝑢 -𝑢 ℎ , 𝑣 ℎ ) -𝑏(𝑣 ℎ , 𝑝 -𝑝 ℎ ) + ⟨𝜆 𝑡 -𝜆 𝑡 ℎ , 𝑣 ℎ ⟩ 𝜕𝒯 = 0 ∀𝑣 ℎ ∈ 𝑋 ℎ , ( 25
)
𝑏(𝑢 -𝑢 ℎ , 𝑞 ℎ ) = 0, ∀𝑞 ℎ ∈ Ψ 0ℎ . ( 26
)
For an arbitrary

𝑧 ℎ ∈ 𝑍 ℎ (𝑔), set 𝑣 ℎ = 𝑧 ℎ -𝑢 ℎ ∈ 𝑍 ℎ (0). Using 𝑣 ℎ as test function in (25), we obtain 𝑎 𝒯 ℎ (𝑢 -𝑢 ℎ , 𝑣 ℎ ) + ⟨𝜆 𝑡 -𝜆 𝑡 ℎ , 𝑣 ℎ ⟩ 𝜕𝒯 = 0, which can be expressed as 𝑎 𝒯 ℎ (𝑧 ℎ -𝑢, 𝑣 ℎ ) + ⟨𝜆 𝑡 ℎ -𝜆 𝑡 , 𝑣 ℎ ⟩ 𝜕𝒯 = 𝑎 𝒯 ℎ (𝑣 ℎ , 𝑣 ℎ ) = ||𝑣 ℎ || 2 𝑎 𝒯 ℎ . Therefore, ||𝑣 ℎ || 2 𝑎 𝒯 ℎ ≤ ||𝑧 ℎ -𝑢|| 𝑎 𝒯 ℎ || ||𝑣 ℎ || 𝑎 𝒯 ℎ + ⟨𝜆 𝑡 ℎ -𝜆 𝑡 , 𝑣 ℎ ⟩ 𝜕𝒯 . Thus, for 𝑣 ℎ ̸ = 0 we obtain ||𝑣 ℎ || 𝑎 𝒯 ℎ ≤ ||𝑧 ℎ -𝑢|| 𝑎 𝒯 ℎ + ⟨𝜆 𝑡 ℎ -𝜆 𝑡 , 𝑣 ℎ ⟩ 𝜕𝒯 ||𝑣 ℎ || 𝑎 𝒯 ℎ .
Consequently,

||𝑢 -𝑢 ℎ || 𝑎 𝒯 ℎ ≤ ||𝑢 -𝑧 ℎ || 𝑎 𝒯 ℎ + ||𝑧 ℎ -𝑢 ℎ || 𝑎 𝒯 ℎ = ||𝑢 -𝑧 ℎ || 𝑎 𝒯 ℎ + ||𝑣 ℎ || 𝑎 𝒯 ℎ ≤ 2||𝑧 ℎ -𝑢|| 𝑎 𝒯 ℎ + ⟨𝜆 𝑡 ℎ -𝜆 𝑡 , 𝑣 ℎ ⟩ 𝜕𝒯 ||𝑣 ℎ || 𝑎 𝒯 ℎ ,
and we obtain estimate [START_REF] Igreja | Stabilized velocity and pressure mixed hybrid DGFEM for the Stokes problem[END_REF] Let Π 𝑝 ℎ 𝑝 ∈ Ψ 0ℎ be the 𝐿 2 -projection of 𝑝 onto Ψ 0ℎ . Thus, from [START_REF] Araya | Multiscale hybrid-mixed method for the Stokes and Brinkman equations -the method[END_REF], 𝑏(𝑣 ℎ , 𝑝 -Π 𝑝 ℎ 𝑝) = 0 for arbitrary 𝑣 ℎ ̸ = 0 ∈ 𝑋 ℎ . By the inf-sup condition [START_REF] Arnold | Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates[END_REF] and equations ( 25)-( 26) we obtain

||Π 𝑝 ℎ 𝑝 -𝑝 ℎ || 𝐿 2 𝑏(𝑣 ℎ , Π 𝑝 ℎ 𝑝 -𝑝 ℎ ) ||𝑣 ℎ || 𝑎 𝒯 ℎ = 𝑏(𝑣 ℎ , 𝑝 -𝑝 ℎ ) ||𝑣 ℎ || 𝑎 𝒯 ℎ = 𝑎 𝒯 ℎ (𝑢 ℎ -𝑢, 𝑣 ℎ ) + ⟨𝜆 𝑡 ℎ -𝜆 𝑡 , 𝑣 ℎ ⟩ 𝜕𝒯 ||𝑣 ℎ || 𝑎 𝒯 ℎ ≤ ||𝑢 ℎ -𝑢|| 𝑎 𝒯 ℎ + ⟨𝜆 𝑡 ℎ -𝜆 𝑡 , 𝑣 ℎ ⟩ 𝜕𝒯 ||𝑣 ℎ || 𝑎 𝒯 ℎ .
Thefore, estimate [START_REF] Cockburn | A comparison of HDG methods for Stokes flow[END_REF] holds by using estimate [START_REF] Igreja | Stabilized velocity and pressure mixed hybrid DGFEM for the Stokes problem[END_REF] and

||𝑝 -Π 𝑝 ℎ 𝑝|| 𝐿 2 = inf 𝑞 ℎ ∈Ψ ℎ ||𝑝 -𝑞 ℎ || 𝐿 2 .
Notice that:

(i) the best approximation error in the first term on the right-hand side of ( 23) measures how well the exact velocity 𝑢 can be approximated by nonconforming functions in 𝑍 ℎ (𝑔). The second term, which would be zero for 𝐻 1 -conforming velocity approximations, is known as consistency error, measuring the nonconformity of 𝑢 ℎ . One would expect to have a method giving a consistency error of at least the same order of the first optimal error;

(ii) the best approximation of pressure 𝑝 in Ψ 0ℎ , expressed in the first term on the righthand side of [START_REF] Cockburn | A comparison of HDG methods for Stokes flow[END_REF], may be dominated by the second an third terms concening velocity best approximation and consistency error.

Examples of HM-H(div)(𝒮 ℎ ) formulations

The discussion of the previous section reveals that for the solvability of the HM-H(div)(𝒮 ℎ ) formulation the FE space setting 𝒮 ℎ = 𝑋 ℎ × Ψ ℎ × Λ 𝑡 ℎ is required to verify the following compatibility conditions:

1. The pair {𝑋 ℎ , Ψ ℎ } is divergence-compatible according to [START_REF] Araya | Multiscale hybrid-mixed method for the Stokes and Brinkman equations -the method[END_REF]. This is the property required for FE pairs applied for Darcy's flow simulations. 2. The pair {𝑋 ℎ , Λ 𝑡 ℎ } is trace compatible under the circunstances of Theorem 3. The verification of trace compatibility conditions occuring in hybrid methods is usually a rather technical issue, requiring an extra effort, as revealed in the analysis of the early primal hybrid method for Poisson problems considered in [START_REF] Raviart | Primal hybrid finite element methods for 2nd order elliptic equations[END_REF].

Our purpose here is to discuss solvability aspects for some examples where 𝑋 ℎ × Ψ ℎ are classic divergence-compatible FE pairs {𝑋 ℎ , Ψ ℎ } based on partitions 𝒯 ℎ whose elements 𝐾 are geometricaly mapped from a master element K of one of the usual geometry:

• For triangular elements 𝐾: T = {(x, ŷ); x ≥ 0, ŷ ≥ 0, x + ŷ ≤ 1}.

• For quadrilateral elements

𝐾: R = [-1, 1] × [-1, 1]. • For tetrahedra 𝐾: T 𝑒 = {(x, ŷ, ẑ); x ≥ 0, ŷ ≥ 0, ẑ ≥ 0, x + ŷ + ẑ ≤ 1}. • For hexahedra 𝐾: Ĥ = [-1, 1] × [-1, 1] × [-1, 1].
Approximation spaces 𝑋( K) and Ψ( K) defined in the master element are supposed to be divergence-compatible polynomial spaces, such that ∇•𝑋( K) = Ψ( K). Moreover, the direct sum decomposition 𝑋( K) = 𝑋 𝜕 ( K)⊕ X( K) is considered in terms of internal velocity fields X( K), with vaninshing normal components over 𝜕 K, and the external velocities 𝑋 𝜕 ( K) otherwise. Given the geometric invertible map 𝐹 𝐾 : K → 𝐾, the FE pair 𝑋( K) × Ψ( K) is mapped to 𝐾 to form FE spaces 𝑋(𝐾) × Ψ(𝐾) by the following pullback operators:

• Scalar FE functions 𝑝 ∈ Ψ(𝐾): 𝑝 = ℱ(p) := p ∘ 𝐹 -1 𝐾 , p ∈ Ψ( K).
Recall that, ℱ is an isomorphism from 𝐿 2 ( K) onto 𝐿 2 (𝐾) (and also from 𝐻 1 ( K) onto 𝐻 1 (𝐾)).

• Vector FE functions 𝑣 ∈ 𝑋(𝐾) ⊂ 𝐻(div, 𝐾): 𝑣

= 𝒢 𝐾 (v) := [ 1 J 𝐾 𝐷𝐹 𝐾 v] ∘ 𝐹 -1 𝐾 , v ∈ 𝑋( K)
, where 𝐷𝐹 𝐾 is the Jacobian matrix of 𝐹 𝐾 , and J 𝐾 = |det(𝐷𝐹 𝐾 )| (constant for affine elements).

The mapping 𝒢 𝐾 , known as Piola transformation, is an isomorphism from 𝐻(div; K) onto 𝐻(div; K), verifying (∇

• 𝒢 𝐾 (v), ℱ(p)) 𝐾 = (∇ • v, p) K , and ⟨𝒢 𝐾 (v) • 𝑛 𝐾 , ℱ(p)⟩ 𝜕𝐾 = ⟨v • 𝑛 K , p⟩ 𝜕 K .
Accordingly, the action of 𝒢 𝐾 preserves the divergence compatibility property and the direct sum decomposition, so that ∇ • 𝑋(𝐾) = Ψ(𝐾) and 𝑋(𝐾) = 𝑋 𝜕 (𝐾) ⊕ X(𝐾).

Let Γ ℎ be the skeleton mesh formed by the facets 𝐸 ⊂ 𝜕𝐾 (faces in 3D elements or edges in 2D cases). Similarly, given geometric invertible maps 𝐹 𝐸 : F → 𝐸 and vector trace spaces 𝑊 𝑡 ( Ê), vector trace FE spaces 𝑊 𝑡 (𝐸) are defined in 𝐸 by the pullback operator:

• Vector FE functions 𝜂 𝑡 ∈ 𝑊 𝑡 (𝐸), 𝐸 ∈ Γ ℎ : 𝜂 𝑡 = 𝒢 𝐸 ( η𝑡 ) := [ 1 J 𝐸 𝐷𝐹 𝐸 η𝑡 ] ∘ 𝐹 -1 𝐸 , η𝑡 ∈ 𝑊 𝑡 ( Ê).
Finally, finite dimensional FE approximation spaces 𝑋 ℎ ⊂ 𝑋, Ψ ℎ ⊂ Ψ are piecewise defined over 𝒯 ℎ by the usual assembly process of the local spaces 𝑋(𝐾), Ψ(𝐾). Similarly, 1 for master elements of triangular ( T ), quadrilateral ( R), tetrahedral ( T 𝑒) or hexahedral ( Ĥ) geometry. The following polynomial spaces are used in K: scalar polynomials P 𝑘 ( K) of total degree at most 𝑘, P𝑘 ( K) denoting the homogeneous ones; vector versions P 𝑘 ( K, R 𝑑 ) are also used; the polynomials in Q 𝑘,𝑚 ( K) or Q 𝑘,𝑚,𝑛 ( K) have maximum degree 𝑘, 𝑚 or 𝑛 in each coordinate.

Λ 𝑡 ℎ ⊂ Λ 𝑡 is piecewise defined over Γ ℎ as Λ 𝑡 ℎ = {𝜂 𝑡 ∈ Λ 𝑡 ; 𝜂 𝑡 | 𝐸 ∈ 𝑊 𝑡 (𝐸), 𝐸 ∈ Γ ℎ }. Some classic examples of divergence-compatible FE pairs 𝑋( K) × Ψ( K) are displayed in Table
K Method X( K) Ψ( K)
T 𝐵𝐷𝑀 (𝑘) [START_REF] Brezzi | Two families of mixed finite elements for second order elliptic problems[END_REF] P 𝑘 ( K, R 2 ) P 𝑘-1 ( K) 𝑅𝑇 (𝑘) [START_REF] Raviart | A mixed finite element method for 2-nd order elliptic problems[END_REF] P 𝑘 ( K, R 2 ) ⊕ x P𝑘 ( K) To be combined with X( K) × Ψ( K), let us consider multiplier vector polynomials of the following types:

P 𝑘 ( K) 𝐵𝐷𝑀 + (𝑘) [33] P 𝜕 𝑘 ( K, R 2 ) ⊕ P𝑘+1 ( K, R 2 ) P 𝑘 ( K) R 𝑅𝑇 (𝑘) [32] Q 𝑘+1,𝑘 × Q 𝑘,𝑘+1 ( K) Q 𝑘,𝑘 ( K) 𝑅𝑇 + (𝑘) [34] 𝑋 𝜕 𝑅𝑇 (𝑘) ( K) ⊕ X𝑅𝑇 (𝑘+1) ( K) Q 𝑘+1,𝑘+1 ( K) T 𝑒 𝐵𝐷𝑀 (𝑘) [35] P 𝑘 ( K, R 3 ) P 𝑘-1 ( K) 𝑅𝑇 (𝑘) [36] P 𝑘 ( K, R 3 ) ⊕ x P𝑘 ( K) P 𝑘 ( K) 𝐵𝐷𝑀 + (𝑘) [33] P 𝜕 𝑘 ( K, R 3 ) ⊕ P𝑘+1 ( K, R 3 ) P 𝑘 ( K) Ĥ 𝑅𝑇 (𝑘) [36] Q 𝑘+1,𝑘,𝑘 × Q 𝑘,𝑘+1,𝑘 × Q 𝑘,𝑘,𝑘+1 ( K) Q 𝑘,𝑘,𝑘 ( K) 𝑅𝑇 + (𝑘) [37] 𝑋 𝜕 𝑅𝑇 (𝑘) ( K) ⊕ X𝑅𝑇 (𝑘+1) ( K) Q 𝑘+1,𝑘+1,𝑘+1 ( K)
• For edges: 𝑊 𝑡 ( Ê) = P 𝑚 ( Ê)𝜏 Ê on edges Ê;

• For triangular facets: 𝑊 𝑡 ( Ê) = n ∧ P 𝑚 ( Ê, R 3 );

• For quadrilateral facets: 𝑊 𝑡

( Ê) = n ∧ Q 𝑚,𝑚 ( Ê, R 3 ).
The parameter 𝑚 has to be be chosen properly for the verification of the trace compatibility property [START_REF] John | On the divergence constraint in mixed finite element methods for incompressible flows[END_REF]. In the sequence we discuss a numerical test for the verification of the trace compatibility condition for the FE spaces based on triangular an quadrilateral elements.

Verification of the trace compatibility condition (21)

The purpose is to verify in which extent the divergence-free subspaces 𝑍 ℎ (0) ⊂ 𝑋 ℎ associated to some examples shown in Table 1 are rich enough to enforce a trace compatibility condition [START_REF] John | On the divergence constraint in mixed finite element methods for incompressible flows[END_REF] with respect to the multiplier space Λ 𝑡 ℎ , by properly choosing the polynomial degree 𝑚. As remarked in the proof of Theorem 3, for shape functions {𝜂 𝑡 ℓ } of Λ 𝑡 ℎ and {𝑣 𝑗 } of 𝑍 ℎ (0), the matrix C 𝑇 = [⟨𝜂 𝑡 ℓ , 𝑣 𝑗 ⟩ 𝜕𝒯 ℎ ] is required to be injective, which can be expresses by the rank condition rank

C 𝑇 = dim Λ 𝑡 ℎ .
For that, we propose a numerical test of rank condition verification for regions Ω = 𝐾 0 ∪𝐾 1 , the elements in 𝒯 ℎ = {𝐾 0 , 𝐾 1 } being triangles or rectangles, as represented in Figure 1. In each case, the characterization of 𝑍 ℎ (0) is crucial, which can be obtained in the context of the exact FE sequences shown in Table 2, where for scalar functions ∇ × 𝜙 = [𝜕𝜙/𝜕𝑦, -𝜕𝜙/𝜕𝑥] 𝑇 . They are discrete versions of the exact differential sequence 𝐻 1

( K) ∇× -→ H(div, K) ∇•
-→ 𝐿 2 ( K), meaning that each space is mapped into the succeeding space by the given differential operator, the composition of two consecutive maps being zero. In such context, if the subspace 𝑈 ( K) ⊂ 𝐻 1 ( K) combines with the space 𝑋( K) ⊂ H(div, K), then the corresponding divergence-free functions are in ∇×𝑈 ( K) ⊂ 𝑋( K). Thus, the cardinality of a given basis ℬ grad ℎ of 𝑈 ℎ is one unit more than the dimension of 𝑍 ℎ (0). Under this point of view, we define shape functions ℬ grad,filt ℎ by removing one arbitrary vertex function from a basis ℬ grad ℎ for 𝑈 ℎ and set ℬ div,0 ℎ = ∇ × ℬ grad,filt ℎ for a basis for 𝑍 ℎ (0). Taking the traction spaces Λ 𝑡 𝛾 with 𝑚 = 𝑘 -1, the rank of C 𝑇 and the dimension of Λ 𝑡 ℎ are presented in Table 3. Comparing the results, trace compatibility ( 21) is confirmed for 𝑅𝑇 (𝑘), 𝑅𝑇 + (𝑘) and 𝐵𝐷𝑀 + (𝑘) cases. However, it is not verified for 𝐵𝐷𝑀 (𝑘).

Space

𝐵𝐷𝑀 (𝑘) 𝐵𝐷𝑀 + (𝑘) 𝑈 ( T ) ⊂ 𝐻 1 ( T ) P 𝑘+1 ( T ) P 𝜕 𝑘+1 ( T ) ⊕ P𝑘+2 ( T ) ↓ ∇× ↓ ∇× ↓ ∇× 𝑋( T ) ⊂ H(div, T ) P 𝑘 ( T , R 2 ) P 𝜕 𝑘 ( T , R 2 ) ⊕ P𝑘+1 ( T , R 2 ) ↓ ∇• ↓ ∇• ↓ ∇• Ψ( T ) ⊂ 𝐿 2 ( T ) P 𝑘-1 ( T ) P 𝑘 ( T ) Space 𝑅𝑇 (𝑘) 𝑅𝑇 + (𝑘) 𝑈 ( R) ⊂ 𝐻 1 ( R) Q 𝑘+1,𝑘+1 ( R) Q 𝜕 𝑘+1,𝑘+1 ( R) ⊕ Q𝑘+2,𝑘+2 ( R) ↓ ∇× ↓ ∇× ↓ ∇× 𝑋( R) ⊂ H(div, R) 𝑋 𝑅𝑇 (𝑘) ( R) 𝑋 𝜕 𝑅𝑇 (𝑘) ( R) ⊕ X𝑅𝑇 (𝑘+1) ( R) ↓ ∇• ↓ ∇• ↓ ∇• Ψ( R) ⊂ 𝐿 2 ( R) Q 𝑘,𝑘 ( R) Q 𝑘+1,𝑘+1 ( R)

Hybridization of nonconforming methods

Another approach for the Stokes-Brinkman problem (1)-( 3) based on divergence-compatible FE pairs 𝑋 ℎ × Ψ ℎ considers the nonconforming formulation (without hybridization).

Problem (P h ): find ( ũℎ , pℎ ) ∈ 𝑋 ℎ (𝑢 𝐷 ) × Ψ 0ℎ such that 𝑎 ℎ ( ũℎ , 𝑣 ℎ ) -𝑏(𝑣 ℎ , pℎ ) = (𝑓 , 𝑣 ℎ ) ∀𝑣 ℎ ∈ 𝑋 ℎ (0), ( 27
)
𝑏( ũℎ , 𝑞 ℎ ) = (𝑔, 𝑞 ℎ ) ∀𝑞 ℎ ∈ Ψ 0ℎ . ( 28
)
We already know that if (𝑢 ℎ , 𝑝 ℎ , 𝜆 𝑡 , 0) ∈ 𝒮 ℎ × R solves of the HM-H(div)(𝒮 ℎ ) formulation then 𝑢 ℎ ∈ 𝑋 ℎ (𝑢 𝐷 ) and 𝑝 ℎ ∈ Ψ 0ℎ . Moreover, by testing equation ( 14) with 𝑣 ℎ ∈ 𝑋 ℎ (0) and equation [START_REF] Guzmán | A family of nonconforming elements for the Brinkman problem[END_REF] with 𝑞 ℎ ∈ Ψ 0ℎ , we conclude that ( ũℎ , pℎ ) = (𝑢 ℎ , 𝑝 ℎ ) solves Problem (P h ). In this section we show that the converse is also true for solutions of some known well posed cases of Problem (P h ).

Having in mind the construction of well posed FE settings 𝑋 ℎ (𝑢 𝐷 ) × Ψ 0ℎ for Problem (P h ), some first order methods based on simplicial meshes (for triangles or tetrahedra) were derived in [START_REF] Mardal | A robust finite element method for Darcy-Stokes flow[END_REF][START_REF] Tai | A discrete de Rham complex with enhanced smoothness[END_REF][START_REF] Xie | Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models[END_REF], and extended in [START_REF] Guzmán | A family of nonconforming elements for the Brinkman problem[END_REF] to higher order cases. First order examples for affine quadrilateral and hexaheral elements were proposed in [START_REF] Chen | Uniformly convergent H(div)-conforming rectangular elements for Darcy-Stokes problem[END_REF][START_REF] Chen | Uniformly convergent cubic nonconforming element for Darcy-Stokes problem[END_REF]. In general terms, the approach of these methods is to combine the following tools:

1. There is a basic divergence-compatible FE pair 𝑉 (𝐾)×Ψ(𝐾). For simplicial elements, the examples in [START_REF] Tai | A discrete de Rham complex with enhanced smoothness[END_REF][START_REF] Xie | Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models[END_REF] are for the FE pair 𝐵𝐷𝑀 (1), and in [START_REF] Guzmán | A family of nonconforming elements for the Brinkman problem[END_REF] they are for 𝐵𝐷𝑀 (𝑘) or 𝑅𝑇 (𝑘), 𝑘 ≥ 1. The cases in [START_REF] Chen | Uniformly convergent H(div)-conforming rectangular elements for Darcy-Stokes problem[END_REF][START_REF] Chen | Uniformly convergent cubic nonconforming element for Darcy-Stokes problem[END_REF] are for the 𝐵𝐷𝑀 (1) elements in quadrilateral or hexaedral elements (see [START_REF] Boffi | Mixed and hybrid finite element methods[END_REF]). In all these cases, the velocity normal traces are of degree 𝑘. 2. There is a tangential FE trace space 𝑊 𝑡 (𝜕𝐾) defined by piecewise vector polynomials of degree 𝑚 = 𝑘 -1.

A complementary divergence-free bubble polynomial space is constructed in the form 𝐻(𝐾) = curl (𝑏 𝐾 𝑄(𝐾)).

The function 𝑏 𝐾 is the usual scalar bubble function associated to 𝐾, and 𝑄(𝐾) are properly chosen spaces such that the functions 𝑣 ∈ 𝐻(𝐾) are characterized by the degrees of freedom defined by the functionals ⟨ η𝑡 , 𝑣⟩ 𝜕𝐾 , 𝜂 𝑡 ∈ 𝑊 𝑡 (𝜕𝐾).

A new pair 𝑋(𝐾) × Ψ(𝐾) is defined such that 𝑋(𝐾) = 𝑉 (𝐾) ⊕ 𝐻(𝐾).

As discussed in [START_REF] Xie | Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models[END_REF], uniform stability of these FE pairs 𝑋 ℎ (𝑢 𝐷 ) × Ψ 0ℎ holds both with respect to the parameters 𝛼 and ℎ for Problem (P h ) with respect to the norm

|||𝑢||| 2 ℎ := 𝑎 ℎ (𝑢, 𝑣) + 𝑀 ||∇ • 𝑢|| 2 𝐿 2 (Ω) ,
where 𝑀 = max{𝛼, 𝜇}.

Lemma 2. Let 𝒮 ℎ = 𝑋 ℎ × Ψ ℎ × Λ 𝑡 ℎ ⊂ 𝐻(div, Ω) × Ψ × Λ t be

a FE setting constructed by assembling local spaces 𝑋(𝐾) × Ψ(𝐾) × 𝑊 𝑡 (𝜕𝐾) verifying the aforementioned properties.

In this way, the FE pair 𝑋 ℎ × Ψ ℎ is divergence-compatible and 𝑋 ℎ × Λ 𝑡 ℎ verifies the trace compatibility condition [START_REF] John | On the divergence constraint in mixed finite element methods for incompressible flows[END_REF]. Thus, the corresponding hybridized HM-H(div)(𝒮 ℎ ) version for this method has unique solution.

Proof. In fact, the enrichment procedure 𝑋(𝐾) = 𝑉 (𝐾) ⊕ 𝐻(𝐾) by divergence-free bubble functions 𝐻(𝐾) does not affect the divergence-compatibility property and ∇•𝑋(𝐾) Ψ(𝐾) still holds. Moreover, since the degrees of freedom of 𝐻 ℎ are characterized by the functionals defined by 𝜂 𝑡 ℎ ∈ Λ 𝑡 ℎ , and recalling that 𝐻 ℎ ⊂ 𝑍 ℎ (0), then the property ⟨𝜂 𝑡 ℎ , 𝑣 ℎ ⟩ 𝜕𝒯 ℎ = 0, ∀𝑣 ℎ ∈ 𝑍 ℎ (0) implies 𝜂 𝑡 ℎ = 0, meaning that Λ 𝑡 ℎ (0) = {0}. The existence of unique solution for the corresponding hybrid-mixed version Problem HM-H(div)(𝒮 ℎ ) for 𝒮 ℎ is consequence of Theorem 3.

Static condensation: reduced global system

Assembly and resolution of the full linear system associated to mixed methods demand elevated computational cost as compared with classic 𝐻 1 -conforming methods, asking for optimized procedures to make the method competitive. Hybridization is a typical strategy adopted for application of a static condensation scheme to the system of equations. The degrees of freedom (DoF) of the problem are classifyied according to two sets of fields: the primal variables (solved globally) and the secondary variables (recovered locally from the primal ones). Usually, in FE methods the secondary variables correspond to terms that can be resolved independently within single elements (that can be condensed), and the primal ones involve contributions spread in more elements (they incorporate the global system to be solved).

For that, we introduce two new variables pℎ ∈ P and ḡℎ ∈ P, playing the roles of element average pressure and element distributed flux, repectively. This way, the HM-H(div)(𝒮 ℎ ) formulation can be expressed in the an equivalent form:

HM-H(div)(𝒮 ℎ ): find (𝑢 ℎ , 𝑝 ℎ , pℎ , ḡℎ , 𝜆 𝑡 ℎ , 𝜌 ℎ ) ∈ 𝑋 ℎ × Ψ ℎ × Ψ × Ψ × Λ 𝑡 ℎ × R satisfying 𝑎 𝒯 (𝑢 ℎ , 𝑣) + 𝑏(𝑣, 𝑝 ℎ ) + ⟨𝜆 𝑡 ℎ , 𝑣⟩ 𝜕𝒯 ℎ = (𝑓 , 𝑣), ( 29
)
𝑏(𝑢 ℎ , 𝑞) + (ḡ ℎ , 𝑞) + (𝜌 ℎ , 𝑞) = 0, ( 30
) (κ, 𝑝 ℎ -pℎ ) = 0, ( 31
)
⟨𝜂 𝑡 , 𝑢 ℎ ⟩ 𝜕𝒯 = ⟨𝜂 𝑡 , 𝑢 𝐷 ⟩, (32)
(ḡ ℎ , q) = 0, ( 33
) (𝜉, 𝑝 ℎ ) = 0, (34)
for all (𝑣, 𝑞, q, κ, 𝜂 𝑡 , 𝜉)

∈ 𝑋 ℎ × Ψ ℎ × Ψ × Ψ × Λ 𝑡 × R.
Notice that the new variable ḡℎ is a multiplier introduced to enforce the solvability constraint 𝑝 ℎ -pℎ ∈ Ψ ⊥ . For mixed formulations using 𝐻(𝑑𝑖𝑣)-conforming approximations, the decomposition 𝑋 ℎ = 𝑋 𝜕 ℎ ⊕ Xℎ is helpful. This way, we express the velocity variable in the form 𝑢 ℎ = 𝑢 𝑒 ℎ + 𝑢 𝑖 ℎ , where the internal component 𝑢 𝑖 is of secondary type, being localized in single elements, whilst 𝑢 𝑒 is a primal variable, receiving influence for both sides of an element interface. Similarly, the pressure variable can be decomposed as 𝑝 ℎ = pℎ + 𝑝 ⊥ ℎ ∈ Ψ ⊕ Ψ ⊥ ℎ , where 𝑝 ⊥ ℎ is also a secondary component and pℎ a primary variable. After rearrangements of ( 29)- [START_REF] Farias | Two dimensional mixed finite element approximations for elliptic problems with enhanced accuracy for the potential and flux divergence[END_REF], splitting secondary and primary variables, the structure of the resulting local linear system can be symbolically represented in each element 𝐾 as (where the subscript h has been omitted for simplicity):

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐴 𝑖𝑖 𝐵 𝑖𝑝 0 𝐵 𝑇 𝑖𝑝 0 M 0 M 𝑇 0 𝐴 𝑖𝑒 C 𝑖𝜂 0 0 𝐵 𝑇 𝑒𝑝 0 0 𝐷 0 0 G 0 𝐴 𝑇 𝑖𝑒 𝐵 𝑒𝑝 0 C 𝑇 𝑖𝜂 0 0 0 0 G 0 𝐷 𝑇 0 𝐴 𝑒𝑒 C 𝑒𝜂 0 0 C 𝑇 𝑒𝜂 0 0 0 0 0 0 0 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ û𝑖 p ḡ û𝑒 λ𝑡 p ρ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑅 𝑖 0 0 𝑅 𝑒 𝑅 𝜌 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ( 35
)
where 𝑉 𝑠 = ( û𝑖 , p, ḡ) and 𝑉 𝑝 = ( û𝑒 , λ𝑡 , p, ρ) represent the DoF associated to secondary and primary variables, respectively. The block-matrices are

• A 𝑖𝑖 ∈ R 𝑛 𝑖 ×𝑛 𝑖 , A 𝑖𝑒 ∈ R 𝑛 𝑖 ×𝑛 𝑒 , A 𝑒𝑒 ∈ R 𝑛 𝑒 ×𝑛 𝑒 , B 𝑖𝑝 ∈ R 𝑛 𝑖 ×𝑛 𝑝 , B 𝑒𝑝 ∈ R 𝑛 𝑒 ×𝑛 𝑝 , • B 𝑒p ∈ R 𝑛 𝑒 ×n , G ∈ R n×n , C 𝑖𝜂 ∈ R 𝑛 𝑖 ×𝑛 𝜂 , C 𝑒𝜂 ∈ R 𝑛 𝑒 ×𝑛 𝜂 , M ∈ R 𝑛 𝑝 ×n , • 𝐷 ∈ R 𝑛 𝑝 , 𝑅 𝑖 ∈ R 𝑛 𝑖 , 𝑅 𝑒 ∈ R 𝑛 𝑒 , 𝑅 𝜌 ∈ R 𝑛 𝜂 ,
where: 𝑛 𝑖 = number of internal velocity shape functions, 𝑛 𝑒 = number of border velocity shape functions, 𝑛 𝑝 = number of pressure shape functions, n = number of local average pressure shape functions, and 𝑛 𝜂 =number of multiplier shape functions. By expressing [START_REF] Brezzi | Mixed finite elements for second order elliptic problems in three variables[END_REF] in the compact block form

⎡ ⎣ 𝐾 11 𝐾 12 𝐾 21 𝐾 22 ⎤ ⎦ ⎡ ⎣ 𝑉 𝑠 𝑉 𝑝 ⎤ ⎦ = ⎡ ⎣ 𝑅 𝑠 𝑅 𝑝 ⎤ ⎦ ,
the primary variable 𝑉 𝑝 is the solution of the reduced system

[ 𝐾 22 -𝐾 21 𝐾 -1 11 𝐾 12 ] 𝑉 𝑝 = 𝑅 𝑝 -𝐾 21 𝐾 -1 11 𝑅 𝑠 .
After the computation of the global variable 𝑉 𝑝 solving the system 𝐾𝑉 𝑝 = 𝑅 obtained by the assembly of these local primary systems, the secondary variable 𝑉 𝑠 is recovered by indenpendent local solvers 𝑉 𝑠 = 𝐾 -1 11 (𝑅 𝑠 -𝐾 𝑉 𝑝 ) in each element 𝐾.

Computational implementation

Some aspects of the computational implementation for a mesh partition 𝒯 ℎ are displayed in Figure 2. The hybridization strategy requires the construction of approximation spaces for velocity, pressure and traction multipliers, including the respective DoF. And for each element 𝐾, we also include the DoF associated with the average pressure p and distributed flux ḡ. Over the mesh skeleton Γ ℎ , interface elements are created linking the DoF of 𝜆 𝑡 and 𝑢 (for velocities in each element 𝐾 side), which weakly impose the terms ⟨𝜆 𝑡 ℎ , 𝑣⟩ 𝜕𝒯 and ⟨𝜂 𝑡 , 𝑢 𝑒 ℎ + 𝑢 𝑖 ℎ ⟩ 𝜕𝒯 , respectively in Eqs. ( 29) and [START_REF] Raviart | A mixed finite element method for 2-nd order elliptic problems[END_REF]. For the boundary condition, the same procedure is adopted. A robust and flexible computational code is developed using object-oriented programming principles and tools of the open-source NeoPZ library.

Numerical verification tests

In this section, we present and discuss some numerical test models to illustrate the approximation properties of the hybrid-mixed formulation analyzed in the previous sections. The simulations were implemented in the computational framework NeoPZ1 , where tools for the construction of the required divergence consistent FE pairs are available, e .g., a hierarchy of shape functions of high degree a variety of element geometry, and data structure allowing the identification of face and internal shape functions of different degrees. The integrals required to form the linear systems are computed exactly by proper high order choice of Gauss integration rules.

Affine 2D models

The Brinkman-Stokes test problems proposed in [START_REF] Botti | A hybrid high-order discretisation of the Brinkman problem robust in the Darcy and Stokes limits[END_REF] defined in the rectangular domain Ω = (0, 2)×(-1, 1) is considered. Taking 𝐶 = 𝐶(𝛼, 𝜇) = 𝛼𝜇 -1 , and 𝜒(𝜉) = exp(𝜉) -1 , 𝜉 ∈ R + (with 𝜒(+∞) = 0), the exact solutions 𝑢 = 𝜒(𝐶) 𝑢 𝒮 + (1 -𝜒(𝐶)) 𝑢 𝒟 , 𝑝 = cos 𝑥 sin 𝑦 -𝑝 0 hold for the problems in the whole range of parameters, where 𝑝 0 ∈ R is such that the zero average 𝑝 ∈ 𝐿 2 0 (Ω) is verified. We present numerical results for the Brinkman problem with (𝜇, 𝛼) = (1, 1), and for the limit cases (𝜇, 𝛼) = (1, 0) (Stokes) and (𝜇, 𝛼) = (0, 1) (Darcy), with corresponding velocity fields

𝑢 𝒮 = ∇ × 𝜓, where 𝜓 = -sin 𝑥 cos 𝑦, ( 36
)
𝑢 𝒟 = -∇𝑝, (37)
𝑢 ℬ = 𝑒 -1 𝑢 𝒮 + (1 -𝑒 -1 ) 𝑢 𝒟 . ( 38
)
Boundary and source data are obtained from the expressions of these exact solutions. Uniform conformal meshes 𝒯 ℎ = {𝐾} of triangular or quadrilateral elements 𝐾 are considered, with mesh size ℎ = 2/𝑁, 𝑁 = 2 𝑖 , 𝑖 = 2, • • • 6. Based on them, we consider FE spaces of type H(div)(ℎ, 𝒮 ℎ,𝑅𝑇 (𝑘) ) for quadrilaterals and 𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) for triangles, 𝑘 = 1, 2, and 3, with polynomial trace spaces having degree 𝑚 = 𝑘 -1.

The 𝐿 2 -errors for 𝑢, ∇•𝑢, and 𝑝 are presented in Table 4 (top side) for the Stokes problem obtained by the formulations HM-H(div)(ℎ, 𝒮 ℎ,𝑅𝑇 (𝑘) ) (left side) and HM-H(div)(ℎ, 𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) ) (right side). Optimal convergence rates of order 𝑘 + 1 are reached for the velocity variable for all FE space configurations. The observed pressure rates are of order 𝑘 for the 𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) model, but an improved order 𝑘 + 1 2 occurs for the 𝒮 ℎ,𝑅𝑇 (𝑘) space setting. For all simulations, the 𝐿 2 -errors of ∇ • 𝑢 illustrate the strongly enforced divergence-free property.

The results for the Brinkman problem presented in Table 4 (bottom side) show similar convergence behaviour for velocity and pressure, excepting that the improved rates of order 𝑘 + 1 2 for pressure using the 𝒮 ℎ,𝑅𝑇 (𝑘) space configuration is only observed for odd 𝑘. For the Brinkman problem, optimal convergence rate of order 𝑘 + 1 is verified. In Figure 3 the effect of static condensation is verified in terms of the size of the global system to be solved with respect to the total number of unknowns for 𝑁 = 64. For all the methods, the reduction is more significant with increasing order of approximation and for rectangular meshes, for which the condensed equations amounts to about 80% when 𝑘 = 5.

𝑁 ||𝑢 -𝑢 ℎ || ||𝑝 -𝑝 ℎ || ||∇ • (𝑢 -𝑢 ℎ )|| ||𝑢 -𝑢 ℎ || ||𝑝 -𝑝 ℎ || ||∇ • (𝑢 -𝑢 ℎ )|| 𝑘 =

Comparison between HM-H(div)(𝒮 ℎ ) and DG methods

The errors for 𝑢 and 𝑝 of the Stokes problem shown in Table 4 for the HM-H(div)(𝒮 𝑅𝑇 (3) ) method are plotted in Figure 4 (top side). For comparison, corresponding error curves given by two versions of the symmetric Discontinuos Galerkin method are included as well.

The DG methods are denoted by 𝑅𝑇 𝐷𝐺 (𝑘) and 𝑅𝑇 div,𝐷𝐺 (𝑘), where the indexes 𝐷𝐺 and div, 𝐷𝐺 are used to distinguish the DG methods using fully discontinuous or H(div)conforming velocity spaces (see [START_REF] Carvalho | On the use of divergence balanced H(div)-L2 pair of approximation spaces for divergence-free and robust simulations of Stokes, coupled Stokes-Darcy and Brinkman problems[END_REF] and references therein). Notice that velocity and pressure spaces of all methods coincide on the master element, but they differ in their assembly strategies.

Comparison results for the HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (3) ) formulation and for the corresponding symmetric Discontinuos Galerkin methods 𝐵𝐷𝑀 +,𝐷𝐺 (3) and 𝐵𝐷𝑀 +,div,𝐷𝐺 (3), based on triangular meshes, are plotted in Figure 4 (botom side). The velocity convergence rates of order 𝑘 + 1 occur for all DG and HM-H(div) formulations, the error curves being quite similar. However, convergence histories for pressure improves significantly when hybridization is applied, as compared to the DG methods.

It should be noted that the matrix bandwidth of the HM-H(div) formulation is considerably smaller than the matrix bandwidth of the DG variants, leading to increased numerical efficiency. Moreover, the exact divergence-free property of velocity fields approximated in H(div)-conforming velocity FE spaces is not verified by the fully discontinuous DG methods.

Pressure-robustness verification

Pressure-robustness is an important feature of numerical methods for fluid flows, meaning that an external perturbation in the source term, which influences only the pressure in the

• • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ 2.E-2 5.E-2 1.E-1 2.E-1
1.E-5

1.E-6

1.E-7

1.E-8

1.E-9 h ||u-uh||

4 1 • RT DG (3) ■ RT div,DG (3)
◆ S h,RT (3)
• • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ 2.E-2 5.E-2 1.E-1 2.E-1 1.E-2
1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

1.E-8 h ||p-ph||

3 1 3.5 1 • RT DG (3) ■ RT div,DG (3)
◆ S h,RT (3)
Triangular meshes

• • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ 2.E-2 5.E-2 1.E-1 2.E-1 1.E-4
1.E-5

1.E-6

1.E-7

1.E-8

1.E-9 h ||u-uh|| stress equilibrium, affects only the pressure in the discrete equations as well. In order to exemplify this property for the the HM-H(div)(𝒮 ℎ ) method, the model problem obtained from [START_REF] John | On the divergence constraint in mixed finite element methods for incompressible flows[END_REF] is considered for the domain Ω = (0, 1) × (0, 1), 𝑢 = 0 and no-slip boundary conditions. We set a no-flow problem 𝑢 = 0 with 𝜇 = 1, the pressure exact solution 𝑝 = 𝑚 𝑝𝑟 (𝑦 3 -𝑦 2 /2 + 𝑦 -7/12) and the source therm 𝑓 = (0, 𝑚 𝑝𝑟 (1 -𝑦 + 3𝑦 2 )) using a parameter 𝑚 𝑝𝑟 > 0.

4 1 • BDM +,DG (3)
■ BDM +,div,DG (3)
◆ S h,BDM+ (3)
• • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ 2.E-2 5.E-2 2.E-1 1.E-2 1.E-3 1.E-4 1.E-5 1.E-6 h ||p-ph|| 3 1 • BDM +,DG (3)
■ BDM +,div,DG (3)
◆ S h,BDM+ (3)
For numerical approximations, consider uniform meshes 𝒯 ℎ = {𝐾} of triangular or quadrilateral elements 𝐾, with mesh size ℎ = 1/𝑁, 𝑁 = 2 𝑖 , 𝑖 = 2, • • • 7. The results for velocity and pressure fields obtained by using FE configurations of types 𝒮 ℎ,𝑅𝑇 (𝑘) and 𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) are plotted in Figure 5, for 𝑘 = 2 and 𝑚 = 1. The 𝐿 2 -norm errors of ∇ • 𝑢 are listed in Table 5 for all simulations.

The pressure-robustness is verified when we impose a variation in the parameter 𝑚 𝑝𝑟 ,

• • • • • • ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▲ 8.0E-3 1.6E-2 3.1E-2 6.2E-2 1.3E-1 2.5E-1 1.E-10 1.E-11 1.E-12
1.E-13

1.E-14

1.E-15 1.E-16 h ||u-uh|| RT, k = 2 • mpr = 1 ■ mpr = 100 ◆ mpr = 10⁴ ▲ mpr = 10⁶ • • • • • • ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▲ 8.0E-3 1.6E-2 3.1E-2 6.2E-2 1.3E-1 2.5E-1 1.E+3 1.E+2 1.E+1
1.E+0

1.E-1

1.E-2

1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

1.E-8 h ||p-ph||

3 1 RT, k = 2 • mpr = 1 ■ mpr = 100 ◆ mpr = 10⁴ ▲ mpr = 10⁶ • • • • • • ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▲ 8.0E-3 1.6E-2 3.1E-2 6.2E-2 1.3E-1 2.5E-1 1.E-10 1.E-11 1.E-12 1.E-13 1.E-14 1.E-15 1.E-16 1.E-17 h ||u-uh|| BDM + , k = 2 • mpr = 1 ■ mpr = 100 ◆ mpr = 10⁴ ▲ mpr = 10⁶ • • • • • • ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▲ 8.0E-3 1.6E-2 3.1E-2 6.2E-2 1.3E-1 2.5E-1 1.E+2 1.E+1
1.E+0

1.E-1

1.E-2

1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

1.E-8 h ||p-ph|| which has a direct impact in the pressure error with a significant increase in its magnitude, as expected. Without significant influence of velocity and consistency errors, the pressure errors recover the optimal convergence rate of order 𝑘 + 1. In terms of the velocity error, we observe the divergence-free property and a reduced error magnitude explained by the machine precision in the numerical computations.

3 1 BDM + , k = 2 • mpr = 1 ■ mpr = 100 ◆ mpr = 10⁴ ▲ mpr = 10⁶

Brinkman problems in Darcy's limit

A numerical verification is proposed for the Brinkman problem close the limit Darcy's regime, which can be a critical flow situation.

• • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▼ ▼ ▼ ▼ ▼ ○ ○ ○ ○ ○ □ □ □ □ □ ◇ ◇ ◇ ◇ ◇ △ △ △ △ △ 2.E-2 5.E-2 1.E-1 2.E-1
1.E-2

1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

1.E-8

1.E-9

1.E-10 h ||u-uh||

2 1 3 1 4 1 • k=1, μ=10 -5 ■ k=1, μ=10 -6 ◆ k=1, μ=10 -7 ▲ k=2, μ=10 -5 ▼ k=2, μ=10 -6 ○ k=2, μ=10 -7 □ k=3, μ=10 -5 ◇ k=3, μ=10 -6 △ k=3, μ=10 -7 • • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▼ ▼ ▼ ▼ ▼ ○ ○ ○ ○ ○ □ □ □ □ □ ◇ ◇ ◇ ◇ ◇ △ △ △ △ △ 2.E-2 5.E-2 1.E-1 2.E-1
1.E-2

1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

1.E-8

1.E-9

1.E-10 h ||p-ph||

2 1 3 1 4 1 • k=1, μ=10 -5 ■ k=1, μ=10 -6 ◆ k=1, μ=10 -7 ▲ k=2, μ=10 -5 ▼ k=2, μ=10 -6 ○ k=2, μ=10 -7 □ k=3, μ=10 -5 ◇ k=3, μ=10 -6 △ k=3, μ=10 -7
Triangular meshes

• • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▼ ▼ ▼ ▼ ▼ ○ ○ ○ ○ ○ □ □ □ □ □ ◇ ◇ ◇ ◇ ◇ △ △ △ △ △ 2.E-2 5.E-2 1.E-1 2.E-1
1.E-2

1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

1.E-8

1.E-9

1.E-10 h ||u-uh||

2 1 3 1 4 1 • k=1, μ=10 -5 ■ k=1, μ=10 -6 ◆ k=1, μ=10 -7 ▲ k=2, μ=10 -5 ▼ k=2, μ=10 -6 ○ k=2, μ=10 -7 □ k=3, μ=10 -5 ◇ k=3, μ=10 -6 △ k=3, μ=10 -7 • • • • • ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▼ ▼ ▼ ▼ ▼ ○ ○ ○ ○ ○ □ □ □ □ □ ◇ ◇ ◇ ◇ ◇ △ △ △ △ △ 2.E-2 5.E-2 1.E-1 2.E-1
1.E-2

1.E-3

1.E-4

1.E-5

1.E-6

1.E-7

1.E-8

1.E-9 h ||p-ph|| The Brinkman problem with the exact solution [START_REF] Tai | A discrete de Rham complex with enhanced smoothness[END_REF] for 𝛼 = 1 and viscosity values 𝜇 = 10 -5 , 10 -6 and 10 -7 are simulated with the methods HM-H(div)(𝒮 ℎ,𝑅𝑇 (𝑘) ) and HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) ), based on uniform quadrilateral and triangular meshes, respectively.

2 1 3 1 4 1 • k=1, μ=10 -5 ■ k=1, μ=10 -6 ◆ k=1, μ=10 -7 ▲ k=2, μ=10 -5 ▼ k=2, μ=10 -6 ○ k=2, μ=10 -7 □ k=3, μ=10 -5 ◇ k=3, μ=10 -6 △ k=3, μ=10 -7
The 𝐿 2 -errors in 𝑢 and 𝑝 are represented graphically in Figs. [START_REF] Brenner | Forty years of the Crouzeix-Raviart element[END_REF] for the different viscosity cases. For each 𝑘 = 1, 2, and 3, the error curves for the different viscosity values are quite similar, showing convergence rates of order 𝑘 + 1, which are typical of the mixed formulation in the Darcy limit corresponding to the parameters (0, 𝛼).

Application to coupled Stokes-Darcy problems

Let us consider a two-dimensional Stokes-Darcy model consisting of a Stokes problem for a fluid in the domain Ω 𝑓 that can flow accross an interface Γ 𝑓 𝑝 = Ω 𝑓 ∩ Ω 𝑝 into a saturated porous medium domain Ω 𝑝 , both regions forming the computational domain Ω = Ω 𝑓 ∪ Ω 𝑝 , as graphically represented in Figure 7.

-∇ • 𝜎(𝑢 𝑓 , 𝑝 𝑓 ) = 𝑓 , ∇ • 𝑢 𝑓 = 0, in Ω 𝑓 , ( 39
)
and the Darcy's law equations

K -1 𝑢 𝑝 + ∇𝑝 𝑝 = 0, ∇ • 𝑢 𝑝 = 𝑔 in Ω 𝑝 . ( 40
)
For the examples to be simulated the viscosity 𝜇 > 0 is constant in Ω 𝑓 and constant permeability tensor K = 𝛾I in Ω 𝑝 , 𝛾 > 0 are adopted.

To complete the problem, interface and boundary conditions need to be imposed. The interface conditions at Γ 𝑓 𝑝 = Ω 𝑓 ∩ Ω 𝑝 are

𝑢 𝑓 • 𝑛 𝑓 + 𝑢 𝑝 • 𝑛 𝑝 = 0 Flux continuity (41) -2𝜇 [𝐷(𝑢 𝑓 )𝑛 𝑓 ] • 𝑛 𝑓 + 𝑝 𝑓 = 𝑝 𝑝 Balance of normal forces (42) -2 [𝐷(𝑢 𝑓 )𝑛 𝑓 ] • 𝜏 𝑓 𝑑 = 𝛼 ℬ𝒥 𝑢 𝑓 • 𝜏 𝑓 𝑝 (BJS) condition, (43)
where 𝜏 𝑓 𝑑 is a unit tangent vector field on Γ 𝑓 𝑝 , 𝑛 𝑓 and 𝑛 𝑝 being the unit normal vectors pointing outwards from Ω 𝑓 and Ω 𝑝 . The parameter 𝛼 ℬ𝒥 > 0 is a friction coefficient (in practice, it is obtained from experimental data [START_REF] Beavers | Boundary conditions at a naturally permeable wall[END_REF]). The boundary conditions are:

𝑢 𝑓 = 𝑢 𝐷 on Γ 𝑓 , ( 44
)
𝑢 𝑝 • 𝑛 𝑝 = 𝜃 𝑁 on Γ 𝑝,𝑁 , 𝑝 𝑝 = 𝑝 𝐷 on Γ 𝑝,𝐷 , Γ 𝑓 = 𝜕Ω 𝑓 ∖ Γ 𝑓 𝑝 and Γ 𝑝,𝑁 ∪ Γ 𝑝,𝐷 = Γ 𝑝 = 𝜕Ω 𝑝 ∖ Γ 𝑓 𝑝 .
When Γ 𝑝𝐷 = ∅, these boundary conditions require the compatibility constraints (𝑔, 1) Ω𝑝 = 0, and (𝑝, 1) Ω = 0. In the applications to be simulated, 𝜃 𝑁 is a piecewise constant function.

Hybrid-mixed FE method for the Stokes-Darcy problem

We propose a discretization of the aforementioned coupled Stokes-Darcy problem by adopting a mixed method for the porous media flow in Ω 𝑑 and a hybrid-mixed model for the approximation of the Stokes equations, following similar steps presented in previous sections for Stokes problems. Precisely,

𝑏(𝑢 ℎ , 𝑞) = -(𝑔, 𝑞) Ω𝑝 , ( 48
) ⟨𝜂 𝑡 , 𝑢 𝑓 ⟩ Γ 𝑓 𝑓 -𝛼 -1 ℬ𝒥 ⟨𝜂 𝑡 , 𝜆 𝑡 ℎ ⟩ Γ 𝑓 𝑝 + ⟨𝜂 𝑡 , 𝑢 𝑓 ⟩ Γ 𝑓 𝑝 = ⟨𝜂 𝑡 , 𝑢 𝐷 ⟩ Γ 𝑓 , (49)

Convergence verification

Consider the Stoke-Darcy model proposed in [START_REF] Chen | Weak Galerkin method for the coupled Darcy-Stokes flow[END_REF] for a free-flow domain Ω 𝑓 = (0, 𝜋) × (0, 1) on top of the porous medium domain Ω 𝑝 = (0, 𝜋) × (-1, 0) with exact solution: The sources f , 𝑔, and the boundary data 𝑢 𝐷 on Γ 𝑓 , and 𝜃 𝑁 on Γ 𝑝,𝑁 = Γ 𝑝 are extracted from the solution. In this section, the new the hybrid-mixed formulation (47)-( 49) is verified for this problem using FE space settings of types 𝒮 𝑓 𝑝 ℎ,𝑅𝑇 (𝑘) (𝜃 𝑁 ) and 𝒮 𝑓 𝑝 ℎ,𝐵𝐷𝑀 + (𝑘) (𝜃 𝑁 ) applied on uniform quadrilateral and triangular meshes 𝒯 ℎ , respectively, with mesh size ℎ = 2 1-𝑁 , 𝑁 = 2 𝑖 , 𝑖 = 2, • • • 6. The resulting 𝐿 2 -errors for 𝑢 = (𝑢 𝑓 , 𝑣 𝑝 ), ∇ • 𝑢 = (∇ • 𝑢 𝑓 , ∇ • 𝑢 𝑝 ), and 𝑝 = (𝑝 𝑓 , 𝑝 𝑝 ) are shown in Tables 6 for each flow regime. Namely, data for the free-flow are listed on the left side, and for the porous media flow they are on the right side.

𝑢 𝑓 = ⎡ ⎣ 𝑣 ′ (𝑦) cos 𝑥 𝑣(𝑦) sin 𝑥 ⎤ ⎦ , 𝑝 𝑓 = sin 𝑥 sin 𝑦, where 𝑣(𝑦) = 1 𝜋 2 sin 2 (𝜋𝑦) -2, ( 50
)
𝑢 𝑝 = ⎡ ⎣ (𝑒 -𝑦 -𝑒 𝑦 ) cos 𝑥 -(𝑒 -𝑦 + 𝑒 𝑦 ) sin 𝑥 ⎤ ⎦ , 𝑝 𝑝 = (-𝑒 -𝑦 + 𝑒 𝑦 ) sin 𝑥. ( 51
)
Optimal convergence rates of order 𝑘 +1 occur for 𝑢 𝑝 , ∇•𝑢 𝑝 and 𝑝 𝑝 in the porous domain Ω 𝑝 , as in usual flux, divergence of the flux, and pressure approximations by the mixed methods for Darcy's flows based on divergence-compatible FE pairs of 𝑅𝑇 (𝑘) or 𝐵𝐷𝑀 + (𝑘) types. In the Stokes domain, convergence for 𝑢 𝑓 and 𝑝 𝑓 is in accordance with the numerical studies for the HM-H(div)(𝒮 ℎ,𝑅𝑇 (𝑘) ) and HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) ) formulations presented in the previous sections for pure Stokes flows, including the divergence-free property for 𝑢 𝑓 .

Application to a carbonate karst reservoir

The numerical study of this section is for a two-dimensional carbonate karst reservoir model composed of multiple vugs with different shapes and sizes in a rectangular domain Ω = (0, 3) × (0, 2), as illustrated in The HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) ) method for 𝑘 = 2 is adopted, based on a triangular mesh created by the finite element mesh generator Gmsh, composed of 11422 nodes and 195665 elements. The resulting numerical results are plotted in Figs. 9 for the velocity magnitude 𝑢(𝑥, 𝑦) and the pressure 𝑝(𝑥, 𝑦) fields. In the velocity field solution, the streamlines show that the vug surfaces affect the flow direction and penetration. Considering this horizontal inlet flow, we observe that velocity magnitudes inside a vug decrease as the vug dimension in the flow direction is larger than the other one.

Conclusions

A numerical scheme for approximating Stokes, Brinkman and/or coupled Stokes-Darcy equations is presented. The core feature of the scheme is the use of divergence-compatible pairs of FE spaces in H(div)-𝐿 2 to approximate the fluid velocity and pressure. The continuity of the tangential velocity is weakly enforced by a Lagrange multiplier variable. In this sense, the method can be seen as a semi-hybrid mixed formulation.

The approximations are pointwise mass conservative. From the theoretical point of view, we were able to demonstrate the solvability of the proposed scheme and to obtain error estimates in terms of optimal erros and consistency errors. We also consider the method for a Stokes-Darcy coupling, where the term modeled by Beaver-Joseph-Saffman is elegantly approximated adding an 𝐿 2 product term in the tangent velocity space.

The focus of the current simulations is on two-dimensional approximations. The results confirm optimal rates of convergence, pressure robustness and a smooth transition for the Darcy limit of the Brinkman problem. When compared to discontinuous Galerkin approximations using the same H(div)-𝐿 2 FE pairs, the proposed scheme gives similar convergence results. Moreover, in the semi-hybrid mixed method a large number of degrees of freedom can be statically condensed resulting in a smaller global system of equations. Efforts to extend the technique to three dimensions are in development.

Figure 1 :

 1 Figure 1: Triangular and rectangular partitions 𝒯 ℎ = {𝐾 0 , 𝐾 1 } used for the rank condition numerical tests.

Figure 2 :

 2 Figure 2: Mesh aspects for computational implementation of hybrid schemes.

Table 4 :

 4 Stokes (𝜇, 𝛼) = (1, 0) and Brinkman (𝜇, 𝛼) = (1, 1) model problems: convergence history for 𝑢, ∇ • 𝑢, and 𝑝, given by the formulations HM-H(div)(𝒮 ℎ,𝑅𝑇 (𝑘) ) (left side) and HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) )(right side), with 𝑚 = 𝑘 -1, and ℎ = 2/𝑁 .

Figure 3 :

 3 Figure 3: Percentage of degrees-of-freedom in the global system at the refinement level 𝑁 = 64 used to solve the Stokes problem by the formulations HM-H(div)(𝒮 ℎ,𝑅𝑇 (𝑘) ) and HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) applied to quadrilateral and triangular meshes.

Figure 4 :

 4 Figure 4: Stokes problem (𝜇, 𝛼) = (1, 0): convergence history for 𝑢 (left side) and 𝑝 (right side), obtained with symmetric DG methods 𝑅𝑇 𝐷𝐺 (3) and 𝑅𝑇 div,𝐷𝐺 (3), and with the HM-H(div)(𝒮 ℎ,𝑅𝑇 (3) formulation based on quadrilateral meshes (top figures); Similar error curves for the symmetric DG methods 𝐵𝐷𝑀 +,𝐷𝐺 (3) and 𝐵𝐷𝑀 +,div,𝐷𝐺 (3), and for the HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (3) ) formulation based on triangular meshes (bottom curves)

Figure 5 :

 5 Figure 5: Pressure-robustness verification: convergence history for 𝑢 (left side) and 𝑝 (right side), obtained with HM-H(div)(𝒮 ℎ,𝑅𝑇 (𝑘) ) method (top) and HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) ) method(bottom), for 𝑘 = 2 and 𝑚 = 1.

Figure 6 :

 6 Figure 6: Brinkman problem (1, 𝜇), 𝜇 = 10 -5 , 10 -6 and 10 -7 : convergence history for 𝑢 (left side) and 𝑝 (right side), obtained with the formulations HM-H(div)(𝒮 ℎ,𝑅𝑇 (𝑘) ) (top side) and HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) )(bottom side), 𝑘 = 1, 2, 3.

Figure 7 :

 7 Figure 7: Coupling problem scheme: fluid domain Ω 𝑓 , porous medium Ω 𝑝 and coupling interface Γ 𝑓 𝑝 .

  ), 𝑝 = (𝑝 𝑓 , 𝑝 𝑝 ), and (∇ • 𝑢 𝑓 , ∇ • 𝑢 𝑝 ) of the solution (50)-(51) for the Stokes-Darcy coupled problem approximated by the HM-H(div)(𝒮 𝑓 𝑝 ℎ,𝐵𝐷𝑀 + (𝑘) (0)) and HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) ^𝑓 𝑝(0) methods based on uniform quadrilateral and triangular meshes 𝒯 ℎ , ℎ = 2 1-𝑁 .

Figure 8 .

 8 The free-flow region Ω 𝑓 is composed by the vugs, which are surrounded by a porous medium Ω 𝑝 . Since 𝜕Ω 𝑓 ∩ 𝜕Ω = ∅, then Γ 𝑆𝐷 = 𝜕Ω 𝑓 . The Neumann boundary datum is: 𝜃 𝑁 = 0 on bottom and top horizontal sides of Ω, inlet 𝜃 𝑁 = 1 on the left vertical side. Dirichlet condition 𝑝 𝑝 = 0 is enforced on the right vertical side. The parameters 𝜇 = 0.1 and 𝛼 ℬ𝒥 = 1 are considered, and the diagonal permeability tensor K = 0.5I is taken in Ω 𝑝 . Vanishing source terms 𝑔 = 0 and f = 0 are applied.

Figure 8 :

 8 Figure 8: Triangular mesh refinement of a vuggy porous media and a zoom in of the marked region.

Figure 9 :

 9 Figure 9: Stokes-Darcy problem with multiple vugs: numerical results for the velocity magnitude field 𝑢(𝑥, 𝑦) and a zoom in view of the marked region (top images); pressure field 𝑝(𝑥, 𝑦) (bottom image).

𝜆 ∈ Λ, the dual norm ||𝜆|| Λ = sup 𝑣∈V ⟨𝜆,𝑣⟩ 𝜕𝒯 ||𝑣|| 1,𝒯 is used. Notice that, when restricted to 𝜆 𝑡 ∈ Λ 𝑡 , this norm can be equivalently stated as ||𝜆 𝑡 || Λ = sup 𝑣∈𝑋 ⟨𝜆 𝑡 ,𝑣⟩ ||𝑣|| 1,𝒯 . The following characterization result for 𝐻 1 0 (Ω, R 𝑑 ) as a subspace of 𝑋 holds as a direct consequence of [20, Lemma 4]. Lemma 1. A continuous

  linear functional 𝐿 on the space 𝑋 vanishes on 𝐻 1 0 (Ω, R 𝑑 ) if and only if there exists a unique element 𝜆 𝑡 ∈ Λ 𝑡 such that 𝐿(𝑣) = ⟨𝜂 𝑡 ,

			1 -norm || • || 1,𝒯 associated
	to the inner product (𝑢, 𝑣) 1,𝒯 = (𝑢, 𝑣) +	∑	𝐾∈𝒯 (∇𝑢, ∇𝑣) 𝐾 .
	As for

𝑣⟩ 𝜕𝒯 , 𝑣 ∈ 𝑋. In

	other
	words, 𝐻 1 0

Table 1 :

 1 Examples of FE pairs X( K) × Ψ( K) in the master element for velocity and pressure fields

Table 2 :

 2 Exact FE sequences for 𝐵𝐷𝑀 (𝑘) in T , 𝑅𝑇 (𝑘) in R, and for their enriched versions 𝐵𝐷𝑀 + (𝑘) and 𝑅𝑇 + (𝑘), respectively.

	𝑘	Rank of C 𝑇	dim Λ 𝑡 ℎ	𝑘	Rank of C 𝑇	)	dim Λ 𝑡
	1	4	5	5	1	7	7		7
	2	9	10	10	2	14	14		14
	3	14	15	15	3	21	21		21

ℎ 𝐵𝐷𝑀 (𝑘) 𝐵𝐷𝑀 + (𝑘) 𝑅𝑇 (𝑘) 𝑅𝑇 + (𝑘

Table 3 :

 3 Rank condition tests

  Quadrilateral elements: 𝒮 ℎ,𝑅𝑇 (𝑘)Triangular elements: 𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) Stokes problem (𝜇, 𝛼) = (1, 0)

	𝑁	||𝑢 -𝑢 ℎ ||	||𝑝 -𝑝 ℎ ||	||∇ • (𝑢 -𝑢 ℎ )|| 𝑘 = 1 ||𝑢 -𝑢 ℎ ||	||𝑝 -𝑝 ℎ ||	||∇ • (𝑢 -𝑢 ℎ )||
	4	1.91𝐸 -02	4.54𝐸 -02	5.57𝐸 -15	5.11𝐸 -02	2.61𝐸 01	1.35𝐸 -14
	8	4.75𝐸 -03	1.61𝐸 -02	1.33𝐸 -14	1.31𝐸 -02	1.35𝐸 -01	3.20𝐸 -14
	16	1.18𝐸 -03	5.68𝐸 -03	2.02𝐸 -14	3.33𝐸 -03	6.83𝐸 -02	6.90𝐸 -14
	32	2.96𝐸 -04	2.01𝐸 -03	4.26𝐸 -14	8.38𝐸 -04	3.43𝐸 -02	1.49𝐸 -13
	64	7.39𝐸 -05	7.10𝐸 -04	1.15𝐸 -14	2.10𝐸 -04	1.72𝐸 -02	1.47𝐸 -14
	Rate	2.00	1.50	-	2.00	1.00	-
	𝑁				𝑘 = 2
	4	1.37𝐸 -03	3.45𝐸 -03	2.93𝐸 -15	3.03𝐸 -03	2.35𝐸 -02	7.10𝐸 -14
	8	1.70𝐸 -04	7.83𝐸 -04	8.76𝐸 -15	3.81𝐸 -04	6.09𝐸 -03	6.45𝐸 -15
	16	2.12𝐸 -05	1.56𝐸 -04	1.49𝐸 -14	4.80𝐸 -05	1.54𝐸 -03	1.41𝐸 -14
	32	2.65𝐸 -06	2.93𝐸 -05	3.05𝐸 -14	6.02𝐸 -06	3.87𝐸 -04	2.83𝐸 -14
	64	3.31𝐸 -07	5.31𝐸 -06	6.03𝐸 -14	7.55𝐸 -07	9.68𝐸 -05	5.54𝐸 -14
	Rate	3.00	2.46	-	3.00	2.00	-
	𝑁				𝑘 = 3
	4	2.83𝐸 -05	1.49𝐸 -04	6.75𝐸 -15	1.61𝐸 -04	1.40𝐸 -03	7.55𝐸 -15
	8	1.77𝐸 -06	1.30𝐸 -05	1.21𝐸 -14	9.64𝐸 -06	1.62𝐸 -04	1.32𝐸 -14
	16	1.10𝐸 -07	1.15𝐸 -06	2.55𝐸 -14	5.88𝐸 -07	1.94𝐸 -05	2.71𝐸 -14
	32	6.87𝐸 -09	1.01𝐸 -07	5.18𝐸 -14	3.63𝐸 -08	2.37𝐸 -06	5.50𝐸 -14
	64	4.29𝐸 -10	8.96𝐸 -09	9.96𝐸 -14	2.26𝐸 -09	2.93𝐸 -07	1.09𝐸 -13
	Rate	4.00	3.50	-	4.01	3.02	-

Brinkman problem (𝜇, 𝛼) = (1, 1)

Table 5 :

 5 Pressure-robustness verification: ||∇ • (𝑢 -𝑢 ℎ )|| for multliplying factors 𝑚 𝑝𝑟 , obtained with the HM-H(div)(𝒮 ℎ,𝑅𝑇 (𝑘) ) and HM-H(div)(𝒮 ℎ,𝐵𝐷𝑀 + (𝑘) ) formulations with 𝑘 = 2, ℎ = 1/𝑁 , and 𝑚 = 1.

	4	𝑚𝑝𝑟 = 10 6

  1. Construction of shape regular and conformal partitions for the flow domains, 𝒯 ℎ,𝑓 = {𝐾 𝑓 } and 𝒯 ℎ,𝑝 = {𝐾 𝑝 }, such that the conglomeration 𝒯 ℎ = 𝒯 ℎ,𝑓 ∪ 𝒯 ℎ,𝑝 is a partition of Ω conformal along the interface Γ 𝑓 𝑝 . 2. Definition of a divergence-compatible FE pair 𝑋 ℎ × Ψ ℎ ⊂ 𝑋 × Ψ based on 𝒯 ℎ . Depending on the boundary conditions, subspaces are considered: (i) 𝑋 ℎ (𝜃 𝑁 ) for fields(𝑣 𝑓 , 𝑣 𝑝 ) ∈ 𝑋 ℎ such that 𝑣 𝑝 • n 𝑝 | Γ 𝑝𝑁 = 𝜃 𝑁 ; (ii) in the case Γ 𝑝𝐷 = ∅, set Ψ ℎ0 = Ψ ℎ ∩ Ψ 0 .3. The Beavers-Joseph-Saffman is weakly enforced using a Lagrange multiplier𝜆 𝑡 ℎ | Γ 𝑓 𝑝 ≈ -[𝜎𝑛 𝑓 ] • 𝜏 𝑓 𝑝 = -2𝜇[𝐷(𝑢)𝑛 𝑓 ] • 𝜏 𝑓 𝑝 = 𝛼 ℬ𝒥 𝑢 𝑓 • 𝜏 𝑓 𝑝 . Moreover,for the hybridization of the Stokes model, let the mesh skeleton Γ 𝑓 𝑓 be formed by the edges of the elements 𝐾 𝑓 ∈ 𝒯 ℎ,𝑓 not included in Γ 𝑓 𝑝 . We consider a FE trace space Λ 𝑡 𝑓,ℎ based on Γ 𝑓 𝑓 ∪ Γ 𝑓 𝑝 verifying the trace compatibility condition with respect to 𝑋 ℎ,𝑓 = 𝑋 ℎ | Ω 𝑓 . For any 𝜂 𝑡 ∈ Λ 𝑡 𝑓,ℎ define ⟨𝜂 𝑡 , 𝑣⟩ Γ 𝑓 𝑓 := ∑ 𝐾 𝑓 ∈𝒯 ℎ,𝑓 ⟨𝜂 𝑡 , 𝑣⟩ 𝜕𝐾 𝑓 ∖Γ 𝑓 𝑝 4. Define the FE setting 𝒮 𝑓 𝑝 ℎ (𝜃 𝑁 ) = 𝑋 ℎ (𝜃 𝑁 )×Ψ ℎ ×Λ 𝑡 𝑓,ℎ (or 𝒮 𝑓 𝑝 ℎ (𝜃 𝑁 ) = 𝑋 ℎ (𝜃 𝑁 )×Ψ ℎ0 ×Λ 𝑡

	𝑓,ℎ
	(46)
	Based on the above setting, we propose a hybrid-mixed FE method for the coupled Stokes-
	Darcy problem (39)-(43) expressed in the form:
	HM-H(div)(𝒮 𝑓 𝑝 ℎ ): find (𝑢

if Γ 𝑝𝐷 = ∅).

Observe that the flux continuity

[START_REF] Beavers | Boundary conditions at a naturally permeable wall[END_REF]

is directly obtained from the H(div)-conforming property of the FE space 𝑋 ℎ . Moreover, the balance of normal forces in

[START_REF] Chen | Weak Galerkin method for the coupled Darcy-Stokes flow[END_REF]

leads us to the definition of the bilinear form in 𝑋 ℎ : 𝑎 𝑓 𝑝 ℎ (𝑢, 𝑣) := 2𝜇 ∑ 𝐾 𝑓 ∈𝒯 ℎ,𝑓 (𝐷(𝑢 𝑓 ), 𝐷(𝑣 𝑓 )) 𝐾 𝑓 + ∑ 𝐾𝑝∈𝒯 ℎ,𝑝 (K -1 𝑢 𝑝 , 𝑣) 𝐾𝑝 . ℎ , 𝑝 ℎ , 𝜆 𝑡 ℎ , 𝜌 ℎ ) ∈ 𝒮 𝑓 𝑝 ℎ (𝜃 𝑁 ) satisfying for all (𝑣, 𝑞, 𝜂 𝑡 ) ∈ 𝒮 𝑓 𝑝 ℎ (0): 𝑎 𝑓 𝑝 ℎ (𝑢 ℎ , 𝑣) + 𝑏(𝑣, 𝑝 ℎ ) + ⟨𝜆 𝑡 ℎ , 𝑣⟩ Γ 𝑓 𝑓 + ⟨𝜆 𝑡 ℎ , 𝑣⟩ Γ 𝑓 𝑝 = ⟨𝑣 • 𝑛 𝑝 , 𝑝 𝐷 ⟩ Γ 𝑝,𝐷 + (𝑓 , 𝑣) Ω 𝑓 , (47)

Table 6 :

 6 Quadrilateral elements:𝒮 𝑓 𝑝 ℎ,𝑅𝑇 (𝑘) (𝜃 𝑁 ) 𝑁 ||𝑢 𝑓 -𝑢 𝑓 ℎ || ||𝑝 𝑓 -𝑝 𝑓 ℎ || ||∇ • (𝑢 𝑓 -𝑢 𝑓 ℎ )|| ||𝑢𝑝 -𝑢𝑝 ℎ || ||𝑝𝑝 -𝑝𝑝 ℎ || ||∇ • (𝑢𝑝 -𝑢𝑝 ℎ )|| 𝐿 2 -errors in 𝑢 = (𝑢 𝑓 , 𝑣 𝑝

	𝑘 = 1

+ (𝑘) (𝜃 𝑁 ) 𝑁 ||𝑢 𝑓 -𝑢 𝑓 ℎ || ||𝑝 𝑓 -𝑝 𝑓 ℎ || ||∇ • (𝑢 𝑓 -𝑢 𝑓 ℎ )|| ||𝑢𝑝 -𝑢𝑝 ℎ || ||𝑝𝑝 -𝑝𝑝 ℎ || ||∇ • (𝑢𝑝 -𝑢𝑝 ℎ )||

NeoPZ open-source platform: http://github.com/labmec/neopz

Acknowledgments

The authors thankfully acknowledge financial support from CAPES -Brazilian Federal Agency for Post-graduate Education -(grant 01P-4376/2015), CNPq -the Brazilian Research Council (grants 305425/2013-7, 304029/2013-0, 305823-2017-5, and 306167/2017-4), and FAPESP -Research Foundation of the State of São Paulo, Brazil (grants 2016/05155, and 2016/13637-5).