# Multiscale hybrid-mixed methods for the stokes and brinkman equations - a priori analysis

Araya Rodolfo, Harder Christopher, Poza Abner, Frederic Valentin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Araya Rodolfo, Harder Christopher, Poza Abner, Frederic Valentin. Multiscale hybrid-mixed methods for the stokes and brinkman equations - a priori analysis. 2024. hal-04425367v2

HAL Id: hal-04425367

https://hal.science/hal-04425367v2

Preprint submitted on 19 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MULTISCALE HYBRID-MIXED METHODS FOR THE STOKES AND BRINKMAN EQUATIONS – A PRIORI ANALYSIS

RODOLFO ARAYA, CHRISTOPHER HARDER, ABNER H. POZA, AND FRÉDÉRIC VALENTIN

**ABSTRACT.** The multiscale hybrid-mixed (MHM) method for the Stokes operator was formally introduced in Araya at al. (2017) and numerically validated. The method has face degrees of freedom associated with multiscale basis functions computed from local Neumann problems driven by discontinuous polynomial spaces on skeletal meshes. The two-level MHM version approximates the multiscale basis using a stabilized finite element method. This work proposes the first numerical analysis for the one- and two-level MHM method applied to the Stokes/Brinkman equations within a new abstract framework relating MHM methods to discrete primal hybrid formulations. As a result, we generalize the two-level MHM method to include general second-level solvers and continuous polynomial interpolation on faces and establish abstract conditions to have those methods well-posed and optimally convergent on natural norms. We apply the abstract setting to analyze the MHM methods using stabilized and stable finite element methods as second-level solvers with (dis)continuous interpolation on faces. Also, we find that the discrete velocity and pressure variables preserve the balance of forces and conservation of mass at the element level. Numerical benchmarks assess theoretical results.

## 1. INTRODUCTION

Let $\Omega \subset \mathbb{R}^d$, $d \in \{2, 3\}$, be an open, bounded polyhedron domain with Lipschitz boundary $\partial\Omega$. We consider the generalized Stokes problem, also called Brinkman model, which corresponds to finding the velocity $\mathbf{u}$ and the pressure $p$ such that

$$ 
\begin{align} 
-\nu \Delta \mathbf{u} + \boldsymbol{\theta} \mathbf{u} + \nabla p &= \mathbf{f} \quad \text{in } \Omega, \\ 
(1.1) \qquad \nabla \cdot \mathbf{u} &= 0 \quad \text{in } \Omega, \\ 
\mathbf{u} &= \mathbf{g} \quad \text{on } \partial\Omega, 
\end{align} 
$$

where $\mathbf{f} \in L^2(\Omega)^d$ and $\mathbf{g} \in H^{1/2}(\partial\Omega)^d$ with $\int_{\partial\Omega} \mathbf{g} \cdot \mathbf{n} \, d\mathbf{s} = 0$, and $\mathbf{n}$ is the outward normal vector on $\partial\Omega$. The viscosity $\nu$ is a positive constant, and the reaction coefficient $\boldsymbol{\theta} = \boldsymbol{\theta}(\mathbf{x})$ is a
semidefinite positive, symmetric tensor which is uniformly elliptic, i.e., there exist constants $c_{min} \ge 0$ and $c_{max} > 0$ such that

$$ (1.2) \qquad 0 \le c_{min} |\xi|^2 \le \theta_{min}(\mathbf{x}) |\xi|^2 \le \xi^T \cdot \theta(\mathbf{x}) \xi \le \theta_{max}(\mathbf{x}) |\xi|^2 \le c_{max} |\xi|^2, $$

for all $\xi \in \mathbb{R}^d$ and $\mathbf{x}$ a.e. in $\Omega$, where $\theta_{min}$ and $\theta_{max}$ are the smallest and largest eigenvalues of $\theta$ ($c_{min} > 0$ if $\theta$ is a definite positive matrix). Also, $\theta$ may contain multiscale geometrical features of the media. We recognize the Stokes problem in (1.1) if $\theta = 0$. The Dirichlet boundary condition in (1.1) is chosen for the sake of the presentation, Neumann or Robin boundary conditions can be easily accommodated in what follows. The Stokes-Brinkman's unique weak solution $(\mathbf{u}, p) \in H^1(\Omega)^d \times L_0^2(\Omega)$, with $\mathbf{u}|_{\partial\Omega} = \mathbf{g}$, satisfies

$$ (1.3) \qquad \int_{\Omega} (\nu \nabla \mathbf{u} : \nabla \mathbf{v} + \theta \mathbf{u} \cdot \mathbf{v}) \, d\mathbf{x} - \int_{\Omega} p \nabla \cdot \mathbf{v} \, d\mathbf{x} = \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \, d\mathbf{x} \quad \text{for all } \mathbf{v} \in H_0^1(\Omega)^d, \\ \qquad \int_{\Omega} q \nabla \cdot \mathbf{u} \, d\mathbf{x} = 0 \quad \text{for all } q \in L_0^2(\Omega), $$

and there exists positive C, dependent on $\nu$, $\theta$ and $\Omega$, such that

$$ \|\mathbf{u}\|_{1,\Omega} + \|p\|_{0,\Omega} \le C (\|\mathbf{f}\|_{0,\Omega} + \|\mathbf{g}\|_{1/2,\partial\Omega}). $$

Given any measurable set $D \subset \bar{\Omega}$ and any integer $m \ge 0$, we respectively denote by $| \cdot |_{m,D}$ and $\| \cdot \|_{m,D}$ the standard semi-norm and norm in $H^m(D)$. The $L_0^2(D)$ space stands for the functions in $L^2(D) := H^0(D)$ with zero mean value in $D$. We use the convention of denoting (matrix-)vector-valued functions and spaces in bold. Hereafter, we shall denote by C a positive constant, which is independent of any partition (mesh) parameter but can change in each occurrence.

When attempting to approximate solutions to problem (1.3), one must deal with some (potential) numerical instabilities. First, one must either choose compatible approximation spaces (in the sense of the inf-sup condition [20]) or use stabilized numerical methods that overcome this incompatibility [54]. Furthermore, in the case of a dominant reaction term, numerical methods must be robust with respect to vanishing diffusion coefficient due to the singularly perturbed nature of the solutions characterized by boundary layers. Without such care, spurious non-physical oscillations may plague the numerical approximation. Beyond these, solutions may exhibit poor approximation properties arising from the insufficient resolution of the physics occurring at multiple scales coming through heterogeneous and or high contrast coefficients. Finally, it is desirable that numerical methods for the Stokes/Brinkman problem produce discrete velocity fields respecting the local mass balance.
Theoretically, the issue involving multiple scales and boundary layers may be overcome by choosing a mesh whose size is smaller than the smallest relevant scale for the physics. Practically, this can lead to problem statements that are intractable computationally. For this reason, multi-scale methods are proposed to correctly incorporate the effect of subscales (including boundary layers) into a numerical scheme proposed on a relatively coarse partition. Such methods use a form of upscaling to incorporate sub-scales within the scheme on the coarse mesh. This idea underlies the Variational Multiscale Method (VMS) [53], Residual-Free Bubbles (RFB) [63] and other enriched finite element methods as the Residual Local Projection (RELP) method [17] and the Petrov-Galerkin Enriched (PGEM) method [4]. See [34] for an overview. Some of these multiscale methods are directly related to stabilized finite element methods, as pointed out in [23, 24, 16], for example, and they deal with issues related to the inf-sup condition and the approximation of boundary layers at the same time.

The idea of multiscale finite element methods goes back to [13] where the multiscale basis concept was introduced and analyzed first in one dimension for a highly oscillatory coefficient problem. This seminal work was later extended to higher dimensions in the form of the MsFEM [52]. Following such an idea, other multiscale methods have been proposed for various operators, including the heterogeneous multiscale method [37], localizable orthogonal decomposition [58], the sub-grid finite element method [10], and the generalized multiscale finite element method [38]. More recently, multiscale methods have been devised for Stokes and Brinkman problems [1, 5, 46], with discontinuous methods which are proving effective in highly heterogeneous media context [25, 55].

The Multiscale Hybrid-Mixed (MHM) finite element method fits with relatively recent attempts to use hybridization to resolve multiple scales [11, 39, 65]. First analyzed in the context of the Poisson problem [61], hybrid methods for various operators are intrinsically related to the domain decomposition methods [3, 26, 27], the discontinuous enriched method (DEM) [42], and strategies to reduce computational cost related to saddle point problems [12, 33]. As for the MHM method, it was proposed for the first time for the Poisson (Darcy) equation in [49] with an a priori and a posteriori error analysis developed in [6, 15, 59, 60]. The methodology was extended to other operators, such as the linear elasticity model in [48, 47], the complete transport equation in [50, 8], and wave problems in [56, 31]. Furthermore, the MHM method is closely related to other multiscale methods, such as the MsHHO [29] and the multiscale version of the Crouzeix-Raviart element [57], and retrieves the lowest-order Raviart-Thomas element in the homogeneous coefficient case [49]. In its two-level version, it is related to the FETI domain decomposition approach [43] and can also be seen as the
dual version of the MsFEM and the multiscale mortar methods [11] (see also [21] for a recent primal version of the multiscale mortar method).

In a companion paper [7] to the present work, the MHM method was extended to the Stokes/Brinkman problem (1.1), and a residual a posteriori error estimator was proposed in [9]. The philosophy underlying the MHM method involves using a hybrid form of (1.1) posed on a coarse partition. The approach decomposes the exact solution into local and global coupled systems, and discretization decouples the system as follows: the global formulation is responsible for the degrees of freedom over the skeleton of the coarse partition, and the local problems driven by upscaling operators $T$ and $\hat{T}$ provide the multiscale basis functions resolving sub-scales. To fix ideas, consider $\mathcal{P}_H$ a general partition of a domain $\Omega$ with characteristic length $\mathcal{H}$ decomposed into elements $K$ with contour $\partial K$. The skeleton of $\mathcal{P}_H$, denoted by $\partial \mathcal{P}_H$, is partitioned with elements of diameter $H$. Remember (c.f. [7]) that the one-level MHM method for the Stokes equation consists of finding a flux $\lambda_H$ in a polynomial space $\Lambda_H$ defined on the skeleton of $\mathcal{P}_H$, a velocity $u_0^H$ in $V_0$, the piecewise constant space over $\mathcal{P}_H$, and $\rho \in \mathbb{R}$ such that

$$ (1.4) \qquad \sum_{K \in \mathcal{P}_H} \int_{\partial K} \lambda_H \cdot v_0 \, ds = \sum_{K \in \mathcal{P}_H} \int_K f \cdot v_0 \, dx \quad \text{for all } v_0 \in V_0, \\ \qquad \sum_{K \in \mathcal{P}_H} \int_{\partial K} \mu_H \cdot u_H \, ds = \int_{\partial \Omega} \mu_H \cdot g \, ds \quad \text{for all } \mu_H \in \Lambda_H, \\ \qquad \xi \int_{\Omega} p_H \, dx = 0 \quad \text{for all } \xi \in \mathbb{R}, $$

where the discrete velocity and pressure variables $(u_H, p_H)$ in (1.4) depend on the coarse scale variables $(u_0^H, \lambda_H, \rho)$ in the form

$$ (1.5) \qquad (u_H, p_H) := (u_0^H, 0) + T(\lambda_H, \rho) + \hat{T}(f). $$

The fine scales are resolved via the locally defined bounded mappings $T$ and $\hat{T}$, which have finite-dimensional images in the local product space $[H^1(K) \cap L_0^2(K)]^d \times L^2(K)$. Specifically, the mappings $T$ and $\hat{T}$ are inverses of well-posed Stokes problems with prescribed Neuman boundary conditions on each $K \in \mathcal{P}_H$. As a result, the MHM method (1.4) is non-conforming in $H^1(\Omega)$ as the discrete velocity in (1.5) does not belong to the $H_0^1(\Omega)$ space. Note that the second equation in (1.4) reinforces the weak continuity of the velocity field in the boundary elements, while the third equation guarantees the uniqueness of the pressure. The first equation in (1.4) imposes that the numerical flux is in local equilibrium with external force, which is a consequence of the spatial decomposition adopted in the construction of the MHM method (see [7] and Section 3 for more details). Furthermore, the
discrete velocity field $u_H$ is locally divergent free, i.e.,

$$ \nabla \cdot u_H = 0 \quad \text{in all } K \in \mathcal{P}_H. $$

In [7], the MHM method was presented along with numerical results involving a stabilized finite element method to approximate local problems, but numerical analysis was absent. The use of these local numerical schemes is viewed as forming approximations of multiscale basis functions and leads to two-level methods. In the context of [7], the choice of the unusual stabilized finite element method (USFEM) [16] intended to take advantage of its well-known robustness to approximate boundary layers that can be presented in the multiscale basis while making the pair of polynomials spaces of the same order available to approximate both pressure and velocity multiscale basis functions. The efficiency of this choice has been extensively verified in [7] numerically.

In the present paper, we fill this theoretical gap and propose the first numerical analysis for the one- and two-level MHM methods applied to the Stokes/Brinkman equations. For example, we prove that (1.4) is well-posed and the error due to exact flow approximation $\lambda_H$ of $\lambda := (\nu \nabla u - p \mathbb{I})n^K$ on $\partial K$ turns out to be an upper bound for the error associated with velocity and pressure approximations, where $n^K$ is the outward normal vector on $\partial K$ for all $K \in \mathcal{P}_H$ and $\mathbb{I}$ is the $d \times d$ identity operator. Notably, we prove (see Theorem 3.2 in the context of the one-level MHM method)

$$ \|u - u_H\|_{1,\mathcal{P}_H} + \|p - p_H\|_{0,\Omega} \le C \inf_{\substack{\mu_H \in \Lambda_H \\ (\lambda - \mu_H, v_0)_{\partial \mathcal{P}_H} = 0 \quad \forall v_0 \in \mathbf{V}_0}} \|\lambda - \mu_H\|_{\Lambda}, $$

where $\| \cdot \|_{\Lambda}$ is a norm in $\Lambda$ and $(\cdot, \cdot)_{\partial \mathcal{P}_H}$ means a dual product on the skeleton $\partial \mathcal{P}_H$ (see Section 2.2 for precisions). So, convergence arises by bringing approximability properties to the definition of the finite-dimensional space $\Lambda_H \subset \Lambda$. Furthermore, under local regularity assumptions for the exact velocity and pressure variables, we prove that the two-level MHM solution in [7] converges in the energy norm with rate $O(H^{\ell+1} + h^k)$, where $\ell \ge 0$ is the polynomial degree of interpolation on faces used in $\Lambda_H$, and $h$ is the local submesh diameter and $k \ge 1$ is the polynomial degree of interpolation on sub-meshes (see Theorem 4.2). This way, convergence can be achieved by keeping the macroelement mesh fixed (i.e., $\mathcal{H}$ remains fixed). This property, in addition to its practical interest in avoiding remeshing complex geometries to improve accuracy, leads to superconvergence with additional convergence order $O(H^{1/2})$ under the smoothest local exact solution assumption (see Theorem 4.2) and assuming that the local multiscale basis functions are “sufficiently” approximate. This is also verified numerically in Section 5.
In a broad sense, the proof strategy is as follows: we first establish a relationship between the primal-hybrid version of (1.3) and the MHM method in an abstract general sense, showing that the well-posedness of both are related at the fully discrete level and also establish how their stability constants are related. Then, the numerical analysis of the MHM method proposed in [7] arises from fulfilling the conditions to have the primal-hybrid version of (1.3) well-posed. Thereby, the overall analysis in this work follows a totally different path from that used in previous MHM works, which is fundamental to prove the existence and uniqueness and the best approximation property of the two-level MHM method with the USFEM method used to approximate $T$ and $\hat{T}$ in (1.5). In a sense, such a perspective is also presented in previous works on MHM-type methods applied to the Poisson problem [51, 29, 14].

The correspondence between the solution of the discrete primal-hybrid methods and the two-level MHM methods also produces an interesting reinterpretation of the two-level MHM methods. When the USFEM is used as a second-level solver on one-element sub-meshes, the corresponding discrete primal-hybrid method can be seen as a *new member* of the class of non-conforming stabilized finite element methods. Stabilization then acts at the element level, which has a global impact by weakly imposing continuity on the boundary elements (see [28, 2] for other examples of non-conforming stabilized methods). As a result, local mass is conserved, which is a property missing from conforming-stabilized methods in general (see [17] for a post-processing strategy to recover this property in the conforming scenario). Furthermore, the use of refined sub-meshes in this context can be interpreted as a *multiscale version* of non-conforming stabilized methods. A similar interpretation arises if we replace the stabilized second-level solver with the Galerkin method based on stable pairs of spaces. For example, when adopting the lowest-degree Taylor-Hood element into the scope of the two-level MHM method with single-element sub-meshes, this corresponds to nothing more than the non-conforming Crouzeix-Raviart method [36]. Thus, with the adoption of refined sub-meshes in the context of the stable two-level MHM method, a non-conforming multiscale version of the Crouzeix-Raviart method emerges.

In addition to providing a numerical analysis for the two-level MHM method introduced in [7], in this work we

(i) establish abstract conditions for the well-posedness and optimality of MHM methods for the Stokes/Brinkman model in natural norms. Such conditions allow extending the analysis to more general second-level solvers than the stabilized method used in [7];
(ii) extend the two-level MHM analysis to include a stable second-level solver and continuous polynomial interpolation for the flux variable. This version of the two-level MHM method is new, and we prove that the method also achieves optimal convergence using the abstract conditions of item (i). This version of the two-level MHM method produces point-wise divergence-free exact discrete velocity when it adopts one-element sub-meshes;

(iii) validate numerically the MHM method using continuous interpolation on faces. The results indicate that interpolating the flux continuously improves convergence when compared to the discontinuous case, at least for regular exact solutions. The use of continuous interpolation within the MHM methodology is new, even in the context of other operators.

We also provide some numerical verification about the dependence of the constant in the error estimates in terms of the physical coefficients, but we leave this theoretical question outside the scope of this work.

## 2. HYBRIDIZATION

The MHM methods are built on reformulating a hybridized version of (1.3). In a broad sense, given spaces $V$, $Q$, and $\Lambda$ and (bi)linear forms $f: V \times Q \to \mathbb{R}$, $g: \Lambda \times \mathbb{R} \to \mathbb{R}$, $a: (V \times Q) \times (V \times Q) \to \mathbb{R}$, and $b: (\Lambda \times \mathbb{R}) \times (V \times Q) \to \mathbb{R}$, the hybrid formulation takes the form: Find $(u, p) \in V \times Q$ and $(\lambda, \rho) \in \Lambda \times \mathbb{R}$ such that

$$ (2.1) \qquad \begin{aligned} a(\boldsymbol{u}, p; \boldsymbol{v}, q) + b(\boldsymbol{\lambda}, \rho; \boldsymbol{v}, q) &= f(\boldsymbol{v}, q), \\ b(\boldsymbol{\mu}, \xi; \boldsymbol{u}, p) &= g(\boldsymbol{\mu}, \xi), \end{aligned} $$

for all $(\boldsymbol{v}, q) \in V \times Q$ and $(\boldsymbol{\mu}, \xi) \in \Lambda \times \mathbb{R}$. The hybrid forms in (2.1) arise working in the spaces $V \times Q$ which contain $H^1(\Omega)^d \times L_0^2(\Omega)$. The burden of enforcing $(\boldsymbol{u}, p) \in H^1(\Omega)^d \times L_0^2(\Omega)$ is then borne by the action of the bilinear form $b(\cdot; \cdot)$ defined over appropriate spaces. The remainder of this section is dedicated to defining this setting precisely and establishing the equivalence of (1.3) and (2.1). This will provide a solid basis for defining the MHM formulation of (1.3) in Section 3.

### 2.1. Partitions.
Consider a family of partitions $\{\mathcal{P}_\mathcal{H}\}_{\mathcal{H}>0}$ of $\Omega$ parameterized by $\mathcal{H}:= \max_{K \in \mathcal{P}_\mathcal{H}} h_K$, where $h_K$ is the diameter of simplex elements $K$. The collection of all element boundaries $\partial K$ is denoted $\partial \mathcal{P}_\mathcal{H}$. Without loss of generality, we shall use hereafter the terminology employed for three-dimensional domains. The collection of all faces $E$ in the triangulations, with diameter $h_E$, is denoted $\mathcal{E}$. This set is decomposed into the set of faces
on $\partial\Omega$ denoted $\mathcal{E}_{\partial}$, and its complement $\mathcal{E}_0$. To each $E \in \mathcal{E}$, a normal $n$ is associated, taking care to ensure this is directed outward on $\partial\Omega$. For each $K \in \mathcal{P}_H$ we collect local faces of $E \subset \partial K$ in the set $\mathcal{E}^K$, while denoting the outward normal on $\partial K$ by $n^K$ and defining $n_E^K := n^K|_E$ for each $E \subset \mathcal{E}^K$. Also, we denote by $\{\mathcal{E}_H\}_{H>0}$ a regular family of simplicial partitions of $\mathcal{E}$, where $H := \max_{F \in \mathcal{E}_H} h_F$ and $h_F$ is the diameter of $F \in \mathcal{E}_H$. We collect the faces of $F \subset \partial K$ in the set $\mathcal{E}_H^K$.

2.2. **Broken spaces and norms.** Given a partition $\mathcal{P}_H$ in $\{\mathcal{P}_H\}_{H>0}$, we adopt the notation

$$ (2.2) \quad Q := L^2(\Omega) \quad \text{and} \quad \mathbf{V} := \{\mathbf{v} \in L^2(\Omega)^d : \mathbf{v}_K \in H^1(K)^d \ \forall K \in \mathcal{P}_H\}, $$

where we denote $\mathbf{v}_D := \mathbf{v}|_D$ with $D \subset \bar{\Omega}$ a measurable set, and

$$ (2.3) \quad \mathbf{\Lambda} := \{\boldsymbol{\sigma}_K \mathbf{n}^K|_{\partial K} \text{ for all } K \in \mathcal{P}_H : \boldsymbol{\sigma} \in \mathbf{H}(\text{div}; \Omega)\}. $$

Also, we define the broken gradient operator $\nabla_H: \mathbf{V} \to L^2(\Omega)^{d \times d}$ as such $(\nabla_H \mathbf{v})|_K := \nabla \mathbf{v}_K$ for all $\mathbf{v} \in \mathbf{V}$ and $K \in \mathcal{P}_H$. Given $\mathbf{w}, \mathbf{v} \in L^2(\Omega)^d$, and owing to the notation $(\mathbf{w}, \mathbf{v})_{\mathcal{P}_H} := \sum_{K \in \mathcal{P}_H} (\mathbf{w}_K, \mathbf{v}_K)_K$, where $(\cdot, \cdot)_D$ stands for the $L^2(D)^d$ inner product, we equip $\mathbf{V}$ with the norm $\| \cdot \|_{\mathbf{V}}$ induced by the inner product

$$ (2.4) \quad (\mathbf{w}, \mathbf{v})_{\mathbf{V}} := d_{\Omega}^{-2} (\mathbf{w}, \mathbf{v})_{\mathcal{P}_H} + (\nabla_{\mathcal{P}_H} \mathbf{w}, \nabla_{\mathcal{P}_H} \mathbf{v})_{\mathcal{P}_H}, $$

where $d_{\Omega}$ is the diameter of $\Omega$. Further, we equip $Q$ with the induced norm $\| \cdot \|_Q^2 := (\cdot, \cdot)_Q := (\cdot, \cdot)_{\mathcal{P}_H}$. We denote $(\cdot, \cdot)_{\mathbf{V} \times Q} := (\cdot, \cdot)_{\mathbf{V}} + (\cdot, \cdot)_Q$ the inner-product on the product space $\mathbf{V} \times Q$, which induces the following norm in $\mathbf{V} \times Q$

$$ (2.5) \quad \|\mathbf{v}, q\|_{\mathbf{V} \times Q}^2 := \|\mathbf{v}\|_{\mathbf{V}}^2 + \|q\|_Q^2 \quad \text{for all } (\mathbf{v}, q) \in \mathbf{V} \times Q. $$

We equip the space $\mathbf{H}(\text{div}; \Omega)$ and $\mathbf{\Lambda}$ with the norms,

$$ (2.6) \quad \|\boldsymbol{\sigma}\|_{\mathbf{H}(\text{div};\Omega)}^2 := \sum_{K \in \mathcal{P}_H} (\|\boldsymbol{\sigma}\|_{0,K}^2 + d_{\Omega}^2 \|\nabla \cdot \boldsymbol{\sigma}\|_{0,K}^2), $$

$$ (2.7) \quad \|\boldsymbol{\mu}\|_{\mathbf{\Lambda}} := \inf_{\substack{\boldsymbol{\sigma} \in \mathbf{H}(\text{div};\Omega) \\ \boldsymbol{\sigma}_{\partial K} \mathbf{n}^K|_{\partial K} = \boldsymbol{\mu}_{\partial K}, K \in \mathcal{P}_H}} \|\boldsymbol{\sigma}\|_{\mathbf{H}(\text{div};\Omega)}. $$

Also, we define the following product norm

$$ (2.8) \quad \|\boldsymbol{\mu}, \boldsymbol{\xi}\|_{\mathbf{\Lambda} \times Q}^2 := \|\boldsymbol{\mu}\|_{\mathbf{\Lambda}}^2 + \|\boldsymbol{\xi}\|_Q^2 \quad \text{for all } (\boldsymbol{\mu}, \boldsymbol{\xi}) \in \mathbf{\Lambda} \times \mathbb{R}. $$

The duality pairing between $H^{-1/2}(\partial K)^d$ and $H^{1/2}(\partial K)^d$ is denoted by $\langle \cdot, \cdot \rangle_{\partial K}$, and we define $(\boldsymbol{\mu}, \boldsymbol{\zeta})_{\partial \mathcal{P}_H} := \sum_{K \in \mathcal{P}_H} \langle \boldsymbol{\mu}_{\partial K}, \boldsymbol{\zeta}_{\partial K} \rangle_{\partial K}$ for $\boldsymbol{\mu} \in \prod_{K \in \mathcal{P}_H} H^{-1/2}(\partial K)^d$ and $\boldsymbol{\zeta} \in \prod_{K \in \mathcal{P}_H} H^{1/2}(\partial K)^d$.
Observe that if $\mu \in \Lambda$ and $\zeta_{\partial K} = v_K|_{\partial K}$ for all $K \in \mathcal{P}_H$ where $v \in H_0^1(\Omega)^d$, then $(\mu, \zeta)_{\partial \mathcal{P}_H} = 0$ (c.f. [7, Lemma 4]).

2.3. (Bi)linear Forms. Given, $w, v \in V$ and $q, r \in Q$, and $(\mu, \xi) \in \Lambda \times \mathbb{R}$, let

$$a(w, r; v, q) := \sum_{K \in \mathcal{P}_H} a_K(w_K, r_K; v_K, q_K) \quad \text{and} \quad b(\mu, \xi; v, q) := \sum_{K \in \mathcal{P}_H} b_K(\mu_{\partial K}, \xi_K; v_{\partial K}, q_K),$$

where

$$a_K(v_K, q_K; w_K, r_K) := (\nu \nabla v_K, \nabla w_K)_K + (\theta v_K, w_K)_K - (q_K, \nabla \cdot w_K)_K \\ (2.9) \qquad \qquad \qquad \qquad \qquad \qquad + (\nabla \cdot v_K, r_K)_K,$$

$$b_K(\mu_{\partial K}, \xi_K; v_{\partial K}, q_K) := \langle \mu_{\partial K}, v_{\partial K} \rangle_{\partial K} + (\xi_K, q_K)_K,$$

and

$$(2.10) \qquad f(v, q) := (f, v)_{\mathcal{P}_H} \quad \text{and} \quad g(\mu, \xi) := (\mu, g)_{\partial \Omega}.$$

Observe that since $f \in L^2(\Omega)^d$ and $g \in H^{1/2}(\partial \Omega)^d$, there exist constants $C > 0$ independent of the mesh $\mathcal{P}_H$, such that $f(v, q) = (f, v)_{\mathcal{P}_H} \le C\|v, q\|_{V \times Q}$ for all $(v, q) \in V \times Q$ and $g(\mu, \xi) = (\mu, g)_{\partial \Omega} \le C\|\mu, \xi\|_{\Lambda \times Q}$ for all $(\mu, \xi) \in \Lambda \times \mathbb{R}$. We denote by $\|f\|$ and $\|g\|$ the smallest possible of such constants. Similarly, $a(\cdot; \cdot)$ is uniformly bounded over all $(v, q)$, $(w, r) \in V \times Q$, and then

$$(2.11) \qquad \|a\| := \sup_{(v,q) \in V \times Q} \sup_{(w,r) \in V \times Q} \frac{a(v,q; w,r)}{\|v,q\|_{V \times Q} \|w,r\|_{V \times Q}} < \infty.$$

Above and hereafter, we lighten notation and understand the supremum to be taken over sets excluding the zero function, even though this is not specifically indicated. The nullspace of $a(\cdot; \cdot)$ will play a critical role in the analysis. We define

$$\mathcal{N}_a := \{(v, q) \in V \times Q : (v_K, q_K) \in \mathcal{N}_a(K) \quad \forall K \in \mathcal{P}_H\},$$

where

$$\mathcal{N}_a(K) := \{(v, q) \in H^1(K)^d \times L^2(K) : a_K(v, q; w, r) = 0 \quad \forall (w, r) \in H^1(K)^d \times L^2(K)\},$$

and their orthogonal complements are

$$(2.12) \qquad \mathcal{N}_a(K)^\perp := \{(v, q) \in H^1(K)^d \times L^2(K) : (v, q; w, r)_{V \times Q} = 0 \quad \forall (w, r) \in \mathcal{N}_a(K)\}, \\ \mathcal{N}_a^\perp := \{(v, q) \in V \times Q : (v_K, q_K) \in \mathcal{N}_a(K)^\perp \quad \forall K \in \mathcal{P}_H\}.$$

We notice that, based on the space $\mathcal{N}_a(K)$, the following characterisations hold

$$\text{Stokes model: } \mathcal{N}_a(K) = \mathbb{P}_0(K)^d \times \{0\} \text{ and } \mathcal{N}_a^\perp(K) = [H^1(K) \cap L_0^2(K)]^d \times L^2(K);$$
**Brinkman model:** $\mathcal{N}_a(K)$ is trivial and $\mathcal{N}_a^{\perp}(K) = H^1(K)^d \times L^2(K)$, where $P_0(K)^d$ is the constant vector function space in $K \in \mathcal{P}_H$.

For simplicity, we assume hereafter that $\mathcal{N}_a(K)$ is of Stokes or Brinkman type for all $K \in \mathcal{P}_H$.

**2.4. The well-posedness of hybrid problem (2.1).** With the previous definitions and the definition of the null space $\mathcal{N}_b$ given by

$$ (2.13) \qquad \mathcal{N}_b := \{(v, q) \in \mathbf{V} \times \mathcal{Q} : b(\boldsymbol{\mu}, \xi; v, q) = 0 \text{ for all } (\boldsymbol{\mu}, \xi) \in \mathbf{\Lambda} \times \mathbb{R}\}, $$

we establish that hybrid formulation (2.1) of (1.1) is well-posed.

**Lemma 2.1.** There exists a positive constant $\alpha_b$, independent of $H$, such that

$$ (2.14) \qquad \alpha_b \|v, q\|_{\mathbf{V} \times \mathcal{Q}} \le \sup_{(w,r) \in \mathcal{N}_b} \frac{a(v, q; w, r)}{\|w, r\|_{\mathbf{V} \times \mathcal{Q}}} \quad \text{for all } (v, q) \in \mathcal{N}_b. $$

Moreover, it holds

$$ (2.15) \qquad \|\boldsymbol{\mu}, \xi\|_{\mathbf{\Lambda} \times \mathcal{Q}} = \sup_{(w,r) \in \mathbf{V} \times \mathcal{Q}} \frac{b(\boldsymbol{\mu}, \xi; w, r)}{\|w, r\|_{\mathbf{V} \times \mathcal{Q}}} \quad \text{for all } (\boldsymbol{\mu}, \xi) \in \mathbf{\Lambda} \times \mathbb{R}. $$

Therefore, hybrid formulation (2.1) is well-posed and

$$ (2.16) \qquad \|\boldsymbol{u}, p\|_{\mathbf{V} \times \mathcal{Q}} \le \frac{1}{\alpha_b} \|f\| + \left(1 + \frac{\|a\|}{\alpha_b}\right) \|g\|, \quad \|\boldsymbol{\lambda}, \rho\|_{\mathbf{\Lambda} \times \mathcal{Q}} \le \left(1 + \frac{\|a\|}{\alpha_b}\right) \left(\|f\| + \|a\|\|g\|\right). $$

**Proof.** First observe that $\mathcal{N}_b = H_0^1(\Omega)^d \times L_0^2(\Omega)$ (c.f. [7, Lemma 4]), and then the first result in (2.14) stems from the classical well-posedness of (1.3) over $H_0^1(\Omega)^d \times L_0^2(\Omega)$. Next, let $(\boldsymbol{\mu}, \xi) \in \mathbf{\Lambda} \times \mathbb{R}$. From Green's Theorem and the Cauchy-Schwarz inequality, we get

$$ b(\boldsymbol{\mu}, \xi; v, q) = \sum_{K \in \mathcal{P}_H} \left[ (\boldsymbol{\sigma}, \nabla v)_K + (\nabla \cdot \boldsymbol{\sigma}, v)_K \right] + (\xi, q)_Q \le \|v, q\|_{\mathbf{V} \times \mathcal{Q}} (\|\boldsymbol{\sigma}\|_{H(\text{div};\Omega)}^2 + \|\xi\|_Q^2)^{1/2}, $$

for all $(v, q) \in \mathbf{V} \times \mathcal{Q}$ and all $\boldsymbol{\sigma} \in \mathbf{H}(\text{div}; \Omega)$ with the property $\boldsymbol{\sigma}_K n^K|_{\partial K} = \boldsymbol{\mu}_{\partial K}$ for each $K \in \mathcal{P}_H$. It follows that $\sup_{(v,q) \in \mathbf{V} \times \mathcal{Q}} \frac{b(\boldsymbol{\mu}, \xi; v, q)}{\|v, q\|_{\mathbf{V} \times \mathcal{Q}}} \le \|\boldsymbol{\mu}, \xi\|_{\mathbf{\Lambda} \times \mathcal{Q}}$ for all $(\boldsymbol{\mu}, \xi) \in \mathbf{\Lambda} \times \mathbb{R}$. Next, for each $K \in \mathcal{P}_H$, note that for each component $\mu_{\partial K,i}$ of $\boldsymbol{\mu}_{\partial K}$, $1 \le i \le d$, there exists a unique $v_{K,i}^* \in H^1(K)$ such that

$$ (\nabla v_{K,i}^*, \nabla z_K)_K + d_{\Omega}^{-2} (v_{K,i}^*, z_K)_K = \langle \mu_{\partial K,i}, z_K \rangle_{\partial K} \quad \text{for all } z_K \in H^1(K). $$
It follows that $\nabla \cdot \nabla v_{K,i}^* = d_{\Omega}^{-2} v_{K,i}^* \in L^2(K)$ and $\nabla v_{K,i}^* \cdot n^K|_{\partial K} = \mu_{\partial K,i}$. Defining $v^* \in V$ by $v^*|_K = v_K^*$, where $v_K^* := (v_{K,1}^*, \dots, v_{K,d}^*)$, it holds that $\nabla_H v^* \in H(\text{div}; \Omega)$ and,

$$ \begin{aligned} \|\boldsymbol{\mu}, \xi\|_{\Lambda \times Q} &\le (\|\nabla_H \boldsymbol{v}^*\|_{0,\Omega}^2 + d_{\Omega}^2 \|\nabla \cdot \nabla_H \boldsymbol{v}^*\|_{0,\Omega}^2 + \|\xi\|_Q^2)^{1/2} \\ &= \frac{\langle \boldsymbol{\mu}, \boldsymbol{v}^* \rangle_{\partial \mathcal{P}_H} + (\xi, \xi)_Q}{\|\boldsymbol{v}^*, \xi\|_{V \times Q}} \le \sup_{(v,q) \in V \times Q} \frac{b(\boldsymbol{\mu}, \xi; v, q)}{\|v, q\|_{V \times Q}}, \end{aligned} $$

thereby proving the second condition in (2.14). Finally, the stability of the solutions follows by classical results (see [41, Theorem 2.34]). □

**Remark 2.2** (Equivalence with primal-hybrid formulation). The solution of the classical weak formulation (1.3) and its hybrid form (2.1) coincides and $\lambda_{\partial K} = (\nu \nabla u_K - p_K) \mathbf{n}^K$ on $\partial K$ and $\rho = 0$ (c.f. [7, Theorem 1]). □

### 3. MHM'S ABSTRACT SETTING

This section presents the MHM formulation in general terms, which will be used to analyze specific cases in the following section. It also establishes necessary and sufficient conditions for the MHM method to be well-posed and with best approximation properties. To this end, consider finite subspaces of spaces $\Lambda$, $V$ and $Q$

$$ (3.1) \quad \Lambda_H := \prod_{K \in \mathcal{P}_H} \Lambda_H(\partial K), \quad V_h := \bigoplus_{K \in \mathcal{P}_H} V_h(K) \quad \text{and} \quad Q_h := \bigoplus_{K \in \mathcal{P}_H} Q_h(K), $$

where $\Lambda_H(\partial K)$, $V_h(K)$ and $Q_h(K)$ are finite subspaces of $H^{-1/2}(\partial K)^d$, $H^1(K)^d$ and $L^2(K)$. Also, let $a_{h,K}(\cdot; \cdot)$ and $f_{h,K}(\cdot)$ be bounded (bi)linears form over $V_h(K) \times Q_h(K)$. For all instances below, we assume the nullspace of $a_{h,K}(\cdot; \cdot)$ equals $\mathcal{N}_a(K)$ and $f_{h,K}(\cdot)$ coincides with $f_K(\cdot)$ over $\mathcal{N}_a(K)$. However, we distinguish its orthogonal complement from (2.12) as follows

$$ \begin{aligned} \mathcal{N}_{a,h}(K)^{\perp} &:= \{ (\mathbf{v}, q) \in \mathbf{V}_h(K) \times Q_h(K) : (\mathbf{v}, q; \mathbf{w}, r)_{\mathbf{V} \times Q} = 0 \quad \forall (\mathbf{w}, r) \in \mathcal{N}_a(K) \}, \\ \mathcal{N}_{a,h}^{\perp} &:= \{ (\mathbf{v}, q) \in \mathbf{V}_h \times Q_h : (\mathbf{v}, q)|_K \in \mathcal{N}_{a,h}(K)^{\perp} \quad \forall K \in \mathcal{P}_H \}. \end{aligned} $$

Next, we define two global mappings $T_h$ and $\hat{T}_h$ from their local counterpart, namely, $T_{h,K} : \Lambda_H(\partial K) \times \mathbb{R} \to \mathcal{N}_{a,h}^{\perp}(K)$ and $\hat{T}_{h,K} : L^2(K)^d \to \mathcal{N}_{a,h}^{\perp}(K)$ defined by

$$ (3.2) \quad a_{h,K}(T_{h,K}(\boldsymbol{\mu}, \xi); \boldsymbol{v}, q) = -b_K(\boldsymbol{\mu}, \xi; \boldsymbol{v}, q) \quad \text{and} \quad a_{h,K}(\hat{T}_{h,K}(\boldsymbol{q}); \boldsymbol{v}, q) = l_{h,K}^{\boldsymbol{q}}(\boldsymbol{v}, q), $$

for all $(\boldsymbol{v}, q) \in \mathcal{N}_{a,h}^{\perp}(K)$. Here, $l_{h,K}^{\boldsymbol{q}}(\cdot)$ is a given bounded linear form over $V_h(K) \times Q_h(K)$ associated with a $\boldsymbol{q} \in L^2(\Omega)^d$ such that $\|l_h^{\boldsymbol{q}}\| \le C\|\boldsymbol{q}\|_{0,\Omega}$, and $b_K(\cdot; \cdot)$ is given in (2.9),
where $l_h^q(\cdot)|_K := l_{h,K}^q(\cdot)$. Global versions of $a_{h,K}(\cdot; \cdot)$, $f_{h,K}(\cdot)$ and $l_{h,K}^q(\cdot)$ are defined, for all $(v, q)$, $(w, r) \in V_h \times Q_h$

$$a_h(v, q; w, r) := \sum_{K \in \mathcal{P}_H} a_{h,K}(v, q; w, r) \quad \text{and} \quad f_h(v, q) := \sum_{K \in \mathcal{P}_H} f_{h,K}(v, q).$$

Given those definitions, the **MHM abstract formulation** is: Find $(u_0^h, 0) \in \mathcal{N}_a$ and $(\lambda_H, \rho_h) \in \Lambda_H \times \mathbb{R}$ such that

$$ (3.3) \qquad \begin{aligned} & b(\boldsymbol{\mu}, \xi; \boldsymbol{\lambda}_H T_h(\boldsymbol{\lambda}_H, \rho_h)) + b(\boldsymbol{\mu}, \xi; \boldsymbol{u}_0^h, 0) = (\boldsymbol{\mu}, \boldsymbol{g})_{\partial\Omega} - b(\boldsymbol{\mu}, \xi; \hat{T}_h(\boldsymbol{f})), \\ & b(\boldsymbol{\lambda}_H, \rho_h; \boldsymbol{v}_0, 0) = f(\boldsymbol{v}_0, 0), \end{aligned} $$

for all $(\boldsymbol{\mu}, \xi) \in \Lambda_H \times \mathbb{R}$ and $(\boldsymbol{v}_0, 0) \in \mathcal{N}_a$. In the case of the Brinkman equation $(\boldsymbol{\theta} \neq \mathbf{0})$, we recall the nullspace $\mathcal{N}_a$ contains only the zero element, and then the MHM formulation (3.3) reduces to find $(\boldsymbol{\lambda}_H, \rho_h) \in \Lambda_H \times \mathbb{R}$ such that

$$b(\boldsymbol{\mu}, \xi; T_h(\boldsymbol{\lambda}_H, \rho_h)) = (\boldsymbol{\mu}, \boldsymbol{g})_{\partial\Omega} - b(\boldsymbol{\mu}, \xi; \hat{T}_h(\boldsymbol{f})) \quad \text{for all } (\boldsymbol{\mu}, \xi) \in \Lambda_H \times \mathbb{R}.$$

**3.1. Well-posedness of the abstract MHM (3.3).** Let $\mathcal{N}_{b,h}$ be the finite space

$$ (3.4) \qquad \mathcal{N}_{b,h} := \{(\boldsymbol{v}, q) \in V_h \times Q_h : b(\boldsymbol{\mu}, \xi; \boldsymbol{v}, q) = 0 \quad \text{for all } (\boldsymbol{\mu}, \xi) \in \Lambda_H \times \mathbb{R}\}, $$

and assume the following conditions hold:

**Assumption.** Given $(v_K, q_K) \in \mathcal{N}_{a,h}^{\perp}(K)$, $(z, m) \in \mathcal{N}_{b,h}$ and $(\boldsymbol{\mu}, \xi) \in \Lambda_H \times \mathbb{R}$, there exist positive constants $\alpha_{a,h}$, $\alpha_{b,h}$ and $\beta_h$ independent of mesh parameters, such that

$$ (H1) \qquad \alpha_{a,h} \|v_K, q_K\|_{H^1(K) \times L^2(K)} \le \sup_{(w_K, r_K) \in \mathcal{N}_{a,h}^{\perp}(K)} \frac{a_{h,K}(v_K, q_K; w_K, r_K)}{\|w_K, r_K\|_{H^1(K) \times L^2(K)}}, $$

$$ (H2) \qquad \alpha_{b,h} \|z, m\|_{V \times Q} \le \sup_{(w,r) \in \mathcal{N}_{b,h}} \frac{a_h(z, m; w, r)}{\|w, r\|_{V \times Q}}, $$

$$ (H3) \qquad \beta_h \|\boldsymbol{\mu}, \xi\|_{\Lambda \times Q} \le \sup_{(w,r) \in V_h \times Q_h} \frac{b(\boldsymbol{\mu}, \xi; w, r)}{\|w, r\|_{V \times Q}}. $$

Condition (H1) ensures computations may be localized to each $K \in \mathcal{P}_H$, thereby ensuring that operators $T_h$ and $\hat{T}_h$ are well-defined. Also, from (H1) the operators $T_h$ and $\hat{T}_h$ are bounded and

$$ (3.5) \qquad \alpha_{a,h} \|T_h(\boldsymbol{\mu}, \xi)\|_{V \times Q} \le \sup_{(w,q) \in \mathcal{N}_{a,h}^{\perp}} \frac{a_h(T_h(\boldsymbol{\mu}, \xi); \boldsymbol{w}, q)}{\|\boldsymbol{w}, q\|_{V \times Q}} = \sup_{(w,q) \in \mathcal{N}_{a,h}^{\perp}} \frac{b(\boldsymbol{\mu}, \xi; \boldsymbol{w}, q)}{\|\boldsymbol{w}, q\|_{V \times Q}} \le \|\boldsymbol{\mu}, \xi\|_{\Lambda \times Q}, $$

$$ \alpha_{a,h} \|\hat{T}_h(\boldsymbol{q})\|_{V \times Q} \le \sup_{(w,q) \in \mathcal{N}_{a,h}^{\perp}} \frac{a_h(\hat{T}_h(\boldsymbol{q}); \boldsymbol{w}, q)}{\|\boldsymbol{w}, q\|_{V \times Q}} = \sup_{(w,q) \in \mathcal{N}_{a,h}^{\perp}} \frac{l_h^q(\boldsymbol{w}, q)}{\|\boldsymbol{w}, q\|_{V \times Q}} \le \|l_h^q\| \le C \|\boldsymbol{q}\|_{0,\Omega}. $$
In (H2)-(H3) we recognize the necessary and sufficient condition for hybrid problem (2.1) to be well-posed over the spaces $V_h \times Q_h$ and $\Lambda_H \times \mathbb{R}$ with $a(\cdot; \cdot)$ replaced by $a_h(\cdot; \cdot)$ and $f(\cdot)$ by $f_h(\cdot)$. We recall that operators $a_h(\cdot; \cdot)$ and $f_h(\cdot)$ are assumed to be continuous on $V_h \times Q_h$ and $b(\cdot; \cdot)$ is continuous on $(\Lambda_H \times \mathbb{R}) \times (V_h \times Q_h)$ since, from the proof of Lemma 2.1, it holds

$$ \sup_{(v,q) \in V_h \times Q_h} \frac{b(\boldsymbol{\mu}, \xi; \boldsymbol{v}, q)}{\|\boldsymbol{v}, q\|_{V \times Q}} \le \sup_{(v,q) \in V \times Q} \frac{b(\boldsymbol{\mu}, \xi; \boldsymbol{v}, q)}{\|\boldsymbol{v}, q\|_{V \times Q}} \le \|\boldsymbol{\mu}, \xi\|_{\Lambda \times Q} \quad \text{for all } (\boldsymbol{\mu}, \xi) \in \Lambda_H \times \mathbb{R}. $$

If $(\mathbf{u}_h, p_h, \tilde{\boldsymbol{\lambda}}_H, \tilde{\rho}_h)$ denotes such a unique solution, from the same argument in Lemma 2.1,

$$ \begin{gathered} \|\mathbf{u}_h, p_h\|_{V \times Q} \le \frac{1}{\alpha_{b,h}} \sup_{(v,q) \in V_h \times Q_h} \frac{f_h(v,q)}{\|v,q\|_{V \times Q}} + \frac{1}{\beta_h} \left(1 + \frac{\|a_h\|}{\alpha_{b,h}}\right) \sup_{\boldsymbol{\mu} \in \Lambda_H} \frac{g(\boldsymbol{\mu}, 0)}{\|\boldsymbol{\mu}, 0\|_{\Lambda \times Q}}, \\ \|\tilde{\boldsymbol{\lambda}}_H, \tilde{\rho}_h\|_{\Lambda \times Q} \le \frac{1}{\beta_h} \left(1 + \frac{\|a_h\|}{\alpha_{b,h}}\right) \left( \sup_{(v,q) \in V_h \times Q_h} \frac{f_h(v,q)}{\|v,q\|_{V \times Q}} + \frac{\|a_h\|}{\beta_h} \sup_{\boldsymbol{\mu} \in \Lambda_H} \frac{g(\boldsymbol{\mu}, 0)}{\|\boldsymbol{\mu}, 0\|_{\Lambda \times Q}} \right). \end{gathered} $$

The following result is central to the analysis of the MHM methods that will follow. We establish that (3.3) is indeed well-posed under the assumptions (H1)-(H3). We also see that the unique solution $(\mathbf{u}_h, p_h, \tilde{\boldsymbol{\lambda}}_H, \tilde{\rho}_h)$ of (2.1) can be constructed using the unique solution of (3.3), establishing thus an equivalence between both formulations.

**Theorem 3.1.** Consider

$$ (3.7) \qquad \mathcal{N}_H := \{ (\boldsymbol{\mu}, \xi) \in \Lambda_H \times \mathbb{R} : b(\boldsymbol{\mu}, \xi; \mathbf{v}_0, 0) = 0 \quad \text{for all } (\mathbf{v}_0, 0) \in \mathcal{N}_a \}. $$

Under assumptions (H1)-(H3), it holds

$$ \begin{gathered} \frac{\alpha_{b,h} \beta_h^2}{(\alpha_{b,h} + \|a_h\|) \|a_h\|} \|\boldsymbol{\gamma}, \tau\|_{\Lambda \times Q} \le \sup_{(\boldsymbol{\mu}, \xi) \in \mathcal{N}_H} \frac{b(\boldsymbol{\mu}, \xi; T_h(\boldsymbol{\gamma}, \tau))}{\|\boldsymbol{\mu}, \xi\|_{\Lambda \times Q}}, \\ \frac{\alpha_{b,h} \beta_h}{\alpha_{b,h} + \|a_h\|} \|\mathbf{v}_0, 0\|_{V \times Q} \le \sup_{(\boldsymbol{\mu}, \xi) \in \Lambda_H \times \mathbb{R}} \frac{b(\boldsymbol{\mu}, \xi; \mathbf{v}_0, 0)}{\|\boldsymbol{\mu}, \xi\|_{\Lambda \times Q}}, \end{gathered} $$

for all $(\boldsymbol{\gamma}, \tau) \in \mathcal{N}_H$ and $(\mathbf{v}_0, 0) \in \mathcal{N}_a$. Thus, the abstract MHM formulation (3.3) admits a unique solution, and

$$ \begin{gathered} \|\boldsymbol{\lambda}_H, \rho_h\|_{\Lambda \times Q} \le C \left( \sup_{(\mathbf{v}_0, 0) \in \mathcal{N}_a} \frac{f(\mathbf{v}_0, 0)}{\|\mathbf{v}_0, 0\|_{V \times Q}} + \|\hat{T}_h(\boldsymbol{f})\|_{V \times Q} + \sup_{\boldsymbol{\mu} \in \Lambda_H} \frac{g(\boldsymbol{\mu}, 0)}{\|\boldsymbol{\mu}, 0\|_{\Lambda \times Q}} \right), \\ \|\mathbf{u}_0^h, 0\|_{V \times Q} \le C \left( \sup_{(\mathbf{v}_0, 0) \in \mathcal{N}_a} \frac{f(\mathbf{v}_0, 0)}{\|\mathbf{v}_0, 0\|_{V \times Q}} + \|\hat{T}_h(\boldsymbol{f})\|_{V \times Q} + \sup_{\boldsymbol{\mu} \in \Lambda_H} \frac{g(\boldsymbol{\mu}, 0)}{\|\boldsymbol{\mu}, 0\|_{\Lambda \times Q}} \right). \end{gathered} $$
Furthermore, $(\tilde{\lambda}_H, \tilde{\rho}_h) = (\lambda_H, \rho_h)$ and the solution $(\mathbf{u}_h, p_h)$ of hybrid problem (2.1) in $\mathbf{V}_h \times Q_h$ can be written in terms of the solution $(\mathbf{u}_0^h, 0)$ and $(\lambda_H, \rho_h)$ of (3.3) as follows

$$ (3.9) \qquad (\mathbf{u}_h, p_h) = (\mathbf{u}_0^h, 0) + T_h(\lambda_H, \rho_h) + \hat{T}_h(\mathbf{f}). $$

*Proof.* Let $(\gamma, \tau) \in \mathcal{N}_H$ and define $T_h(\gamma, \tau)$ by the first equation of (3.2), observing the need for assumption (H1). Then, define $(\mathbf{w}, r) \in \mathbf{V}_h \times Q_h$ and $(\boldsymbol{\mu}, \boldsymbol{\xi}) \in \mathbf{\Lambda}_H \times \mathbb{R}$ the unique solution (from (H2)-(H3)) of

$$ (3.10) \qquad \begin{aligned} & a_h^T(\mathbf{w}, r; \mathbf{v}, q) + b(\boldsymbol{\mu}, \boldsymbol{\xi}; \mathbf{v}, q) = (T_h(\gamma, \tau); \mathbf{v}, q)_{\mathbf{V} \times Q} \quad \text{for all } (\mathbf{v}, q) \in \mathbf{V}_h \times Q_h, \\ & b(\boldsymbol{\zeta}, \boldsymbol{\phi}; \mathbf{w}, r) = 0 \quad \text{for all } (\boldsymbol{\zeta}, \boldsymbol{\phi}) \in \mathbf{\Lambda}_H \times \mathbb{R}, \end{aligned} $$

where $a_h^T(\cdot, \cdot)$ stands for the adjoint operator of $a_h(\cdot, \cdot)$. Picking $(\mathbf{v}, q) := T_h(\gamma, \tau)$ in (3.10),

$$ \begin{aligned} \|T_h(\gamma, \tau)\|_{\mathbf{V} \times Q}^2 &= a^T(\mathbf{w}, r; T_h(\gamma, \tau)) + b(\boldsymbol{\mu}, \boldsymbol{\xi}; T_h(\gamma, \tau)) = a(T_h(\gamma, \tau); \mathbf{w}, r) + b(\boldsymbol{\mu}, \boldsymbol{\xi}; T_h(\gamma, \tau)) \\ &= -b(\gamma, \tau; \mathbf{w}, r) + b(\boldsymbol{\mu}, \boldsymbol{\xi}; T_h(\gamma, \tau)) = b(\boldsymbol{\mu}, \boldsymbol{\xi}; T_h(\gamma, \tau)), \end{aligned} $$

where we used the second equation in (3.10). The term $\|T_h(\gamma, \tau)\|_{\mathbf{V} \times Q}^2$ may be bounded from below. On one hand, from (3.6) it holds

$$ \|\boldsymbol{\mu}, \boldsymbol{\xi}\|_{\boldsymbol{\Lambda} \times Q} \le \frac{1}{\beta_h} \left(1 + \frac{\|a_h\|}{\alpha_{b,h}}\right) \sup_{(\mathbf{v}, q) \in \mathbf{V}_h \times Q_h} \frac{(T_h(\gamma, \tau); \mathbf{v}, q)_{\mathbf{V} \times Q}}{\|\mathbf{v}, q\|_{\mathbf{V} \times Q}} \le \frac{1}{\beta_h} \left(1 + \frac{\|a_h\|}{\alpha_{b,h}}\right) \|T_h(\gamma, \tau)\|_{\mathbf{V} \times Q}. $$

On the other hand, from (H3) and the first equation in (3.2) (which we use owing to the definition of $\mathcal{N}_H$), and the boundedness of $a_h(\cdot, \cdot)$ that $\frac{\beta_h}{\|a_h\|} \|\gamma, \tau\|_{\boldsymbol{\Lambda} \times Q} \le \|T_h(\gamma, \tau)\|_{\mathbf{V} \times Q}$ for all $(\gamma, \tau) \in \mathcal{N}_H$. As a result, the first inequality of (3.8) holds since the above equations imply

$$ \frac{\beta_h^2 \alpha_{b,h}}{\|a_h\| (\alpha_{b,h} + \|a_h\|)} \|\boldsymbol{\mu}, \boldsymbol{\xi}\|_{\boldsymbol{\Lambda} \times Q} \|\gamma, \tau\|_{\boldsymbol{\Lambda} \times Q} \le b(\boldsymbol{\mu}, \boldsymbol{\xi}; T_h(\gamma, \tau)). $$

We follow the same strategy to establish the second inequality of (3.8), this time using $(\mathbf{v}_0, 0) \in \mathcal{N}_a$ and defining $(\mathbf{w}, r) \in \mathbf{V}_h \times Q_h$ and $(\boldsymbol{\mu}, \boldsymbol{\xi}) \in \mathbf{\Lambda}_H \times \mathbb{R}$ as the unique solution of

$$ \begin{aligned} & a_h^T(\mathbf{w}, r; \mathbf{v}, q) + b(\boldsymbol{\mu}, \boldsymbol{\xi}; \mathbf{v}, q) = (\mathbf{v}_0, 0; \mathbf{v}, q)_{\mathbf{V} \times Q} \quad \text{for all } (\mathbf{v}, q) \in \mathbf{V}_h \times Q_h, \\ & b(\boldsymbol{\zeta}, \boldsymbol{\tau}; \mathbf{w}, r) = 0 \quad \text{for all } (\boldsymbol{\zeta}, \boldsymbol{\tau}) \in \mathbf{\Lambda}_H \times \mathbb{R}. \end{aligned} $$

Taking $(\mathbf{v}, q) := (\mathbf{v}_0, 0)$ we have

$$ \|\mathbf{v}_0, 0\|_{\mathbf{V} \times Q}^2 = a_h^T(\mathbf{w}, r; \mathbf{v}_0, 0) + b(\boldsymbol{\mu}, \boldsymbol{\xi}; \mathbf{v}_0, 0) = a_h(\mathbf{v}_0, 0; \mathbf{w}, r) + b(\boldsymbol{\mu}, \boldsymbol{\xi}; \mathbf{v}_0, 0) = b(\boldsymbol{\mu}, \boldsymbol{\xi}; \mathbf{v}_0, 0). $$

By (3.6), we get $\|\boldsymbol{\mu}, \boldsymbol{\xi}\|_{\boldsymbol{\Lambda} \times Q} \le \frac{1}{\beta_h} \left(1 + \frac{\|a_h\|}{\alpha_{b,h}}\right) \|\mathbf{v}_0, 0\|_{\mathbf{V} \times Q}$, which establishes the second inequality in (3.8) since it implies

$$ \frac{\beta_h \alpha_{b,h}}{\alpha_{b,h} + \|a_h\|} \|\mathbf{v}_0, 0\|_{\mathbf{V} \times Q} \|\boldsymbol{\mu}, \boldsymbol{\xi}\|_{\boldsymbol{\Lambda} \times Q} \le b(\boldsymbol{\mu}, \boldsymbol{\xi}; \mathbf{v}_0, 0). $$
Hence, the abstract MHM method (3.3) has a unique $(\mathbf{u}_0^h, \lambda_H, \rho_h)$ from standard saddle-point theory. Finally, it is straightforward to verify that $(\mathbf{u}_0^h, 0) + T_h(\lambda_H, \rho_h) + \hat{T}_h(\mathbf{f})$ and $(\lambda_H, \rho_h)$ satisfy the hybrid formulation (2.1) when restricted to subspaces $\mathbf{V}_h \times Q_h \times \Lambda_H \times \mathbb{R}$. Then, (3.9) follows from the uniqueness of solution of (2.1) on $\mathbf{V}_h \times Q_h \times \Lambda_H \times \mathbb{R}$ using (H2)-(H3). □

**Remark 3.1** (The continuous MHM formulation). Consider the MHM formulation (3.3) with spaces $\mathbf{V} \times Q \times \Lambda$ in place of $\mathbf{V}_h \times Q_h \times \Lambda_H$, i.e., we seek $(\mathbf{u}_0, 0) \in \mathcal{N}_a$ and $(\lambda, \rho) \in \Lambda \times \mathbb{R}$

$$ (3.11) \qquad \begin{aligned} & b(\boldsymbol{\mu}, \xi; T(\boldsymbol{\lambda}, \rho)) + b(\boldsymbol{\mu}, \xi; \mathbf{u}_0, 0) = (\boldsymbol{\mu}, \boldsymbol{g})_{\partial\Omega} - b(\boldsymbol{\mu}, \xi; \hat{T}(\boldsymbol{f})) \quad \forall (\boldsymbol{\mu}, \xi) \in \Lambda \times \mathbb{R}, \\ & b(\boldsymbol{\lambda}, \rho; \mathbf{v}_0, 0) = f(\mathbf{v}_0, 0) \quad \forall (\mathbf{v}_0, 0) \in \mathcal{N}_a, \end{aligned} $$

where (3.2) is solved using $a_h(\cdot; \cdot) := a(\cdot; \cdot)$ and $l_h^f(\mathbf{v}, q) := (\mathbf{f}, \mathbf{v})_{\mathcal{P}_H}$, and then, we adopt above the notation $T_h = T: \Lambda \times \mathbb{R} \to \mathcal{N}_a^\perp$ and $\hat{T}_h = \hat{T}: L^2(\Omega)^d \to \mathcal{N}_a^\perp$. Note that Lemma 2.1 establishes directly (H2)-(H3), and (H1) holds. Thereby, Theorem 3.1 asserts (3.11) is well-posed, and we can characterize the (unique) exact solution to the hybrid formulation (2.1) as follows

$$ (3.12) \qquad (\mathbf{u}, p) = (\mathbf{u}_0, 0) + T(\boldsymbol{\lambda}, \rho) + \hat{T}(\boldsymbol{f}) \quad \text{and} \quad \boldsymbol{\lambda}_{\partial K} = (\nu \nabla \mathbf{u}_K - p_K) \mathbf{I} \mathbf{n}^K \text{ on } \partial K. $$

Note that $\rho$ vanishes (see [7, Theorem 1] for details), and following (3.5), $T$ and $\hat{T}$ are bounded

$$ (3.13) \qquad \alpha_a \|T(\boldsymbol{\mu}, \xi)\|_{\mathbf{V} \times Q} \le \|\boldsymbol{\mu}, \xi\|_{\Lambda \times Q} \quad \text{and} \quad \alpha_a \|\hat{T}(\boldsymbol{q})\|_{\mathbf{V} \times Q} \le C \|\boldsymbol{q}\|_{0, \Omega}. $$

**3.2. Best approximation results.** Having established conditions for the MHM formulation to be well-posed, we can explore best approximation properties. Let

$$ (3.14) \qquad \Lambda_H^* := \{\boldsymbol{\mu} \in \Lambda_H : b(\boldsymbol{\mu}, 0; \mathbf{v}_0, 0) = f(\mathbf{v}_0, 0) \quad \forall (\mathbf{v}_0, 0) \in \mathcal{N}_a\}, $$

and note that $\Lambda_H^* = \Lambda_H$ in the Brinkman case.

**Theorem 3.2** (Best approximation). Assume (H1)-(H3) hold. Let $(\lambda_H, \rho_h) \in \Lambda_H \times \mathbb{R}$ and $(\mathbf{u}_0^h, 0) \in \mathcal{N}_a$ be the solution of (3.3), $(\lambda, \rho) \in \Lambda \times \mathbb{R}$ and $(\mathbf{u}_0, 0) \in \mathcal{N}_a$ the exact solution of (3.11), and $(\mathbf{u}, p)$ and $(\mathbf{u}_h, p_h)$ given in (3.12) and (3.9), respectively. Then, there exist positive constants $C$, depending only on constants in (H1)-(H3), such that

$$ \begin{aligned} \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_H\|_{\Lambda} + \|\mathbf{u}_0 - \mathbf{u}_0^h\|_{\mathbf{V}} &\le C \left( \inf_{\boldsymbol{\mu}_H \in \Lambda_H^*} \|\boldsymbol{\lambda} - \boldsymbol{\mu}_H\|_{\Lambda} + \|(T - T_h)(\boldsymbol{\lambda}, 0) + (\hat{T} - \hat{T}_h)(\boldsymbol{f})\|_{\mathbf{V} \times Q} \right), \\ \|\mathbf{u} - \mathbf{u}_h, p - p_h\|_{\mathbf{V} \times Q} &\le C \left( \inf_{\boldsymbol{\mu}_H \in \Lambda_H^*} \|\boldsymbol{\lambda} - \boldsymbol{\mu}_H\|_{\Lambda} + \|(T - T_h)(\boldsymbol{\lambda}, 0) + (\hat{T} - \hat{T}_h)(\boldsymbol{f})\|_{\mathbf{V} \times Q} \right). \end{aligned} $$
Proof. First, recall that $\rho = \rho_h = 0$. Given $\boldsymbol{\mu}_H \in \Lambda_H^*$, and using $\boldsymbol{\lambda}_H$ satisfies the second equation in (3.3), it holds $b(\boldsymbol{\mu}_H - \boldsymbol{\lambda}_H, 0; \boldsymbol{v}_0, 0) = 0$ for all $(\boldsymbol{v}_0, 0) \in \mathcal{N}_a$, and then $(\boldsymbol{\mu}_H - \boldsymbol{\lambda}_H, 0) \in \mathcal{N}_H$ where $\mathcal{N}_H$ is defined in (3.7). Therefore, from Theorem 3.1 there is $\alpha := \frac{\alpha_{b,h} \beta_h^2}{(\alpha_{b,h} + \|a_h\|) \|a_h\|}$ such that

$$
\begin{align*}
\alpha \|\boldsymbol{\mu}_H - \boldsymbol{\lambda}_H, 0\|_{\Lambda \times Q} &\le \sup_{(\boldsymbol{\mu}, \xi) \in \mathcal{N}_H} \frac{b(\boldsymbol{\mu}, \xi; T_h(\boldsymbol{\mu}_H - \boldsymbol{\lambda}_H, 0))}{\|\boldsymbol{\mu}, \xi\|_{\Lambda \times Q}} \\
&\le \sup_{(\boldsymbol{\mu}, \xi) \in \mathcal{N}_H} \frac{-b(\boldsymbol{\mu}, \xi; T(\boldsymbol{\lambda} - \boldsymbol{\mu}_H, 0)) - b(\boldsymbol{\mu}, \xi; (T - T_h)(\boldsymbol{\mu}_H, 0) + (\hat{T} - \hat{T}_h)(\boldsymbol{f}))}{\|\boldsymbol{\mu}, \xi\|_{\Lambda \times Q}} \\
&\le \frac{1}{\alpha_a} \|\boldsymbol{\lambda} - \boldsymbol{\mu}_H\|_{\Lambda} + \|(T - T_h)(\boldsymbol{\mu}_H, 0) + (\hat{T} - \hat{T}_h)(\boldsymbol{f})\|_{\mathbf{V} \times Q} \\
&\le \frac{\alpha_a + 2\alpha_{a,h}}{\alpha_a \alpha_{a,h}} \|\boldsymbol{\lambda} - \boldsymbol{\mu}_H\|_{\Lambda} + \|(T - T_h)(\boldsymbol{\lambda}, 0) + (\hat{T} - \hat{T}_h)(\boldsymbol{f})\|_{\mathbf{V} \times Q},
\end{align*}
$$

where we used (3.13), (3.5) and $\|(T - T_h)(\boldsymbol{\lambda} - \boldsymbol{\mu}_H, 0)\|_{\mathbf{V} \times Q} \le \frac{\alpha_a + \alpha_{a,h}}{\alpha_a \alpha_{a,h}} \|(\boldsymbol{\lambda} - \boldsymbol{\mu}_H, 0)\|_{\Lambda \times Q}$.
Therefore, the triangle inequality yields

$$
\|\boldsymbol{\lambda} - \boldsymbol{\lambda}_H\|_{\Lambda} \le \|\boldsymbol{\lambda} - \boldsymbol{\mu}_H\|_{\Lambda} + \|\boldsymbol{\mu}_H - \boldsymbol{\lambda}_H\|_{\Lambda} \le C \left( \inf_{(\boldsymbol{\mu}_H) \in \Lambda_H^*} \|\boldsymbol{\lambda} - \boldsymbol{\mu}_H\|_{\Lambda} + \|(T - T_h)(\boldsymbol{\lambda}, 0) + (\hat{T} - \hat{T}_h)(\boldsymbol{f})\|_{\mathbf{V} \times Q} \right).
$$

Next, observe that from Theorem 3.1 there exists $\beta := \frac{\alpha_{b,h}\beta_h}{\alpha_{b,h} + \|a_h\|}$ such that

$$
\begin{align*}
\beta \|\mathbf{u}_0 - \mathbf{u}_0^h, 0\|_{\mathbf{V} \times Q} &\le \sup_{(\boldsymbol{\mu}, \xi) \in \boldsymbol{\lambda}_H^* \times \mathbb{R}} \frac{b(\boldsymbol{\mu}, \xi; \mathbf{u}_0 - \mathbf{u}_0^h, 0)}{\|\boldsymbol{\mu}, \xi\|_{\boldsymbol{\Lambda} \times Q}} \\
&\le \sup_{(\boldsymbol{\mu}, \xi) \in \boldsymbol{\lambda}_H^* \times \mathbb{R}} \frac{-b(\boldsymbol{\mu}, \xi; (\hat{T} - \hat{T}_h)(\boldsymbol{f}) + (T - T_h)(\boldsymbol{\lambda}_H, 0) + T(\boldsymbol{\lambda} - \boldsymbol{\lambda}_H, 0))}{\|\boldsymbol{\mu}, \xi\|_{\boldsymbol{\Lambda} \times Q}} \\
&\le \|(\hat{T} - \hat{T}_h)(\boldsymbol{f}) + (T - T_h)(\boldsymbol{\lambda}_H, 0)\|_{\mathbf{V} \times Q} + \frac{1}{\alpha_a} \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_H\|_{\boldsymbol{\Lambda}} \\
&\le \|(\hat{T} - \hat{T}_h)(\boldsymbol{f}) + (T - T_h)(\boldsymbol{\lambda}, 0)\|_{\mathbf{V} \times Q} + \frac{\alpha_a + 2\alpha_{a,h}}{\alpha_a \alpha_{a,h}} \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_H\|_{\boldsymbol{\Lambda}},
\end{align*}
$$

we used again (3.13), (3.5), and the first estimate follows. The second estimate is a simple consequence of the characterization of $(\mathbf{u}, p)$ in (3.12) and $(\mathbf{u}_h, p_h)$ in (3.9), using the first estimate. $\square$

The estimates in Theorem 3.2 highlight how the choice of the second-level solver impacts approximability. It corresponds to the consistency errors $\|T - T_h\|$ and $\|\hat{T} - \hat{T}_h\|$, which vanishes in the case of the one-level MHM method (c.f. [7]). This is precisely the reason for employing the inf-sup conditions (3.8) in the proof of Theorem 3.2 instead of those of (H2)–(H3), which could be also used to estimate $\|\mathbf{u} - \mathbf{u}_h, p - p_h\|_{\mathbf{V} \times Q}$ (probably with a different error constant).
## 4. TWO-LEVEL MHM METHODS

In this section, we introduce practical strategies to approximate the solutions of the local Stokes/Brinkman problems (3.2) in order to seed the global problem of the form (3.3). The first level of a two-level MHM method is built on partitions $\mathcal{P}_H$ of $\Omega$ and $\mathcal{E}_H$ of $\mathcal{E}$. Specifically, let $0 \le \ell \in \mathbb{N}$ and define the space of (dis)continuous piecewise polynomial functions on $\mathcal{E}_H$ of degree up to $\ell$ given by

$$ \mathbb{P}_{\ell}^{\text{dis}}(\mathcal{E}_H^K)^d := \{\boldsymbol{\mu} \in L^2(\partial K)^d : \boldsymbol{\mu}_F \in \mathbb{P}_{\ell}(F)^d \text{ for all } F \in \mathcal{E}_H^K\}, $$

and

$$ \mathbb{P}_{\ell}^{\text{con}}(\mathcal{E}_H^K)^d := \{\boldsymbol{\mu} \in C^0(E)^d \text{ for all } E \subset \partial K : \boldsymbol{\mu}_F \in \mathbb{P}_{\ell}(F)^d \text{ for all } F \in \mathcal{E}_H^K\}. $$

Their global counterparts are

$$ \mathbb{P}_{\ell}^{\text{dis}}(\mathcal{E}_H)^d := \{\boldsymbol{\mu} \in L^2(\mathcal{E})^d : \boldsymbol{\mu}_F \in \mathbb{P}_{\ell}(F)^d \text{ for all } F \in \mathcal{E}_H\}, $$

and

$$ \mathbb{P}_{\ell}^{\text{con}}(\mathcal{E}_H)^d := \{\boldsymbol{\mu} \in C^0(E)^d \text{ for all } E \subset \mathcal{E} : \boldsymbol{\mu}_F \in \mathbb{P}_{\ell}(F)^d \text{ for all } F \in \mathcal{E}_H\}. $$

Then, we consider the following finite-dimensional subspace of $\Lambda$

$$ (4.1) \qquad \Lambda_H = \Lambda_H^\ell := \{\boldsymbol{\mu} \in \Lambda : \boldsymbol{\mu}_{\partial K} \in \mathbb{P}_\ell^{\text{dis}}(\mathcal{E}_H^K)^d \text{ or } \mathbb{P}_\ell^{\text{con}}(\mathcal{E}_H^K)^d \text{ for all } K \in \mathcal{P}_H\}. $$

In light of Theorem 3.2, we focus first on convergence with respect to the mesh parameters $\mathcal{H}$, $H$, and the key convergence properties will follow from the interpolation result for $\Lambda_H^\ell$. For that, the next result requires a sub-mesh $\Xi_H^K$, which denotes a regular simplicial partition for each $K \in \mathcal{P}_H$ and extends the partition on $\mathcal{E}_H$ to the interior of $K$. That is, for each $F \in \mathcal{E}_H$, there exists $\tau_F \in \Xi_H^K$ with diameter $h_{\tau_F}$ such that $\partial\tau_F \cap \partial K = F$, and for simplicity, we assume that for two different $F, F' \in \mathcal{E}_H$, we have $\tau_F \neq \tau_{F'}$.

**Lemma 4.1.** Suppose $(\mathbf{v}, q) \in H^{m+1}(\mathcal{P}_H)^d \times H^m(\mathcal{P}_H)$, $1 \le m \le \ell+1$, and $(\nu\nabla_H\mathbf{v} - q\mathbb{I}) \in \mathbf{H}(\text{div}; \Omega)$, and $\ell \ge 1$ if $\Lambda_H^\ell$ is the continuous space and $\ell \ge 0$ otherwise, and let $\boldsymbol{\mu} \in \Lambda$ be such that $\boldsymbol{\mu}_{\partial K}|_E := (\nu\nabla_K\mathbf{v}_K - q_K\mathbb{I})\boldsymbol{n}_E^K|_E$ for each $E \in \mathcal{E}$. Then, there exists $\boldsymbol{\mu}_\ell \in \Lambda_H^\ell$ and $C$ such that

$$ (4.2) \qquad \|\boldsymbol{\mu} - \boldsymbol{\mu}_\ell\|_{\Lambda} \le C H^m \left( |\nu\nabla_H \mathbf{v}|_{m, \mathcal{P}_H} + |q|_{m, \mathcal{P}_H} \right), $$

and $b(\boldsymbol{\mu} - \boldsymbol{\mu}_\ell, 0; \mathbf{w}_0, 0) = 0$ for all $\mathbf{w}_0 \in \mathbb{P}_0(\mathcal{P}_H)^d$. Moreover, assuming the quasi-uniformity of $\{\mathcal{P}_H\}_{H>0}$ and $(\mathbf{v}, q) \in H^{m+2}(\mathcal{P}_H)^d \times H^{m+1}(\mathcal{P}_H)$, such a $\boldsymbol{\mu}_\ell \in \Lambda_H^\ell$ satisfies

$$ (4.3) \qquad \|\boldsymbol{\mu} - \boldsymbol{\mu}_\ell\|_{\Lambda} \le C \frac{H^{m+1/2}}{\mathcal{H}^{1/2}} \left( \|\mathbf{v}\|_{m+2, \mathcal{P}_H} + \|q\|_{m+1, \mathcal{P}_H} \right). $$
Proof. We follow closely [47] for the discontinuous case and [60] for the continuous one using the existence of $\Xi_H^K$, for all $K \in \mathcal{P}_H$, to prove (4.2). Also, such a function $\boldsymbol{\mu}_\ell \in \Lambda_H^\ell$ satisfies $(\boldsymbol{\mu} - \boldsymbol{\mu}_\ell, \boldsymbol{v}_0)_{\partial \mathcal{P}_H} = 0$, for all $\boldsymbol{v}_0 \in \mathbb{P}_0(\mathcal{P}_H)^d$ and $\|\boldsymbol{\mu} - \boldsymbol{\mu}_\ell\|_{\partial \mathcal{P}_H} \le CH^m \|\boldsymbol{\mu}\|_{m, \partial \mathcal{P}_H}$ with $1 \le m \le \ell + 1$ and $\ell \ge 0$ in the discontinuous case and $\ell \ge 1$ otherwise. The proof of (4.3) follows the strategy in [30]. Notably, let $\Pi_K^0 : L^2(K)^d \to \mathbb{P}_0(K)^d$ be the $L^2$ projection operator on the space of constant functions on $K \in \mathcal{P}_H$, and $\Pi^0(\cdot)$ its global counterpart defined by $\Pi^0(\cdot)|_K = \Pi_K^0(\cdot)$. Then,

$$ \begin{aligned} \|\boldsymbol{\mu} - \boldsymbol{\mu}_\ell\|_{\Lambda} &= \sup_{(w,r) \in \mathbf{V} \times Q} \frac{b(\boldsymbol{\mu} - \boldsymbol{\mu}_\ell, 0; \boldsymbol{w}, r)}{\|\boldsymbol{w}, r\|_{\mathbf{V} \times Q}} = \sup_{(w,r) \in \mathbf{V} \times Q} \frac{(\boldsymbol{\mu} - \boldsymbol{\mu}_\ell, \boldsymbol{w} - \Pi^0(\boldsymbol{w}))_{\partial \mathcal{P}_H}}{\|\boldsymbol{w}, r\|_{\mathbf{V} \times Q}} \\ &\le \sup_{(w,r) \in \mathbf{V} \times Q} \frac{\|\boldsymbol{\mu} - \boldsymbol{\mu}_\ell\|_{\partial \mathcal{P}_H} \|\boldsymbol{w} - \Pi^0(\boldsymbol{w})\|_{\partial \mathcal{P}_H}}{\|\boldsymbol{w}, r\|_{\mathbf{V} \times Q}} \le CH^{1/2} \|\boldsymbol{\mu} - \boldsymbol{\mu}_\ell\|_{\partial \mathcal{P}_H} \le CH^{m+1/2} \|\boldsymbol{\mu}\|_{m, \partial \mathcal{P}_H}, \end{aligned} $$

where we used Cauchy-Schwarz inequality and the estimate $\|\boldsymbol{w}-\Pi^0(\boldsymbol{w})\|_{\partial \mathcal{P}_H} \le C H^{1/2}|\boldsymbol{w}|_{1, \mathcal{P}_H}$. Now, using $\boldsymbol{\mu}_{\partial K} = (\nu \nabla \boldsymbol{v}_K - q_K \mathbf{I})\boldsymbol{n}^K|_{\partial K}$ for all $K \in \mathcal{P}_H$, and the regularity of $\boldsymbol{v}$ and $q$ and trace inequalities (see [30] for details), we get (4.3). $\square$

The second level of the MHM method involves constructing numerical solutions to local problems using simplicial and shape-regular triangulations $\{\mathcal{T}_h^K\}_{h>0}$ of each $K \in \mathcal{P}_H$. They are construct as red refinements of the baseline mesh $\Xi_H^K$, and we denote the collection of such triangulations $\mathcal{T}_h := \bigcup_{K \in \mathcal{P}_H} \mathcal{T}_h^K$. Here, $h = \max_{\tau \in \mathcal{T}_h} h_\tau$ and $h_\tau$ is the diameter of an element $\tau \in \mathcal{T}_h$. We denote by $\mathcal{E}_h^K$ the set of faces on $\mathcal{T}_h^K$, and $\mathcal{E}_0^K$ the set of internal faces. To each face $\gamma \in \mathcal{E}_h^K$, we associate a normal vector $\boldsymbol{n}_\gamma^T$, taking care to ensure this is facing outward on $\partial \tau$. We illustrate the different partitions of the global and local domains in Figure 1.
<figure><img src="image_3.png" /><figcaption>FIGURE 1. A sketch of the global and local meshes with notations. The red dots represent the degrees of freedom associated with the $\Lambda_H$ space and the blue dots with the $V_h$ (or $Q_h$) space.</figcaption></figure>

In principle, second-level solvers may be quite general. Here, we consider two approaches:

(i) a stabilized finite element method with equal-order nodal approximation for velocity and pressure (the USFEM, [16]);

(ii) the Galerkin method using stable pairs of spaces (Taylor Hood elements, c.f. [20]).

Owing to those discrete (local) operators, we form the two-level MHM discrete solution $(\mathbf{u}_{H,h}, p_{H,h}) \in \mathbf{V}_h \times Q_h$ as follows

$$ (4.4) \qquad (\mathbf{u}_{H,h}, p_{H,h}) := (\mathbf{u}_0^{H,h}, 0) + T_h(\boldsymbol{\lambda}_{H,h}, \rho_{H,h}) + \hat{T}_h(\mathbf{f}), $$

where $(\mathbf{u}_0^{H,h}, 0) \in \mathcal{N}_a$ and $(\boldsymbol{\lambda}_{H,h}, \rho_{H,h}) \in \Lambda_H^\ell \times \mathbb{R}$ is the the solution of (3.3) with mappings $T_h$ and $\hat{T}_h$ using the methods cited in items (i) or (ii) above. Also, we recall that $(\mathbf{u}_{H,h}, p_{H,h})$ corresponds to the solution to the hybrid formulation (2.1) over $\mathbf{V}_h \times Q_h$ and $\Lambda_H^\ell \times \mathbb{R}$ according Theorem 3.1.

### 4.1. Analysis of a stabilized two-level MHM method.
The unusual stabilized finite element method (USFEM) proposed in [16] was first adopted in [7] as a second-level solver for local problems (3.2) to solve MHM formulation (3.3). For $K \in \mathcal{P}_H$, we take

$$ \mathbf{V}_h(K) := \{ \mathbf{v} \in C^0(K)^d : \mathbf{v}_\tau \in \mathbb{P}_k(\tau)^d \quad \forall \tau \in \mathcal{T}_h^K \}, $$

$$ Q_h(K) := \{ q \in C^0(K) : q_\tau \in \mathbb{P}_k(\tau) \quad \forall \tau \in \mathcal{T}_h^K \} \quad \text{with } k \ge 1. $$
The global operators $T_h \in \mathcal{L}(\Lambda_H^l \times \mathbb{R}, \mathcal{N}_{a,h}^{\perp})$ and $\hat{T}_h \in \mathcal{L}(L^2(\Omega)^d, \mathcal{N}_{a,h}^{\perp})$ are defined locally by using the USFEM in (3.2), i.e., taking $(\mathbf{u}, p), (\mathbf{v}, q) \in \mathbf{V}_h(K) \times Q_h(K)$ and the (bi)linear forms

$$ (4.5) \qquad a_{h,K}(\mathbf{u}, p; \mathbf{v}, q) := (\nu \nabla \mathbf{u}, \nabla \mathbf{v})_K + (\boldsymbol{\theta} \mathbf{u}, \mathbf{v})_K - (p, \nabla \cdot \mathbf{v})_K + (q, \nabla \cdot \mathbf{u})_K \\ \qquad \qquad - \sum_{\tau \in \mathcal{T}_h^K} \delta_{\tau}(-\nu \Delta \mathbf{u} + \boldsymbol{\theta} \mathbf{u} + \nabla p, -\nu \Delta \mathbf{v} + \boldsymbol{\theta} \mathbf{v} - \nabla q)_{\tau}, $$

$$ (4.6) \qquad l_{h,K}^{f}(\mathbf{v}, q) := (\mathbf{f}, \mathbf{v})_K - \sum_{\tau \in \mathcal{T}_h^K} \delta_{\tau}(\mathbf{f}, -\nu \Delta \mathbf{v} + \boldsymbol{\theta} \mathbf{v} - \nabla q)_{\tau}. $$

The stabilization parameter $\delta_{\tau}$ reads

$$ \delta_{\tau} := \frac{h_{\tau}^2}{\theta_{\max}^{\tau} h_{\tau}^2 \max\{1, Pe_{\tau}\} + \frac{4\nu}{m_k}} \quad \text{with} \quad Pe_{\tau} := \frac{4\nu}{\theta_{\max}^{\tau} h_{\tau}^2 m_k}, $$

where $m_k := \min\{\frac{1}{3}, C_k\}$, the non-negative constant $\theta_{\max}^{\tau} := \max_{x \in \tau} \theta_{\max}(x)$ and $C_k$ is a positive constant, independent of $h_{\tau}$, such that,

$$ C_k h_{\tau} \|\Delta \mathbf{v}\|_{0,\tau} \le \|\nabla \mathbf{v}\|_{0,\tau} \quad \text{for all } \mathbf{v} \in \mathbf{V}_h(K). $$

The USFEM has the property (see [16, Lemma 4.2] for details), given $(\mathbf{v}, q) \in \mathbf{V}_h(K) \times Q_h(K)$,

$$ (4.7) \qquad \alpha_{a,h} \|\mathbf{v}, q\|_{H^1(K) \times L^2(K)} \le \sup_{(\mathbf{w}, r) \in \mathcal{N}_{a,h}^{\perp}(K)} \frac{a_{h,K}(\mathbf{v}, q; \mathbf{w}, r)}{\|\mathbf{w}, r\|_{H^1(K) \times L^2(K)}}, $$

where $\alpha_{a,h}$ depends on physical parameters but is independent of $h_K$ and $h_{\tau}$. Furthermore, following closely the proof of Theorem 4.1 in [16] using the Cauchy-Schwartz inequality, the definition of the stabilized parameter $\delta_{\tau}$ and inverse inequalities, we obtain that $a_{h,K}(\cdot, \cdot)$ is bounded in $\mathbf{V}_h(K) \times Q_h(K)$ with a constant independent of the mesh parameters. Following closely (3.5) and using (4.7), mappings $T_h$ and $\hat{T}_h$ are bounded

$$ (4.8) \qquad \|T_h(\boldsymbol{\mu}, \boldsymbol{\xi})\|_{\mathbf{V} \times Q} \le \alpha_{a,h}^{-1} \|\boldsymbol{\mu}, \boldsymbol{\xi}\|_{\Lambda \times \mathbb{R}} \quad \text{and} \quad \|\hat{T}_h(\boldsymbol{q})\|_{\mathbf{V} \times Q} \le C \alpha_{a,h}^{-1} \|l_h^q\|. $$

Furthermore, there exists $C$ dependent only on physical coefficients, such that $\|l_h^q\| \le C \|\boldsymbol{q}\|_{0,\Omega} \Rightarrow \|\hat{T}_h(\boldsymbol{q})\|_{\mathbf{V} \times Q} \le C \alpha_{a,h}^{-1} \|\boldsymbol{q}\|_{0,\Omega}$ as a result of definition of the dual norm (2.11), the linear form $l_{h,K}^q(\cdot)$ in (4.6), and inverse inequalities. The remainder of this subsection is dedicated to the analysis of the method (3.3) with $T_h$ and $\hat{T}_h$ defined through (4.5)-(4.6). Before these may be tackled, some preliminary results must be established.
4.1.1. Preliminary results. The analysis for the stabilized two-level MHM method requires a geometrical assumption between meshes. Specifically, letting $\kappa_F$ be an element of minimal partition $\Xi_H^K$ of $K \in \mathcal{P}_H$ such that $\partial\kappa_F \cap \partial K = F$, and choose $\tau \in \mathcal{T}_h^K$ such that $\tau \subset \kappa_F$. We assume that there exists a positive constant $C$ such that

$$ (4.9) \qquad h_{\kappa_F} \le C h_{\tau}. $$

Since $h_F \le h_{\kappa_F}$, we see that the assumption (4.9) is equivalent to the requirement that the collection $\mathcal{T}_h$ cannot contain elements of size “too small” relative to the elements in $\mathcal{E}_H$. This indicates the need for a sufficiently small partition of faces in areas of the domain with strong small-scale physics, which fits in with practice. Let $\mathcal{C}_h : H^1(\mathcal{P}_H)^d \to \mathbf{V}_h$ be the Clément interpolation operator defined locally. In other words, for every $\mathbf{v} \in \mathbf{V}$ we define $\mathcal{C}_h(\mathbf{v})|_K := \mathcal{C}_h^K(\mathbf{v})$ where $\mathcal{C}_h^K : H^1(K)^d \to \mathbf{V}_h(K)$ is the usual Clément interpolation operator. It is well-known that the operator $\mathcal{C}_h^K$ satisfies the following two properties (see [41, Lemma 1.127]):

(i) there exists $C > 0$ such that

$$ (4.10) \qquad \|\mathcal{C}_h^K(\mathbf{v}_K)\|_{1,K} \le C \|\mathbf{v}_K\|_{1,K}, $$

(ii) for $m$ and $s$ satisfying $0 \le m \le s$, with $s = 0, 1$, there is $C > 0$ such that

$$ (4.11) \qquad \|\mathbf{v}_{\tau} - \mathcal{C}_h^K(\mathbf{v}_K)\|_{m,\tau} \le C h_{\tau}^{s-m} \|\mathbf{v}_{\tau}\|_{s,\omega_{\tau}^K}, $$

for all $\mathbf{v}_{\tau} \in H^s(\omega_{\tau}^K)^d$ and all $\tau \in \mathcal{T}_h^K$, where $\omega_{\tau}^K := \{\tau \in \mathcal{T}_h^K : \tau \cap \tau' \ne \emptyset\}$. The constants $C$ depend only on $k$ and $d$. The next lemma recalls the Fortin operator proposed in [47], (which is based on an argument in [15]), adding some additional information about it.

**Lemma 4.2.** Assume integers $k \ge 1$ and $\ell \ge 0$ satisfy $k - \ell \ge d$, and let $F \in \mathcal{E}_H$. Then, there exists a mapping $\Pi_h : \mathbf{V} \to \mathbf{V}_h$ such that, for all $\mathbf{v} \in \mathbf{V}$, it holds

$$ (4.12) \qquad (\Pi_h(\mathbf{v})|_F, \mu_\ell)_F = (\mathbf{v}_F, \mu_\ell)_F \quad \text{for all } \mu_\ell \in \mathbb{P}_\ell(F)^d, \\ \int_K \nabla \cdot \Pi_h(\mathbf{v})|_K d\mathbf{x} = \int_K \nabla \cdot \mathbf{v}_K d\mathbf{x} \quad \text{for all } K \in \mathcal{P}_H, $$

and $(\Pi_h(\mathbf{v}), 0) \in \mathcal{N}_{b,h}$ given in (3.4) if $\mathbf{v} \in H_0^1(\Omega)^d$. Moreover, there exist constants $C$ such that

$$ (4.13) \qquad \|\Pi_h(\mathbf{v})\|_{\mathbf{V}} \le C \|\mathbf{v}\|_{\mathbf{V}} \quad \text{and} \quad \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} h_{\tau}^{-2} \|\mathbf{v}_{\tau} - \Pi_h(\mathbf{v})|_{\tau}\|_{0,\tau}^2 \le C \|\mathbf{v}\|_{\mathbf{V}}^2. $$

*Proof.* From [47], it follows there exists a mapping $\Pi_h : \mathbf{V} \to \mathbf{V}_h$ that satisfies the first equality in (4.12) and the left inequality in (4.13). The second equality in (4.12) follows from
the first one after integration by parts and $n^K|_F \in \mathbb{P}_0(F)^d$. As for the second inequality in (4.13), we recall the definition of $\Pi_h(\cdot)$ defined locally on each $K \in \mathcal{P}_H$ as follows (see [47, Lemma 4.2])

$$ (4.14) \qquad \Pi_h(\mathbf{v})|_K := \mathcal{C}_h^K(\mathbf{v}_K) + \sum_{F \subset \partial K} \rho_F^K(\mathbf{v}_K - \mathcal{C}_h^K(\mathbf{v}_K)), $$

where $\rho_F^K$ is a mapping from $H^1(K)^d$ to $V_h(K)$ (see [47, Lemma 4.1] for details). Mapping $\rho_F^K$ is bounded, for each $\tau \in \mathcal{T}_h^K$, as follows $\|\rho_F^K \mathbf{v}_K|_\tau\|_{0,\tau} + h_\tau |\rho_F^K \mathbf{v}_K|_\tau|_{1,\tau} \le C (\|\mathbf{v}_\tau\|_{0,\tau} + h_\tau |\mathbf{v}_\tau|_{1,\tau})$. Using this, mesh regularity, (4.14) and (4.10)-(4.11), we get

$$ \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} h_\tau^{-2} \|\mathbf{v}_\tau - \Pi_h(\mathbf{v})\|_{0,\tau}^2 \le C \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} \left( h_\tau^{-2} \|\mathbf{v}_\tau - \mathcal{C}_h^K(\mathbf{v}_K)\|_{0,\tau}^2 + \|\mathbf{v}_\tau - \mathcal{C}_h^K(\mathbf{v}_K)\|_{1,\tau}^2 \right) \\ \le C \|\mathbf{v}\|_{V}^2. $$

Finally, if one restricts $\Pi_h(\cdot)$ to $H_0^1(\Omega)^d$, then from (4.12) it holds $\sum_{K \in \mathcal{P}_H} \int_{\partial K} \Pi_h(\mathbf{v})|_{\partial K} \mu_\ell|_{\partial K} \mathrm{d}\mathbf{s} = \sum_{K \in \mathcal{P}_H} \int_{\partial K} \mathbf{v}_{\partial K} \mu_\ell|_{\partial K} \mathrm{d}\mathbf{s} = 0$, for all $\mu_\ell \in \Lambda_H^\ell$, and then $(\Pi_h(\mathbf{v}), 0)$ belongs to $\mathcal{N}_{b,h}$. $\square$

The next result generalizes [45, Lemma 3.3].

**Lemma 4.3.** Let $(0, q) \in \mathcal{N}_{b,h}$ given in (3.4). There exist positive constants $C_1$ and $C_2$, independent of mesh parameters, such that for all $q \in Q$ we have

$$ (4.15) \qquad \sup_{(\mathbf{v},0) \in \mathcal{N}_{b,h}} \frac{(\nabla \cdot \mathbf{v}, q)_{\mathcal{P}_H}}{\|\mathbf{v}, 0\|_{V \times Q}} \ge C_1 \|q\|_Q - C_2 |q|_h, $$

where $|q|_h^2 := \sum_{K \in \mathcal{P}_H} |q_K|_{h,K}^2$ with $|q_K|_{h,K} := (\sum_{\tau \in \mathcal{T}_h^K} \delta_\tau \|\nabla q_\tau\|_{0,\tau}^2)^{1/2}$ is a semi-norm on $Q$.

*Proof.* Let $(0, q) \in \mathcal{N}_{b,h}$ and $\mathbf{v} \in H_0^1(\Omega)^d$ be such that $\nabla \cdot \mathbf{v} = q$ and $\|\mathbf{v}\|_{1,\Omega} \le C\|q\|_Q$. Define $(\tilde{v}, 0) := (\Pi_h(\mathbf{v}), 0) \in \mathcal{N}_{b,h}$ where $\Pi_h(\cdot)$ is the Fortin operator available from Lemma 4.2. Given $K \in \mathcal{P}_H$, let $\kappa_F$ be an element of $\Xi_H^K$ such that $\partial\kappa_F \cap \partial K = F$. Set $q_0 := \frac{1}{|\kappa_F|} \int_{\kappa_F} q_{\kappa_F} \mathrm{d}\mathbf{x}$, and define $h_{\kappa_F}$ the diameter of $\kappa_F$. From the trace and Poincaré inequalities (4.16)

$$ \|q_F - q_0|_F\|_{0,F} \le C \left( \frac{1}{h_{\kappa_F}^{1/2}} \|q_{\kappa_F} - q_0\|_{0,\kappa_F} + h_{\kappa_F}^{1/2} \|\nabla(q_{\kappa_F} - q_0)\|_{0,\kappa_F} \right) \le C h_{\kappa_F}^{1/2} \|\nabla q_{\kappa_F}\|_{0,\kappa_F}. $$
Then, by integration by parts, using (4.12), the trace inequality (4.16), (4.9) and (4.13), and the fact that $\delta_{\tau}^{-1} \le C \frac{\theta_{\max}^{\tau} + 4\nu}{h_{\tau}^2}$ (c.f. [16, estimate (39)]), we obtain

$$
\begin{align*}
(\nabla \cdot \tilde{\mathbf{v}}, q)_{\mathcal{P}_H} &= (\nabla \cdot \mathbf{v}, q)_{\mathcal{P}_H} + (\nabla \cdot (\tilde{\mathbf{v}} - \mathbf{v}), q)_{\mathcal{P}_H} \\
&= \|q\|_Q^2 + (\nabla \cdot (\tilde{\mathbf{v}} - \mathbf{v}), q)_{\mathcal{P}_H} \\
&= \|q\|_Q^2 - \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} (\tilde{\mathbf{v}}_{\tau} - \mathbf{v}_{\tau}, \nabla q_{\tau})_{\tau} + \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} (\tilde{\mathbf{v}}_{\partial\tau} - \mathbf{v}_{\partial\tau}, q_{\partial\tau} \mathbf{n}_{\tau})_{\partial\tau} \\
&= \|q\|_Q^2 - \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} (\tilde{\mathbf{v}}_{\tau} - \mathbf{v}_{\tau}, \nabla q_{\tau})_{\tau} + \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} \sum_{\gamma \subset \partial K \cap \partial \tau} (\tilde{\mathbf{v}}_{\gamma} - \mathbf{v}_{\gamma}, (q_{\gamma} - q_{0|\gamma}) \mathbf{n}_{\gamma}^K)_{\gamma} \\
&\ge \|q\|_Q^2 - \left( \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} \delta_{\tau}^{-1} \|\mathbf{v}_{\tau} - \tilde{\mathbf{v}}_{\tau}\|_{0,\tau}^2 \right)^{1/2} \left( \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} \delta_{\tau} \|\nabla q_{\tau}\|_{0,\tau}^2 \right)^{1/2} \\
&\quad - \left( \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} \sum_{\gamma \subset \partial K \cap \partial \tau} h_{\gamma}^{-1} \|\tilde{\mathbf{v}}_{\gamma} - \mathbf{v}_{\gamma}\|_{0,\gamma}^2 \right)^{1/2} \left( \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} \sum_{\gamma \subset \partial K \cap \partial \tau} h_{\gamma} \|q_{\gamma} - q_{0|\gamma}\|_{0,\gamma}^2 \right)^{1/2} \\
&\ge \|q\|_Q^2 - C \left( \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} h_{\tau}^{-2} \|\mathbf{v}_{\tau} - \tilde{\mathbf{v}}_{\tau}\|_{0,\tau}^2 \right)^{1/2} \left( \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} \delta_{\tau} \|\nabla q_{\tau}\|_{0,\tau}^2 \right)^{1/2} \\
&\quad - C \left( \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} (h_{\tau}^{-2} \|\mathbf{v}_{\tau} - \tilde{\mathbf{v}}_{\tau}\|_{0,\tau}^2 + \|\nabla(\mathbf{v}_{\tau} - \tilde{\mathbf{v}}_{\tau})\|_{0,\tau}^2) \right)^{1/2} \left( \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} \sum_{\gamma \subset \partial K \cap \partial \tau} h_{\gamma} \|q_{\gamma} - q_{0|\gamma}\|_{0,\gamma}^2 \right)^{1/2} \\
&\ge \|\tilde{\mathbf{v}}\|_{\mathbf{v}} [C_1 \|q\|_Q - C \left( \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} \delta_{\tau} \|\nabla q_{\tau}\|_{0,\tau}^2 \right)^{1/2} - C \left( \sum_{K \in \mathcal{P}_H} \sum_{F \subset \partial K} \sum_{\tau \in \mathcal{T}_h^K} h_{\tau}^2 \|\nabla q_{\tau}\|_{0,\tau}^2 \right)^{1/2}] \\
&\ge \|\tilde{\mathbf{v}}\|_{\mathbf{v}} [C_1 \|q\|_Q - C_2 \left( \sum_{K \in \mathcal{P}_H} \sum_{\tau \in \mathcal{T}_h^K} \delta_{\tau} \|\nabla q_{\tau}\|_{0,\tau}^2 \right)^{1/2}],
\end{align*}
$$

and the result follows.

<figure><img src="image_4.png" /></figure>

4.1.2. *Well-posedness.* We prove that the two-level MHM method (3.3) with the stabilized finite element method (4.5)-(4.6) is well-posed in the following theorem.

**Theorem 4.1** (Well-posedness of (3.3) with (4.5)-(4.6)). Assume that $k - \ell \ge d$ with $\ell \ge 0$ if $\Lambda_h^\ell$ is discontinuous and $\ell \ge 1$ otherwise, and assume that (4.9) is satisfied. Then, there exists a unique solution of (3.3) with the stabilized method (4.5)-(4.6), denoted by $(\mathbf{u}_0^{H,h}, 0) \in$
$\mathcal{N}_a$ and $(\boldsymbol{\lambda}_{H,h}, \rho_{H,h}) \in \boldsymbol{\Lambda}_H^\ell \times \mathbb{R}$, and

$$ \|\boldsymbol{\lambda}_{H,h}, \rho_{H,h}\|_{\boldsymbol{\Lambda}\times Q} \le C \left( \sup_{(\boldsymbol{v}_{0,0}) \in \mathcal{N}_a} \frac{(\boldsymbol{f}; \boldsymbol{v}_{0,0})_{\boldsymbol{v}\times Q}}{\|\boldsymbol{v}_{0,0}\|_{\boldsymbol{v}\times Q}} + \|\hat{T}_h(\boldsymbol{f})\|_{\boldsymbol{v}\times Q} + \sup_{\boldsymbol{\mu} \in \boldsymbol{\Lambda}_H^\ell} \frac{g(\boldsymbol{\mu}, 0)}{\|\boldsymbol{\mu}, 0\|_{\boldsymbol{\Lambda}\times Q}} \right), $$

$$ \|\boldsymbol{u}_0^{H,h}, 0\|_{\boldsymbol{v}\times Q} \le C \left( \sup_{(\boldsymbol{v}_{0,0}) \in \mathcal{N}_a} \frac{(\boldsymbol{f}; \boldsymbol{v}_{0,0})_{\boldsymbol{v}\times Q}}{\|\boldsymbol{v}_{0,0}\|_{\boldsymbol{v}\times Q}} + \|\hat{T}_h(\boldsymbol{f})\|_{\boldsymbol{v}\times Q} + \sup_{\boldsymbol{\mu} \in \boldsymbol{\Lambda}_H^\ell} \frac{g(\boldsymbol{\mu}, 0)}{\|\boldsymbol{\mu}, 0\|_{\boldsymbol{\Lambda}\times Q}} \right). $$

*Proof.* We prove conditions (H2)–(H3) as (H1) holds from (4.7). First, condition (H3) follows from the Fortin operator in Lemma 4.2 under the condition $k \ge \ell + d$, and using the $L^2$ orthogonal projection $\Pi_0(r) := \frac{1}{|\Omega|} \int_{\Omega} r \, d\boldsymbol{x}$. In fact, let $(\boldsymbol{\mu}, \boldsymbol{\xi}) \in \boldsymbol{\Lambda}_H^\ell \times \mathbb{R}$ then

$$ \|\boldsymbol{\mu}, \boldsymbol{\xi}\|_{\boldsymbol{\Lambda}\times\mathbb{R}} = \sup_{(\boldsymbol{w},r)\in\boldsymbol{V}\times Q} \frac{b(\boldsymbol{\mu}, \boldsymbol{\xi}; \boldsymbol{w}, r)}{\|\boldsymbol{w}, r\|_{\boldsymbol{V}\times Q}} \le C \sup_{(\boldsymbol{w},r)\in\boldsymbol{V}\times Q} \frac{b(\boldsymbol{\mu}, \boldsymbol{\xi}; \Pi_h(\boldsymbol{w}), \Pi_0(r))}{\|\Pi_h(\boldsymbol{w}), \Pi_0(r)\|_{\boldsymbol{V}\times Q}} \le C \sup_{(\boldsymbol{w},r)\in\boldsymbol{V}_h \times Q_h} \frac{b(\boldsymbol{\mu}, \boldsymbol{\xi}; \boldsymbol{w}, r)}{\|\boldsymbol{w}, r\|_{\boldsymbol{V}\times Q}}, $$

where we used (2.14). As for condition (H2), we adapt the strategy of the proof of [16, Lemma 4.2]. Let $(\boldsymbol{v}, q) \in \mathcal{N}_{b,h}$ and set $\boldsymbol{w} := \boldsymbol{v} - \boldsymbol{\varepsilon}\boldsymbol{z}$, where $\boldsymbol{z}$ achieves the supremum in Lemma 4.3 with properties $\|\boldsymbol{z}\|_{\boldsymbol{v}} = \|q\|_{0,\Omega}$, and $(\boldsymbol{z}, 0) \in \mathcal{N}_{b,h}$. First, [16] establishes $a_{h,K}(\boldsymbol{v}_K, q_K; \boldsymbol{v}_K, q_K) \ge \frac{3\nu}{4}\|\nabla\boldsymbol{v}_K\|_{0,K}^2 + |q_K|_{h,K}^2$ for all $K \in \mathcal{P}_H$. Further, given $K \in \mathcal{P}_H$ define $\theta_{max}^K := \max_{\tau \in \mathcal{T}_h^K} \theta_{max}^\tau$. Since $\delta_\tau \theta_{max}^\tau \le 1$ and $\delta_\tau \nu \|\Delta\boldsymbol{v}_\tau\|_{0,\tau}^2 \le \frac{1}{4}\|\nabla\boldsymbol{v}_\tau\|_{0,\tau}^2$, and using $ab \le \frac{1}{2\gamma}a^2 + \frac{\gamma}{2}b^2$ we get

$$ a_{h,K}(\boldsymbol{v}_K, q_K; -\boldsymbol{z}_K, 0) \ge -\theta_{max}^K \|\boldsymbol{v}_K\|_{0,K} \|\boldsymbol{z}_K\|_{0,K} - \nu \|\nabla\boldsymbol{v}_K\|_{0,K} \|\nabla\boldsymbol{z}_K\|_{0,K} + (\nabla \cdot \boldsymbol{z}_K, q_K)_K \\ - \sum_{\tau \in \mathcal{T}_h^K} \delta_\tau (\|\boldsymbol{\theta}\boldsymbol{v}_\tau - \nu\Delta\boldsymbol{v}_\tau\|_{0,\tau} + \|\nabla q_K\|_{0,\tau}) \|\boldsymbol{\theta}\boldsymbol{z}_\tau - \nu\Delta\boldsymbol{z}_\tau\|_{0,\tau} $$

$$ \ge -\theta_{max}^K \|\boldsymbol{v}_K\|_{0,K} \|\boldsymbol{z}\|_{0,K} - \nu \|\nabla\boldsymbol{v}_K\|_{0,K} \|\nabla\boldsymbol{z}_K\|_{0,K} + (\nabla \cdot \boldsymbol{z}_K, q_K)_K $$

$$ - \left[ \left( \sum_{\tau \in \mathcal{T}_h^K} \delta_\tau \|\boldsymbol{\theta}\boldsymbol{v}_\tau - \nu\Delta\boldsymbol{v}_\tau\|_{0,\tau}^2 \right)^{1/2} + |q_K|_{h,K} \right] \left( \sum_{\tau \in \mathcal{T}_h^K} \delta_\tau \|\boldsymbol{\theta}\boldsymbol{z}_\tau - \nu\Delta\boldsymbol{z}_\tau\|_{0,\tau}^2 \right)^{1/2} $$

$$ \ge -3 (\theta_{max}^K \|\boldsymbol{v}_K\|_{0,K}^2 + \nu \|\nabla\boldsymbol{v}_K\|_{0,K}^2)^{1/2} (\theta_{max}^K \|\boldsymbol{z}_K\|_{0,K}^2 + \nu \|\nabla\boldsymbol{z}_K\|_{0,K}^2)^{1/2} + (\nabla \cdot \boldsymbol{z}_K, q_K)_K $$

$$ - |q_K|_{h,K} \left( 2 \sum_{\tau \in \mathcal{T}_h^K} \theta_{max}^\tau \|\boldsymbol{z}_\tau\|_{0,\tau}^2 + \frac{\nu}{4} \|\nabla_H \boldsymbol{z}_\tau\|_{0,\tau}^2 \right)^{1/2} $$

$$ \ge -\frac{3}{2\gamma_1} (\theta_{max}^K \|\boldsymbol{v}_K\|_{0,K}^2 + \nu \|\nabla\boldsymbol{v}_K\|_{0,K}^2) - \frac{3\gamma_1 \tilde{C}}{2} \left( \frac{1}{d_\Omega^2} \|\boldsymbol{z}_K\|_{1,K}^2 + \|\nabla\boldsymbol{z}_K\|_{0,K}^2 \right) + (\nabla \cdot \boldsymbol{z}_K, q_K)_K $$

$$ - \frac{1}{2\gamma_2} |q_K|_{h,K}^2 - \frac{2\gamma_2 \tilde{C}}{2} \left( \frac{1}{d_\Omega^2} \|\boldsymbol{z}_K\|_{1,K}^2 + \|\nabla\boldsymbol{z}_K\|_{0,K}^2 \right), $$
where $\tilde{C} := \max\{\theta_{max}^K d_{\Omega}^2, \nu\}$. Therefore, using the lower bound on $(\nabla \cdot z_K, q_K)_K$ from Lemma 4.3, and $\|z\|_{\mathbf{V}} = \|q\|_{Q}$, we get

$$
\begin{align*}
a_h(\mathbf{v}, q; \mathbf{w}, q) &\ge \frac{3\nu}{2} \left(\frac{1}{2} - \frac{\varepsilon}{\gamma_1}\right) \|\nabla_H \mathbf{v}\|_{\mathcal{P}_H}^2 - \frac{3 c_{max} \varepsilon}{2\gamma_1} \|\mathbf{v}\|_{\mathcal{P}_H}^2 + \varepsilon(C_1 \|q\|_Q^2 - C_2 |q|_h \|q\|_Q) \\
&\quad - \frac{\tilde{C} \varepsilon}{2}(3\gamma_1 + 2\gamma_2) \|q\|_Q^2 + \left(1 - \frac{\varepsilon}{2\gamma_2}\right) |q|_h^2 \\
&\ge \frac{3\nu}{2} \left(\frac{1}{2} - \frac{\varepsilon}{\gamma_1}\right) \|\nabla_H \mathbf{v}\|_{\mathcal{P}_H}^2 - \frac{3 c_{max} \varepsilon}{2\gamma_1} \|\mathbf{v}\|_{\mathcal{P}_H}^2 + \frac{1}{2} \left(2 - \frac{\varepsilon}{\gamma_2} - \frac{C_2 \varepsilon}{\gamma_3}\right) |q|_h^2 \\
&\quad + \frac{\varepsilon}{2} \left(2C_1 - \tilde{C}(3\gamma_1 + 2\gamma_2) - C_2\gamma_3\right) \|q\|_Q^2 \\
&= \frac{3\nu}{2} \left(\frac{1}{2} - 12\tilde{C} \frac{\varepsilon}{C_1}\right) \|\nabla_H \mathbf{v}\|_{\mathcal{P}_H}^2 - 18 c_{max} \tilde{C} \frac{\varepsilon}{C_1} \|\mathbf{v}\|_{\mathcal{P}_H}^2 + \frac{\varepsilon C_1}{2} \|q\|_Q^2 \\
&\quad + \left(1 - \frac{\varepsilon}{C_1}(4\tilde{C} + C_2^2)\right) |q|_h^2,
\end{align*}
$$

where we used $\gamma_1 = \frac{C_1}{12\tilde{C}}$, $\gamma_2 = \frac{C_1}{8\tilde{C}}$, and $\gamma_3 = \frac{C_1}{2C_2}$. Now, since $(\mathbf{v}, q) \in \mathcal{N}_{b,h}$, we have the Poincaré inequality [22]

$$ \frac{C_P}{d_{\Omega}^2} \|\mathbf{v}\|_{\mathcal{P}_H}^2 \le \|\nabla_H \mathbf{v}\|_{\mathcal{P}_H}^2, $$

where $C_P$ is a positive constant independent of mesh parameters, from which we find

$$
a_{h,K}(\mathbf{v}, q; \mathbf{w}, q) \ge \frac{3\nu}{4} \left(\frac{1}{2} - 12\tilde{C} \frac{\varepsilon}{C_1}\right) \|\nabla_H \mathbf{v}\|_{\mathcal{P}_H}^2 + \frac{3}{8} \left(\frac{C_P \nu}{d_{\Omega}^2} - 48 c_{max} \tilde{C} \frac{\varepsilon}{C_1}\right) \|\mathbf{v}\|_{\mathcal{P}_H}^2 + \frac{\varepsilon C_1}{2} \|q\|_Q^2 \\
+ \left(1 - \frac{\varepsilon}{C_1}(4\tilde{C} + C_2^2)\right) |q|_h^2 \ge \frac{3\nu}{16} \left(\|\nabla_H \mathbf{v}\|_{\mathcal{P}_H}^2 + \frac{C_P}{d_{\Omega}^2} \|\mathbf{v}\|_{\mathcal{P}_H}^2\right) + \frac{\varepsilon C_1}{2} \|q\|_Q^2 \ge C \|\mathbf{v}, q\|_{\mathbf{V} \times Q}^2,
$$

where we took $\frac{\varepsilon}{C_1} = \min\{\frac{1}{48\tilde{C}}, \frac{C_P \nu}{96c_{max}\tilde{C}d_{\Omega}^2}, \frac{1}{4\tilde{C}+C_2^2}\}$. Finally, condition (H2) follows observing that

$$ \|\mathbf{w}, q\|_{\mathbf{V} \times Q} \le \|\mathbf{v}, q\|_{\mathbf{V} \times Q} + \varepsilon \|\mathbf{z}, 0\|_{\mathbf{V} \times Q} \le \|\mathbf{v}, q\|_{\mathbf{V} \times Q} + C \varepsilon \|\mathbf{0}, q\|_{\mathbf{V} \times Q} \le C \|\mathbf{v}, q\|_{\mathbf{V} \times Q}, $$

and the result follows from Theorem 3.1. $\square$

<figure><img src="image_5.png" /></figure>

**Remark 4.4** (Relaxing $k - \ell \ge d$ constraint). The constraint $k - d \ge \ell \ge 0$ in Theorem 4.1 leads to the well-posedness of eq. (3.3) in two and three-dimensional problems on simplicial meshes. We restrict our analysis in this work to simplicial $K$ elements to avoid unnecessary technical details, but Theorem 4.1 also holds for polytopal elements under the mesh conditions given in [15, Section 1.2] and [8, Section 2.1]), for instance. In the two-dimensional case, other possibilities exist to ensure the well-posedness of eq. (3.3), such as
• $T_h = \mathcal{P}_H$, and even polynomial degree $\ell$ under the constraint $k = \ell + 1$. The proof follows from [62, Lemma 4] applied to each component of the vector function $\boldsymbol{\mu}_H^\ell \in \Lambda_H^\ell$;

• meshes with “enough” refined local meshes $T_h^K$ (see [15] for precisions) under the conditions $k = \ell + 1$ with $\ell \ge 0$ or $k = \ell$ with $\ell \ge 1$. The proof follows from [15, Lemmas 4 and 5] applied to each component of the vector function $\boldsymbol{\mu}_H^\ell \in \Lambda_H^\ell$.

Finally, we note the existence of multiscale methods that avoid the above space restriction, such as the MsHDG method [40] applied to the mixed form of the Poisson model. Such an attractive feature comes at the price of designing a priori stabilization terms on faces balanced by positive stabilization parameters.

4.1.3. **Convergence.** We recall from [16, Theorem 4.1] that the mapping $T_h$ and $\hat{T}_h$ carry approximation properties. In other words, adopted element-wisely in each $K \in \mathcal{P}_H$, they approximate $T$ and $\hat{T}$ with sharp constants which are independent of mesh parameters. Specifically, assume that $(\mathbf{u}, p) \in H^{m+1}(\mathcal{P}_H)^d \times H^m(\mathcal{P}_H)$, for $1 \le m \le k$ and $k \ge 1$, it holds

$$ (4.17) \qquad \|(T - T_h)(\boldsymbol{\lambda}, \rho) + (\hat{T} - \hat{T}_h)(\boldsymbol{f})\|_{\mathbf{V} \times Q} \le C h^m \left( |\mathbf{u}|_{m+1, \mathcal{P}_H} + \|p\|_{m, \mathcal{P}_H} \right), $$

where $T_h$ and $\hat{T}_h$ are defined from (3.2) with (4.5)-(4.6). We are ready to present the main convergence result.

**Theorem 4.2** (Convergence of (3.3) with (4.5)-(4.6)). Assume $k - \ell \ge d$, with $\ell \ge 0$ ($\ell \ge 1$) if $\Lambda_H^\ell$ is discontinuous (continuous), and $(\mathbf{u}, p) \in H^{m+1}(\mathcal{P}_H)^d \times H^m(\mathcal{P}_H)$, with $1 \le m \le \ell + 1$, and $(\nu \nabla \mathbf{u} - p\mathbb{I}) \in \mathbf{H}(\mathrm{div}; \Omega)$. Then, there exist positive constants $C_i$, $i = 1, \dots, 5$, independent of $\mathcal{H}$, $H$, $h$, such that

$$ \begin{aligned} & \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_{H,h}\|_{\boldsymbol{\Lambda}} + \|\mathbf{u}_0 - \mathbf{u}_0^{H,h}\|_{\mathbf{V}} \le C_1 H^m \left( |\nu \nabla_{\mathcal{H}} \mathbf{u}|_{m, \mathcal{P}_{\mathcal{H}}} + |p|_{m, \mathcal{P}_{\mathcal{H}}} \right) + C_2 h^m \left( |\mathbf{u}|_{m+1, \mathcal{P}_{\mathcal{H}}} + \|p\|_{m, \mathcal{P}_{\mathcal{H}}} \right), \\ & \|\mathbf{u} - \mathbf{u}_{H,h}, p - p_{H,h}\|_{\mathbf{V} \times Q} \le C_3 H^m \left( |\nu \nabla_{\mathcal{H}} \mathbf{u}|_{m, \mathcal{P}_{\mathcal{H}}} + |p|_{m, \mathcal{P}_{\mathcal{H}}} \right) + C_4 h^m \left( |\mathbf{u}|_{m+1, \mathcal{P}_{\mathcal{H}}} + \|p\|_{m, \mathcal{P}_{\mathcal{H}}} \right), \\ & \text{where } (\boldsymbol{\lambda}_{H,h}, \mathbf{u}_0^{H,h}) \text{ solves (3.3) with (4.5)-(4.6), and } (\mathbf{u}_{H,h}, p_{H,h}) \text{ is given in (4.4). In addition, if } (\mathbf{u}, p) \in H^{m+2}(\mathcal{P}_{\mathcal{H}})^d \times H^{m+1}(\mathcal{P}_{\mathcal{H}}), \text{ then} \end{aligned} $$

$$ \|\mathbf{u} - \mathbf{u}_{H,h}, p - p_{H,h}\|_{\mathbf{V} \times Q} \le C_5 \frac{H^{m+1/2}}{\mathcal{H}^{1/2}} \left( \|\mathbf{u}\|_{m+2, \mathcal{P}_H} + \|p\|_{m+1, \mathcal{P}_H} \right) + C_4 h^m \left( |\mathbf{u}|_{m+1, \mathcal{P}_H} + \|p\|_{m, \mathcal{P}_H} \right). $$
*Proof.* The result is a direct consequence of Theorem 3.2, Lemma 4.1, and (4.17). $\square$

**Remark 4.5** (Local conservation). The discrete velocity $\mathbf{u}_{H,h}$ built according to (4.4) preserves the divergence-free constraint weakly. Indeed, taking $(\mathbf{w}_h, q_h)|_K = (0, 1_K) \in \mathcal{N}_{a,h}^\perp(K)$
in (4.5)-(4.6) and using $\rho_{H,h} = 0$, it holds $\int_K \nabla \cdot \mathbf{u}_{H,h} \, d\mathbf{x} = 0$ for all $K \in \mathcal{P}_H$. In addition, the local equilibrium between the numerical traction $\lambda_{H,h}$ and the external force is satisfied, i.e.,

$$ (4.18) \qquad \int_K \theta \mathbf{u}_{H,h} \, d\mathbf{x} + \int_{\partial K} \lambda_{H,h} \, ds = \int_K \mathbf{f} \, d\mathbf{x} \quad \text{for all } K \in \mathcal{P}_H, $$

where we used the second equation in (3.3) if $\theta = 0$, and the definition of $\mathbf{u}_{H,h}$ in (4.4) and local problems (3.2) with (4.5)-(4.6) tested with $(\mathbf{w}_h, q_h)|_K = (\mathbf{1}_K, 0) \in \mathbf{V}_h(K) \times Q_h(K)$ otherwise.

**4.2. Analysis of a stable two-level MHM method.** Stable pairs of spaces may be used in local problems (3.2) to solve MHM formulation (3.3). Notably, we analyze here a two-level method based on the Taylor-Hood element (c.f. [20]). Recall that the Taylor-Hood element is such that given $k \ge 2$ and $K \in \mathcal{P}_H$,

$$ \begin{align*} \mathbf{V}_h(K) &:= \{ \mathbf{v}_h \in \mathbf{V}(K) \cap C^0(K)^d : \mathbf{v}_h|_{\tau} \in \mathbb{P}_k(\tau)^d \quad \forall \tau \in \mathcal{T}_h^K \}, \\ Q_h(K) &:= \{ q_h \in Q(K) \cap C^0(K) : q_h|_{\tau} \in \mathbb{P}_{k-1}(\tau) \quad \forall \tau \in \mathcal{T}_h^K \}. \end{align*} $$

We then have a fully defined MHM method by letting $a_{h,K}(\cdot; \cdot) := a_K(\cdot; \cdot)$ and $l_{h,K}^f(\cdot) := (\mathbf{f}, \cdot)_K$ as defined in (2.9) and (2.10), respectively. We immediately obtain that $a_{h,K}(\cdot; \cdot)$ is bounded on $\mathbf{V}_h(K) \times Q_h(K)$ from (2.11) with a constant independent of the mesh parameters. To see that the method (3.3) is well-posed and converges optimally using these definitions, we first note that assumption (H1) is met when using spaces $\mathbf{V}_h(K)$ and $Q_h(K)$ (c.f. [18, 19]), with a constant independent of $K$, when the sub-mesh $\mathcal{T}_h^K$ satisfies the following mild conditions:

**Assumption (Assumption (M)).**

2D case: $\mathcal{T}_h^K$ contains at least one internal vertex;

3D case: Every tetrahedron in $\mathcal{T}_h^K$ has at least one internal vertex.

**Lemma 4.6.** Given $K \in \mathcal{P}_H$, assume that $\mathcal{T}_h^K$ satisfies Assumption (M). Then, condition (H1) holds, and

$$ (4.19) \qquad \|(T - T_h)(\boldsymbol{\lambda}, \rho) + (\hat{T} - \hat{T}_h)(\mathbf{f})\|_{\mathbf{V} \times Q} \le C h^m (\|\mathbf{u}\|_{m+1, \mathcal{P}_H} + \|p\|_{m, \mathcal{P}_H}), $$

where $1 \le m \le k$, and $k \ge 2$. Moreover, if $\mathcal{T}_h^K$ is a one-element mesh, then condition (H1) and estimate (4.19) also hold.

**Proof.** We use the standard technique for mixed problems to prove (H1) (see [20]). First, note that $(\nu \nabla \mathbf{v}_h, \nabla \mathbf{v}_h)_K + (\boldsymbol{\theta} \mathbf{v}_h, \mathbf{v}_h)_K \ge C \|\mathbf{v}_h\|_{1,K}^2$ for all $(\mathbf{v}_h, 0) \in \mathcal{N}_{a,h}^\perp(K)$ where we used
Poincaré inequality for the Stokes case ($\theta = 0$). Next, the pair of spaces $V_h(K) \times Q_h(K)$ (or $[V_h(K) \cap L_0^2(K)^d] \times Q_h(K)$) is compatible in the sense of Babuska-Brezzi using [18, Theorem 4.1] and [19, Theorem 3.1] under the assumption (M) on the partition $\mathcal{T}_h^K$. When $\mathcal{T}_h^K$ is a mesh of one element, the corresponding inf-sup conditions follows noting that the image of the Raviart-Thomas space $RT_k(K) \subset V_h(K)$ by the divergence operator coincides with $\mathbb{P}_{k-1}(K)$. Therefore, we recover the condition (H1) using [66, Theorem 3]. Finally, classical arguments based on Galerkin orthogonality, continuity and interpolation results applied to each $K \in \mathcal{P}_H$ result in (4.19). $\square$

The two-level MHM method (3.3) with the stable second-level solver (3.2) is well-posed and convergent. This is established in the next theorem.

**Theorem 4.3.** Assume that $\mathcal{T}_h^K$ satisfies Assumption (M). Then, the method (3.3) with the local problems (3.2) defined using (2.9)-(2.10) is well-posed with $k - \ell \ge d$ and $\ell \ge 0$ if $\Lambda_H^\ell$ is discontinuous and $\ell \ge 1$ otherwise. Moreover, if $(u, p) \in H^{m+1}(\mathcal{P}_H)^d \times H^m(\mathcal{P}_H)$, with $1 \le m \le \ell + 1$, and $(\nu \nabla u - p\mathbb{I}) \in \mathbf{H}(\mathrm{div}; \Omega)$ then there exist $C$ such that

$$ (4.20) \qquad \begin{aligned} & \|\lambda - \lambda_{H,h}\|_{\Lambda} + \|u_0 - u_0^{H,h}\|_{\mathbf{V}} \le C(H^m + h^m) \left( \|u\|_{m+1,\mathcal{P}_H} + \|p\|_{m,\mathcal{P}_H} \right), \\ & \|u - u_{H,h}\|_{\mathbf{V}} + \|p - p_{H,h}\|_{Q} \le C(H^m + h^m) \left( \|u\|_{m+1,\mathcal{P}_H} + \|p\|_{m,\mathcal{P}_H} \right), \end{aligned} $$

where $(\lambda_{H,h}, u_0^{H,h})$ solves (3.3) with (2.9)-(2.10), and $(u_{H,h}, p_{H,h})$ is given in (4.4).

*Proof.* First, observe that (H3) holds by the arguments used in the proof of Theorem 4.1 and (H1) follows from Lemma 4.6. It remains to show condition (H2), which is a consequence of the local inf-sup condition of the generalized Taylor-Hood element presented in the proof of Lemma 4.6, and the Fortin mapping proposed in Lemma 4.2. Specifically, note that given $(0, q) \in \mathcal{N}_{b,h}$, there exists $v \in H_0^1(\Omega)^d$ such that

$$ (4.21) \qquad \nabla \cdot v = q \quad \text{in } \Omega \quad \text{and} \quad \|v\|_{1,\Omega} \le C \|q\|_Q. $$

Then, we adapt the technique of [46, Lemma 3.5], and define the following “corrector” mapping $\mathcal{M}_h: H_0^1(\Omega)^d \to V_h \cap H_0^1(\mathcal{P}_H)^d$ as

$$ (4.22) \qquad \int_K r_h \nabla \cdot \mathcal{M}_h(v)|_K d\mathbf{x} = \int_K r_h \nabla \cdot (v_K - \Pi_h(v)|_K) d\mathbf{x} \quad \text{for all } r_h \in Q_h(K) \cap L_0^2(K), $$

where $H_0^1(\mathcal{P}_H)^d$ denotes the space of functions $v_K \in H_0^1(K)^d$, for all $K \in \mathcal{P}_H$. Observe that $\mathcal{M}_h(v)|_K \in V_h(K) \cap H_0^1(K)^d$ exists under the mesh assumption (M) which allows the
use of [18, Theorem 4.1] and [19, Theorem 3.1], and then

$$ (4.23) \qquad \|\mathcal{M}_h(\mathbf{v})\|_{\mathbf{V}} \le C \|\mathbf{v}\|_{\mathbf{V}}. $$

Next, using $k \ge d + \ell \ge 2$ by assumption, we set $\Theta_h: H_0^1(\Omega)^d \to \mathbf{V}_h$ as $\Theta_h(\mathbf{v}) := \Pi_h(\mathbf{v}) + \mathcal{M}_h(\mathbf{v})$, and then from (4.13), (4.23) and (4.21), we get $\|\Theta_h(\mathbf{v})\|_{\mathbf{V}} \le \|\Pi_h(\mathbf{v})\|_{\mathbf{V}} + \|\mathcal{M}_h(\mathbf{v})\|_{\mathbf{V}} \le C \|\mathbf{v}\|_{\mathbf{V}} \le C\|q\|_Q$. Notice from the properties of operator $\Pi_h(\cdot)$ in (4.12) and $\mathcal{M}_h(\mathbf{v})|_{\partial K} = 0$ for all $K \in \mathcal{P}_H$, that for all $\mu_h \in \Lambda_H^\ell$ it holds $(\mu_h, \Theta_h(\mathbf{v}))_{\partial \mathcal{P}} = (\mu_h, \Pi_h(\mathbf{v}))_{\partial \mathcal{P}} = 0 \Rightarrow (\Theta_h(\mathbf{v}), 0) \in \mathcal{N}_{b,h}$. In addition, for all $q \in Q_h$, $q = \tilde{q} + q^0$ with $\tilde{q} \in L_0^2(\mathcal{P}_H)$ and $q^0 \in \mathbb{P}_0(\mathcal{P}_H)$,

$$ \begin{aligned} \int_K q_K \nabla \cdot \Theta_h(\mathbf{v})|_K d\mathbf{x} &= \int_K \tilde{q}_K \nabla \cdot \Theta_h(\mathbf{v})|_K d\mathbf{x} + \int_K q_K^0 \nabla \cdot \Theta_h(\mathbf{v})|_K d\mathbf{x} \\ &= \int_K \tilde{q}_K \nabla \cdot \mathbf{v}_K d\mathbf{x} + \int_K q_K^0 \nabla \cdot \Pi_h(\mathbf{v})|_K d\mathbf{x} + \int_K q_K^0 \nabla \cdot \mathcal{M}_h(\mathbf{v})|_K d\mathbf{x} \\ &= \int_K q_K \nabla \cdot \mathbf{v}_K d\mathbf{x}, \end{aligned} $$

where we used the definition of operator $\mathcal{M}_h(\cdot)$ in (4.22), Lemma 4.2 and $\int_K q_K^0 \nabla \cdot \mathcal{M}_h(\mathbf{v})|_K d\mathbf{x} = 0$ since $\mathcal{M}_h(\mathbf{v})|_K \in H_0^1(K)^d$. We conclude that, given $(0, q) \in \mathcal{N}_{b,h}$, there exists $(\tilde{v}, 0) := (\Theta_h(\mathbf{v}), 0) \in \mathcal{N}_{b,h}$ such that

$$ (4.24) \qquad (q, \nabla \cdot \tilde{\mathbf{v}})_\Omega = (q, \nabla \cdot \mathbf{v})_\Omega = \|q\|_Q^2 \quad \text{and} \quad \|\tilde{\mathbf{v}}\|_{\mathbf{V}} \le C\|q\|_Q. $$

Let $(\mathbf{v}, q) \in \mathcal{N}_{b,h}$ and define $\mathbf{w} := \mathbf{v} - \varepsilon \tilde{\mathbf{v}}$. We observe $(\mathbf{w}, q) \in \mathcal{N}_{b,h}$ and then there exists positive constants $C_1$ and $C_2$ independent of mesh parameters, such that

$$ \begin{aligned} a_h(\mathbf{v}, q; \mathbf{w}, q) &\ge C_1 \|\mathbf{v}\|_{\mathbf{V}}^2 - \varepsilon \left[ \nu(\nabla \mathbf{v}, \nabla \tilde{\mathbf{v}})_{\mathcal{P}_H} + (\boldsymbol{\theta} \mathbf{v}, \tilde{\mathbf{v}})_{\mathcal{P}_H} - (q, \nabla \cdot \tilde{\mathbf{v}})_{\mathcal{P}_H} \right] \\ &\ge C_1 \|\mathbf{v}\|_{\mathbf{V}}^2 - \varepsilon C_2 \|\mathbf{v}\|_{\mathbf{V}} \|\tilde{\mathbf{v}}\|_{\mathbf{V}} + \varepsilon \|q\|_Q^2 \\ &\ge C_1 \|\mathbf{v}\|_{\mathbf{V}}^2 - \varepsilon C_2 C \|\mathbf{v}\|_{\mathbf{V}} \|q\|_Q + \varepsilon \|q\|_Q^2 \\ &\ge \left( C_1 - \frac{C_2 C \varepsilon}{2 \gamma} \right) \|\mathbf{v}\|_{\mathbf{V}}^2 + \varepsilon \left( 1 - \frac{C_2 C \gamma}{2} \right) \|q\|_Q^2 \\ &\ge C_3 \|\mathbf{v}, q\|_{\mathbf{V} \times Q}^2, \end{aligned} $$

where we chose the positive constants $\gamma$ and $\varepsilon$ such that $C_3$ is a positive constant independent of mesh parameters. In addition, we have $\|\mathbf{w}, q\|_{\mathbf{V} \times Q} \le \|\mathbf{v}, q\|_{\mathbf{V} \times Q} + \varepsilon \|\tilde{\mathbf{v}}, 0\|_{\mathbf{V} \times Q} \le \|\mathbf{v}, q\|_{\mathbf{V} \times Q} + C \varepsilon \|0, q\|_{\mathbf{V} \times Q} \le C \|\mathbf{v}, q\|_{\mathbf{V} \times Q}$, and the well-posedness follows from Theorem 3.1. The error estimates (4.20) follow from Theorem 3.2, Lemma 4.1 and estimate (4.19) using $\min \{\ell + 1, k\} = \ell + 1$. $\square$

**Remark 4.7** (The one-element sub-mesh case). We note that if $T_h^K$ is a mesh of an element $(\mathcal{H} = H = h)$, then the MHM method with the stable Taylor-Hood element at the second level coincides with the non-conforming Galerkin method with the Crouzeix-Raviart element for
polynomial degrees $k = 1, 2, 3$. In fact, note that the MHM method can be recast as find
$(\mathbf{u}_h, p_h) \in \mathcal{N}_{b,h}$

$$
a_h(\mathbf{u}_h, p_h; \mathbf{v}_h, q_h) = (\mathbf{f}, \mathbf{v}_h)_\Omega \quad \text{for all } (\mathbf{v}_h, q_h) \in \mathcal{N}_{b,h},
$$

where the product space with weakly continuous polynomial velocity and pressure with zero
mean value $\mathcal{N}_{b,h}$ is given in (3.4). Then, we recognize [36] for $\ell = 0$ and $k = 1$ ($d = 2, 3$),
[44] for $\ell = 1$ and $k = 2$ ($d = 2$) and [35] for $\ell = 2$ and $k = 3$ ($d = 2$). Consequently, the
two-level stable MHM method using one-element sub-meshes is well-posed for those cases and
satisfies estimate (4.20) (with $H = h$). The pair of spaces $\mathbf{V}_h \times Q_h$ with polynomial degree
$k \ge \ell + d$ and $\ell \ge 2$ in $2D$ and $\ell \ge 0$ in $3D$ are inf-sup stable as it includes the (stable)
Scott-Vogelius element (c.f. [64, 67]). Then, the two-level MHM method is also well-posed
in those cases and satisfies (4.20).

5. NUMERICAL BENCHMARKS

In this section, we assess the two-level stabilized MHM method (3.3) using USFEM as
a local solver (Section 4.1), performing convergence tests and verifying that this method is
robust when simulating fluid flows in highly heterogeneous porous media. The validations
focus on two-dimensional problems, but since the theoretical estimates of the MHM also
work for three-dimensional problems, equivalent numerical results are expected in this case
(see the validation for three-dimensional linear elasticity in [47] and for a transport model in
[8], for example). Details of the underlying algorithm and its implementation can be found
in [7, Algorithm 1].

5.1. **Convergence studies.** We first focus on the Stokes model ($\boldsymbol{\theta} = \mathbf{0}$) with $\nu = 1$. The domain $\Omega$ is $]0, 1[\times]0, 1[$, the function $\boldsymbol{f}$ is chosen such that the exact solution is given by

$$
u_1(x, y) = -u_2(x, y) = -256 x^2 (x-1)^2 y (y-1) (2y-1), \quad p(x, y) = 150 \left(x - \frac{1}{2}\right) \left(y - \frac{1}{2}\right).
$$

The discontinuous version of $\Lambda_H^\ell$ was validated in [7], and then we assume hereafter that $\Lambda_H^\ell$
is the space of piecewise continuous polynomial functions. We set $\ell = 1$ and first validate the
method using USFEM with the elements $\mathbb{P}_3(K)^2 \times \mathbb{P}_3(K)$ in one-element sub-meshes (i.e.,
$\mathcal{H} = H = h$), respectively. We observe convergence rates concerning $H$ in Figure 2 (left)
that perfectly agree with Theorem 4.2. Using the broken norm in the $\mathbf{H}(\text{div}; \mathcal{P}_{\mathcal{H}})$ space,
we also verify in Figure 2 (left) the convergence of the two-level stress variable $\boldsymbol{\sigma}_{H,h} :=$
$\nu \nabla_{\mathcal{H}} \mathbf{u}_{H,h} - p_{H,h} \mathbf{I}$, and $\|\boldsymbol{\sigma} - \boldsymbol{\sigma}_{H,h}\|_{\mathbf{H}(\text{div}; \mathcal{P}_{\mathcal{H}})}^2 := \sum_{K \in \mathcal{P}_{\mathcal{H}}} \|\boldsymbol{\sigma} - \boldsymbol{\sigma}_{H,h}\|_{\mathbf{H}(\text{div}; K)}^2$, which is not
covered in Theorem 4.2. Owing to the regularity of the exact solution, we also observe that
the error in the $L^2(\Omega)$ norm for the velocity converges of order $O(H^{\ell+2})$ as expected. Next, we
maintain $\mathcal{H}$ fixed and refine the face partitions ($H \to 0$). Since the exact solution is regular, we achieve super-convergence with an extra $O(H^{1/2})$ rate in agreement with Theorem 4.2 (see Figure 2 (right)).

<figure><img src="image_6.png" /><figcaption>FIGURE 2. Convergence curves with respect to $\mathcal{H}$ (left) and $H$ (right). Here $\nu = 1$ and $\theta = 0$.</figcaption></figure>

Table 1 shows convergence when using continuous or piecewise discontinuous space $\Lambda_H^1$. To do so, we keep the number of mesh elements constant at the first level (64 elements) and refine the mesh only at the skeleton. Note that using the piecewise continuous space $\Lambda_H^1$, we obtain better error results for a given number of degrees of freedom than the discontinuous case.

<table><caption>TABLE 1. Convergence history with respect to the number of global degrees of freedom.</caption><thead><tr><th rowspan="2">dof</th><th colspan="2">||<i>p</i> − <i>p</i><sub>H,h</sub>||<sub>0,&#x2126;</sub></th><th colspan="2">||<b>u</b> − <b>u</b><sub>H,h</sub>||<sub>0,&#x2126;</sub></th><th colspan="2">||<b>u</b> − <b>u</b><sub>H,h</sub>||<sub>V</sub></th></tr><tr><th>&Lambda;<sup>1</sup><sub>H</sub> (cont)</th><th>&Lambda;<sup>1</sup><sub>H</sub> (disc)</th><th>&Lambda;<sup>1</sup><sub>H</sub> (cont)</th><th>&Lambda;<sup>1</sup><sub>H</sub> (disc)</th><th>&Lambda;<sup>1</sup><sub>H</sub> (cont)</th><th>&Lambda;<sup>1</sup><sub>H</sub> (disc)</th></tr></thead><tbody><tr><td>545</td><td>0.1978e+00</td><td>0.1978e+00</td><td>0.8244e-02</td><td>0.8244e-02</td><td>0.2857e+00</td><td>0.2857e+00</td></tr><tr><td>1377</td><td>0.3197e-02</td><td>0.1275e-01</td><td>0.2673e-04</td><td>0.1435e-03</td><td>0.5017e-02</td><td>0.1667e-01</td></tr><tr><td>2209</td><td>0.7108e-03</td><td>0.3535e-02</td><td>0.3203e-05</td><td>0.2343e-04</td><td>0.1120e-02</td><td>0.4614e-02</td></tr><tr><td>3873</td><td>0.1421e-03</td><td>0.8087e-03</td><td>0.3322e-06</td><td>0.2940e-05</td><td>0.2246e-03</td><td>0.1053e-02</td></tr></tbody></table>

Next, we are interested in the exact solution of the Brinkman problem ($\theta = \mathbb{I}$) using the same analytical solution above, just changing the definition of the exact pressure by $p(x, y) = (x - y)^6 - 1/28$. We do not show the convergence in relation to $H$ and $\mathcal{H}$ because this is close to that obtained for the Stokes problem. Instead, we show in Figure 3 the behavior of the error in relation to the diffusion coefficient $\nu$ when it tends to zero. We note that the error remains limited for a wide range of $\nu$ values.
<figure><img src="image_7.png" /><figcaption>FIGURE 3. Error sensitivity with respect to $\nu$ on a fixed mesh with 64 elements.</figcaption></figure>

This indicates that the constant in the error estimates in Theorem 4.2 depends only slightly in terms of physical coefficients, which verifies the robustness of the two-level stabilized MHM method. The analysis of such property, which appears closely related to the known accuracy of the USFEM in dealing with vanishing coefficient problems, deserves further theoretical investigation.

5.2. A highly heterogeneous case. We adopt the following version of the Brinkman model

$$-\mu \Delta u + \theta u + \nabla p = 0 \quad \text{in } \Omega, \quad \nabla \cdot u = 0 \quad \text{in } \Omega,$$

with $\theta := \mu K^{-1}$, where $K$ is the permeability tensor and $\mu$ is the fluid viscosity. We adopt a heterogenous isotropic permeability coefficient $K$, which is obtained from layer 36 in SPE10 project [32] (second dataset). The domain $\Omega := ]0, 1200[ \times ]0, 2200[$ and $\mu = 0.3$. The permeability tensor $K$ and the boundary conditions are depicted in Figure 4.

<figure><img src="image_8.png" /><figcaption>FIGURE 4. The statement of the problem (left), the permeability field $K$ in a logarithmic scale (center), and the coarse mesh with 528 triangle elements (right).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIGURE 5. Comparison between the isolines of the velocity magnitude from the reference solution (left) and the MHM solution (center). Profile of the pressure (right) from the reference (red) and the MHM solutions (blue).</figcaption></figure>

A reference solution is calculated by solving the problem using USFEM on a mesh with 1,081,344 elements (1,625,283 degrees of freedom) with $\mathbb{P}_1(K)^2 \times \mathbb{P}_1(K)$ element. We perform the calculations using the MHM method with the continuous space $\Lambda_H^1$ with 10 elements by faces on a coarse structured mesh with 528 elements, as shown in Figure 4 (right). The number of degrees of freedom is 33,040. The coarse partition faces are not aligned with changing coefficients, and multiple scales still persist within the elements. The sub-meshes contain one hundred triangles with $\mathbb{P}_3(K)^2 \times \mathbb{P}_3(K)$ interpolation. We note a good agreement between the reference and the MHM solution through the isolines of $|u_{H,h}|$, and the pressure profiles (see Figure 5).

## 6. CONCLUSION

Hybridization was used to characterize exact velocity and pressure variables in terms of local and global problems, producing new face-based multiscale numerical methods for the Stokes/Brinkman model. The Neumann local problems respond to multiscale basis functions that incorporate physical and geometric aspects of the coarse mesh, while the global formulation responds to the degrees of freedom. The present work is also the first to prove that, under conditions of local regularity, the two-level MHM method originally given in [7] is well-posed and has optimal approximation properties. Such results are previously established in an abstract way, starting at the continuous level and assuming generic second-level solvers that satisfy some local properties. Then, by particularizing the second-level solver via a stabilized and stable finite element method and (dis)continuous polynomial interpolation for the flow variable (e.g., Lagrange multiplier), we demonstrate that the underlying MHM methods are high-order super-convergent for pressure and velocity variables and release local mass conservative numerical velocity fields and numerical flux in local equilibrium with
external forces. These properties are achieved without post-processing. In the case of a sub-mesh composed of only one element, we establish relationships between the MHM method and classical non-conforming methods in the literature. The numerical tests validated the theoretical results and fully complement and support the extensive numerical validation first proposed in [7]. The robustness of the proposed MHM methods with respect to physical coefficients is verified numerically, but its precise demonstration is left for future work.

## FUNDING

The first author was partially founded by ANID/Chile through projects AFB170001 of the PIA Program: Concurso Apoyo a Centros Científicos y Tecnológicos de Excelencia con Financiamiento Basal and FONDECYT project No. 1211649. The third author was partially supported by ANID/Chile through FONDECYT project No. 1211649, and Dirección de Investigación of the Universidad Católica de la Santísima Concepción through project FGII 06/2024. The fourth author was partially supported by CNPq/Brazil No. 309173/2020-5 and FAPERJ E-26/201.182.2021, Project EOLIS (MATH-AMSUD 21-MATH-04) and Inria/France under the Inria International Chair.

## REFERENCES

[1] A. ABDULLE AND O. BUDAC, *An adaptive finite element heterogeneous multiscale method for Stokes flow in porous media*, SIAM Multiscale Model. Simul., 13 (2015), pp. 256–290, https://doi.org/10.1137/130950136. 3

[2] B. ACHCHAB, A. AGOUZAL, K. BOUIHAT, A. MAJDOUBI, AND A. SOUISSI, *Projection stabilized nonconforming finite element methods for the Stokes problem*, Numer. Methods Partial Diff. Eqns, 33 (2016), pp. 218–240, https://doi.org/10.1002/num.22083. 6

[3] A. AGOUZAL AND J. M. THOMAS, *Une méthode d'éléments finis hybrides en décomposition de domaines*, RAIRO Modél. Math. Anal. Numér., 29 (1995), pp. 749–765, https://doi.org/10.1051/m2an/1995290607491. 3

[4] R. ARAYA, G. R. BARRENECHEA, L. P. FRANCA, AND F. VALENTIN, *Stabilization arising from PGEM: A review and further developments*, Appl. Numer. Math., 59 (2009), pp. 2065–2081, https://doi.org/10.1016/j.apnum.2008.12.004. 3

[5] R. ARAYA, G. R. BARRENECHEA, AND F. VALENTIN, *Stabilized finite element methods based on multiscale enrichment for the Stokes problem*, SIAM J. Numer. Anal., 44 (2006), pp. 322–348, https://doi.org/10.1137/050623176. 3

[6] R. ARAYA, C. HARDER, D. PAREDES, AND F. VALENTIN, *Multiscale hybrid-mixed methods*, SIAM J. Numer. Anal., 51 (2013), pp. 3505–3531, https://doi.org/10.1137/120888223. 3

[7] R. ARAYA, C. HARDER, A. POZA, AND F. VALENTIN, *Multiscale hybrid-mixed method for the Stokes and Brinkman equations – The method*, Comput. Methods Appl. Mech. Eng., 324 (2017), pp. 29–53, https://doi.org/10.1016/j.cma.2017.05.027. 4, 5, 6, 9, 10, 11, 15, 16, 19, 30, 33, 34
[8] R. ARAYA, F. JAILLET, D. PAREDES, AND F. VALENTIN, *Generalizing the multiscale hybrid-mixed method for reactive-advection-diffusive equations*, Comput. Methods Appl. Mech. Engrg., 428 (2024), p. 117089. 3, 25, 30

[9] R. ARAYA, R. REBOLLEDO, AND F. VALENTIN, *On a multiscale a posteriori error estimator for the Stokes and Brinkman equations*, IMA J. Numer. Anal., 41 (2019), pp. 344–380, https://doi.org/10.1093/imanum/drz053. 4

[10] T. ARBOGAST AND K. BOYD, *Subgrid upscaling and mixed multiscale finite elements*, SIAM J. Numer. Anal., 44 (2006), pp. 1150–1171, https://doi.org/10.1137/050631811. 3

[11] T. ARBOGAST, G. PENCHEVA, M. F. WHEELER, AND I. YOTOV, *A multiscale mortar mixed finite element method*, SIAM Multiscale Model. Simul., 6 (2007), pp. 319–346, https://doi.org/10.1137/060662587. 3, 4

[12] D. N. ARNOLD AND F. BREZZI, *Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates*, RAIRO Modél. Math. Anal. Numér., 19 (1985), pp. 7–32, https://doi.org/10.1051/m2an/1985190100071. 3

[13] I. BABUŠKA AND E. OSBORN, *Generalized finite element methods: Their performance and their relation to mixed methods*, SIAM J. Num. Anal., 20 (1983), pp. 510–536, https://doi.org/10.1137/0720034. 3

[14] G. R. BARRENECHEA, A. T. A. GOMES, AND D. PAREDES, *A multiscale hybrid method*, SIAM Journal on Scientific Computing, 46 (2024), pp. A1628–A1657, https://doi.org/10.1137/22M1542556. 6

[15] G. R. BARRENECHEA, F. JAILLET, D. PAREDES, AND F. VALENTIN, *The multiscale hybrid mixed method in general polygonal meshes*, Numer. Math., 145 (2020), pp. 197–237, https://doi.org/10.1007/s00211-020-01103-5. 3, 21, 25, 26

[16] G. R. BARRENECHEA AND F. VALENTIN, *An unusual stabilized finite element method for a generalized Stokes problem*, Numer. Math., 92 (2002), pp. 653–677, https://doi.org/10.1007/s002110100371. 3, 5, 19, 20, 23, 24, 26

[17] G. R. BARRENECHEA AND F. VALENTIN, *Consistent local projection stabilized finite element methods*, SIAM J. Numer. Anal., 48 (2010), pp. 1801–1825, https://doi.org/10.1137/090753334. 3, 6

[18] D. BOFFI, *Stability of higher order triangular Hood-Taylor methods for the stationary Stokes equations*, Math. Methods Appl. Sci., 4 (1994), pp. 223–235, https://doi.org/10.1142/S0218202594000133. 27, 28, 29

[19] D. BOFFI, *Three-dimensional finite element methods for the Stokes problem*, SIAM J. Numer. Anal., 34 (1997), pp. 664–670, https://doi.org/10.1137/S0036142994270193. 27, 28, 29

[20] D. BOFFI, F. BREZZI, AND M. FORTIN, *Mixed finite element methods and applications*, vol. 44 of Springer Series in Computational Mathematics, Springer, Heidelberg, 2013, https://doi.org/10.1007/978-3-642-36519-5. 2, 19, 27

[21] W. BOON, D. GLASER, R. HELMING, AND I. YOTOV, *Flux-mortar mixed finite element methods on nonmatching grids*, SIAM J. Numer. Anal., 60 (2022), pp. 1193–1225, https://doi.org/10.1137/20M1361407. 4

[22] S. BRENNER, *Poincaré-Friedrichs inequalities for piecewise H¹ functions*, SIAM J. Numer. Anal., 41 (2003), pp. 306–324, https://doi.org/10.1137/S0036142902401311. 25
[23] F. BREZZI, M. O. BRISTEAU, L. P. FRANCA, M. MALLET, AND G. ROGÉ, *A relationship between stabilized finite element methods and the Galerkin method with bubble functions*, Comput. Methods Appl. Mech. Engrg., 96 (1992), pp. 117-129, https://doi.org/10.1016/0045-7825(92)90102-P. 3

[24] F. BREZZI, L. P. FRANCA, T. J. R. HUGHES, AND A. RUSSO, $b = \int g$, Comput. Methods Appl. Mech. Engrg., 145 (1997), pp. 329-339, https://doi.org/10.1016/S0045-7825(96)01221-2. 3

[25] F. BREZZI, T. J. R. HUGHES, D. MARINI, AND A. MASUD, *Mixed discontinuous Galerkin methods for Darcy flow*, J. Sci. Comput., 22 (2005), pp. 119-145, https://doi.org/10.1007/s10915-004-4150-8. 3

[26] F. BREZZI AND D. MARINI, *Error estimates for the three-field formulation with bubble stabilization*, Math. Comp., 70 (2000), pp. 911-934, https://doi.org/10.1090/S0025-5718-00-01250-3. 3

[27] F. BREZZI AND L. D. MARINI, *A three-field domain decomposition method*, in Domain Decomposition in Science and Engineering, A. Quarteroni et al., ed., Providence, 1994, American Mathematical Society, pp. 27-34, https://doi.org/10.1090/conm/157/01402. 3

[28] E. BURMAN AND P. HANSBO, *A stabilized non-conforming finite element method for incompressible flow*, Comput. Methods Appl. Mech. Engrg., 195 (2006), pp. 2881-2899, https://doi.org/10.1016/j.cma.2004.11.033. 6

[29] T. CHAUMONT-FRELET, A. ERN, S. LEMAIRE, AND F. VALENTIN, *Bridging the multiscale hybrid-mixed and multiscale hybrid high-order methods*, ESAIM: Math. Model. Num. Anal., 56 (2022), pp. 261-285, https://doi.org/10.1051/m2an/2021082. 3, 6

[30] T. CHAUMONT-FRELET, D. PAREDES, AND F. VALENTIN, *Flux approximation on unfitted meshes and application to multiscale hybrid-mixed methods*, tech. report, hal-03834748, Oct. 2022, https://inria.hal.science/hal-03834748. 18

[31] T. CHAUMONT-FRELET AND F. VALENTIN, *A multiscale hybrid-mixed method for the Helmholtz equation in heterogeneous domains*, SIAM J. Numer. Anal., 58 (2020), pp. 1029-1067, https://doi.org/10.1137/19M1255616. 3

[32] M. CHRISTIE AND M. BLUNT, *Tenth SPE comparative solution project: a comparison of upscaling techniques*, PE Reservoir Evaluation and Engineering, 4s (2001), pp. 308-317, https://doi.org/10.2118/72469-PA. 32

[33] B. COCKBURN AND J. GOPALAKRISHNAN, *The derivation of hybridizable discontinuous Galerkin methods for Stokes flow*, SIAM J. Numer. Anal., 47 (2009), pp. 1092-1125, https://doi.org/10.1137/080726653. 3

[34] A. L. G. A. COUTINHO, L. P. FRANCA, AND F. VALENTIN, *Numerical multiscale methods*, Int. J. Numer. Methods Fluids, 70 (2012), pp. 403-419, https://doi.org/10.1002/fld.2727. 3

[35] M. CROUZEIX AND R. FALK, *Nonconforming finite elements for the Stokes problem*, Math. Comp., 52 (1989), pp. 437-456, https://doi.org/10.1090/S0025-5718-1989-0958870-8. 30

[36] M. CROUZEIX AND P. RAVIART, *Conforming and nonconforming finite element methods for solving the stationary Stokes equations*, R.A.I.R.O., 7 (1973), pp. 33-75, https://doi.org/10.1051/m2an/197307R300331. 6, 30

[37] W. E AND B. ENGQUIST, *The heterogeneous multiscale methods*, Comm. Math. Sci., 1 (2003), pp. 87-133, https://doi.org/10.4310/CMS.2003.v1.n1.a8. 3
[38] Y. EFENDIEV, J. GALVIS, AND T. Y. HOU, *Generalized multiscale finite element methods*, J. Comput. Phys., 251 (2013), pp. 116–135, https://doi.org/10.1016/j.jcp.2013.04.045. 3

[39] Y. EFENDIEV, R. LAZAROV, M. MOON, AND K. SHI, *A spectral multiscale hybridizable discontinuous Galerkin method for second order elliptic problems*, Comp. Methods Appl. Mech. Eng., 292 (2015), pp. 243–256, https://doi.org/10.1016/j.cma.2014.09.036. 3

[40] Y. EFENDIEV, R. LAZAROV, AND K. SHI, *A multiscale hdg method for second order elliptic equations. part i. polynomial and homogenization-based multiscale spaces*, SIAM J. Numer. Anal., 53 (2015), pp. 342–369. 26

[41] A. ERN AND J.-L. GUERMOND, *Theory and practice of finite elements*, vol. 159 of Applied Mathematical Sciences, Springer-Verlag, New York, 2004, https://doi.org/10.1007/978-1-4757-4355-5. 11, 21

[42] C. FARHAT, I. HARARI, AND L. P. FRANCA, *The discontinuous enrichment method*, Comput. Methods Appl. Mech. Engrg., 190 (2001), pp. 6455–6479, https://doi.org/10.1016/S0045-7825(01)00232-8. 3

[43] C. FARHAT, J. MANDEL, AND F.-X. ROUX, *Optimal convergence properties of the FETI domain decomposition method*, Comput. Methods Appl. Mech. Engrg., 115 (1994), pp. 365–385, https://doi.org/10.1016/0045-7825(94)90068-X. 3

[44] M. FORTIN AND M. SOULIE, *A non-conforming piecewise quadratic finite element on triangles*, Internat. J. Numer. Methods Engrg., 19 (1983), pp. 505–520, https://doi.org/10.1002/nme.1620190405. 30

[45] L. P. FRANCA AND R. STENBERG, *Error analysis of Galerkin least squares methods for the elasticity equations*, SIAM J. Numer. Anal., 28 (1991), pp. 1680–1697, https://doi.org/10.1137/0728084. 22

[46] V. GIRAULT, D. VASSILEV, AND I. YOTOV, *Mortar multiscale finite element methods for Stokes-Darcy flows*, Numer. Math., 165 (2014), pp. 93–165, https://doi.org/10.1007/s00211-013-0583-z. 3, 28

[47] A. T. A. GOMES, W. PEREIRA, AND F. VALENTIN, *The MHM method for linear elasticity on polytopal meshes*, IMA J. Numer. Anal., 43 (2022), pp. 2265–2298, https://doi.org/10.1093/imanum/drac041. 3, 18, 21, 22, 30

[48] C. HARDER, A. MADUREIRA, AND F. VALENTIN, *A hybrid-mixed method for elasticity*, ESAIM: Math. Model. Num. Anal., 50 (2016), pp. 311–336, https://doi.org/10.1051/m2an/2015046. 3

[49] C. HARDER, D. PAREDES, AND F. VALENTIN, *A family of multiscale hybrid-mixed finite element methods for the Darcy equation with rough coefficients*, J. Comput. Phys., 245 (2013), pp. 107–130, https://doi.org/10.1016/j.jcp.2013.03.019. 3

[50] C. HARDER, D. PAREDES, AND F. VALENTIN, *On a multiscale hybrid-mixed method for advective-reactive dominated problems with heterogenous coefficients*, SIAM Multiscale Model. and Simul., 13 (2015), pp. 491–518, https://doi.org/10.1137/130938499. 3

[51] C. HARDER AND F. VALENTIN, *Foundations of the MHM method*, in Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, G. R. Barrenechea, F. Brezzi, A. Cangiani, and E. H. Georgoulis, eds., Lecture Notes in Computational Science and Engineering, Springer, 2016, https://doi.org/10.1007/978-3-319-41640-3_13. 6

[52] T. Y. HOU, X. WU, AND Z. CAI, *Convergence of a multiscale finite element method for elliptic problems with rapidly oscillating coefficients*, Math. Comp., 68 (1999), pp. 913–943, https://doi.org/10.1090/S0025-5718-99-01077-7. 3
[53] T. J. R. HUGHES, G. R. FEIJOO, L. MAZZEI, AND J. QUINCY, *The variational multiscale method - a paradigm for computational mechanics*, Comput. Methods Appl. Mech. Engrg., 166 (1998), pp. 3-24, https://doi.org/10.1016/S0045-7825(98)00079-6. 3

[54] T. J. R. HUGHES, L. P. FRANCA, AND M. BALESTRA, *A new finite element formulation for computational fluid dynamics. V. Circumventing the Babuška-Brezzi condition: A stable Petrov Galerkin formulation of the Stokes problem accommodating equal-order interpolations*, Comput. Methods Appl. Mech. Engrg., 59 (1986), pp. 85-99, https://doi.org/10.1016/0045-7825(86)90025-3. 2

[55] T. J. R. HUGHES, A. MASUD, AND J. WAN, *A stabilized mixed discontinuos Galerkin method for Darcy flow*, Comput. Methods Appl. Mech. Engrg., 195 (2006), pp. 3347-3381, https://doi.org/10.1016/j.cma.2005.06.018. 3

[56] S. LANTERI, D. PAREDES, C. SCHEID, AND F. VALENTIN, *The multiscale hybrid-mixed method for the Maxwell equations in heterogeneous media*, SIAM Multiscale Model. Simul., 16 (2018), pp. 1648-1683, https://doi.org/10.1137/16M110037X. 3

[57] C. LE BRIS, F. LEGOLL, AND A. LOZINSKI, *MsFEM à la Crouzeix-Raviart for highly oscillatory elliptic problems*, Chin. Ann. Math., 34B (2013), pp. 113-138, https://doi.org/10.1007/s11401-012-0755-7. 3

[58] A. MÅLQVIST AND D. PETERSEIM, *Localization of elliptic multiscale problems*, Math. Comp., 83 (2014), pp. 2583-2603, https://doi.org/10.1090/S0025-5718-2014-02868-8. 3

[59] D. PAREDES, F. VALENTIN, AND H. M. VERSIEUX, *On the robustness of multiscale hybrid-mixed methods*, Math. Comp., 86 (2017), pp. 525-548, https://doi.org/10.1090/mcom/3108. 3

[60] D. PAREDES, F. VALENTIN, AND H. M. VERSIEUX, *Revisiting the robustness of the multiscale hybrid-mixed method: The face-based strategy*, J. Comput. Appl. Math., 436 (2024), p. 115415, https://doi.org/10.1016/j.cam.2023.115415. 3, 18

[61] P. RAVIART AND J. THOMAS, *A mixed finite element method for 2nd order elliptic problems*, Mathematical aspect of finite element methods, no. 606 in Lecture Notes in Mathematics, Springer-Verlag, New York, 1977, pp. 292-315, https://doi.org/10.1007/BFb0064470. 3

[62] P. RAVIART AND J. THOMAS, *Primal hybrid finite element methods for 2nd order elliptic equations*, Math. Comp., 31 (1977), pp. 391-413, https://doi.org/10.1090/S0025-5718-1977-0431752-8. 26

[63] G. SANGALLI, *Capturing small scales in elliptic problems using a Residual-Free Bubbles finite element method*, SIAM Multiscale Model. Simul., 1 (2003), pp. 485-503, https://doi.org/10.1137/S1540345902411402. 3

[64] L. SCOTT AND M. VOGELIUS, *Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials*, ESAIM: M2AN, 19 (1985), pp. 111-143, https://doi.org/10.1051/m2an/1985190101111. 30

[65] M. F. WHEELER, G. XUE, AND I. YOTOV, *A multiscale mortar multipoint flux mixed finite element method*, Math. Models Methods Appl. Sci., 46 (2012), pp. 759-796, https://doi.org/10.1051/m2an/2011064. 3

[66] J. XU AND L. ZIKATANOV, *Some observations on Babuška and Brezzi theories*, Numer. Math., 94 (2003), pp. 195-202, https://doi.org/10.1007/s002110100308. 28
[67] S. ZHANG, *A new family of stable mixed finite elements for the 3D Stokes equations*, Math. Comp., 250 (2005), pp. 543–554, https://doi.org/10.1090/S0025-5718-04-01711-9_30

(R.A.) DEPARTAMENTO DE INGENIERÍA MATEMÁTICA & CI$^2$MA, UNIVERSIDAD DE CONCEPCIÓN, CASILLA 160-C, CONCEPCIÓN, CHILE. rodolfo.araya@udec.cl

(C.H.) METROPOLITAN STATE UNIVERSITY OF DENVER, DEPARTMENT OF MATHEMATICAL AND COMPUTER SCIENCES, P.O. Box 173362, CAMPUS BOX 38, DENVER, CO 80217-3362, harderc@msudenver.edu

(A.H.P.) DEPARTAMENTO DE MATEMÁTICA Y FÍSICA APLICADAS, UNIVERSIDAD CATÓLICA DE LA SANTÍSIMA CONCEPCIÓN, CASILLA 297, CONCEPCIÓN, CHILE, AND, GRUPO DE INVESTIGACIÓN EN ANÁLISIS NUMÉRICO Y CÁLCULO CIENTÍFICO, GIANuC$^2$, CONCEPCIÓN, CHILE, apoza@ucsc.cl

(F.V.) DEPARTMENT OF COMPUTATIONAL AND MATHEMATICAL METHODS, NATIONAL LABORATORY FOR SCIENTIFIC COMPUTING - LNCC, Av. GETÚLIO VARGAS, 333, 25651-070 PETRÓPOLIS - RJ, BRAZIL, valentin@lncc.br