# Revisiting the robustness of the multiscale hybrid-mixed method: the face-based strategy

Diego Paredes, Frederic Valentin, Henrique M Versieux

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

## To cite this version:

Diego Paredes, Frederic Valentin, Henrique M Versieux. Revisiting the robustness of the multiscale hybrid-mixed method: the face-based strategy. 2023. hal-03878574v2

HAL Id: hal-03878574

https://inria.hal.science/hal-03878574v2

Preprint submitted on 6 Jun 2023

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
REVISITING THE ROBUSTNESS OF THE MULTISCALE
HYBRID-MIXED METHOD: THE FACE-BASED STRATEGY

DIEGO PAREDES*, FRÉDÉRIC VALENTIN†, HENRIQUE M. VERSIEUX‡

ABSTRACT. This work proposes a new finite element for the mixed multiscale hybrid method (MHM) applied to the Poisson equation with highly oscillatory coefficients. Unlike the original MHM method, multiscale bases are the solution to local Neumann problems driven by piecewise continuous polynomial interpolation on the skeleton faces of the macroscale mesh. As a result, we prove the optimal convergence of MHM by refining the face partition and leaving the mesh of macroelements fixed. This property allows the MHM method to be resonance free under the usual assumptions of local regularity. The numerical analysis of the method also revisits and complements the original approach proposed by D. Paredes, F. Valentin and H. Versieux (2017). Numerical experiments assess the new theoretical results.

# 1. INTRODUCTION

Multiscale finite element methods approximate the solution of differential equations with heterogeneous coefficients on coarse partitions, becoming an attractive alternative to classical finite element methods requiring very fine meshes. Since the seminal work [4], there has been a vast literature on the subject, such as the variational multiscale method (VMS) [18], the multiscale finite element method (MsFEM) and its generalization (GMsFEM) [10], the (Petrov-)Galerkin enriched method (GEM and PGEM) [1, 12], the heterogeneous multiscale method (HMM) [27], multiscale mortar method [3], the local orthogonal decomposition (LOD) method [21], the hybrid localized spectral decomposition (LSD) method [20], and the hybrid higher order multiscale (MsHHO) method [8], to cite a few. The main idea of multiscale methods is replacing polynomial approximation spaces with space spam by basis function computed from independent partial differential problems at the element level. As a result, the multiscale methods become precise on coarse meshes while parallel facilities handle the overhead computation due to basis functions computations.

<footnote>The first author was supported by Project EOLIS (MATH-AMSUD 21-MATH-04) and ANID-Chile through grant FONDECYT-1181572. The second 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. The third author was partially supported by CNPq/Brazil No. 307873/2019-6.</footnote>
The MHM method is an example of a multiscale finite element method that was originally proposed in [15, 2] (see [16] for an abstract framework). It is a by-product of a hybrid formulation that starts at the continuous level placed in a partition, which characterizes the exact solution in local and global contributions. When discretized, decouples the local and the global problems. The global formulation is responsible for the degrees of freedom over the skeleton of the coarse partition, and the local problems provide the multiscale basis functions. The original MHM method relates to other multiscale methods like the MsHHO method (c.f. [6]), for example, and has variants that extend it to handle polygonal meshes [5] and preserve conservation properties locally [9].

The a priori error analysis provided in the seminal work [2] showed that the MHM method is optimally convergent with respect to the discretization parameter $H$. These error estimates are useful in the $H \le \varepsilon$ regime, where $\varepsilon$ is the oscillatory frequency. Then, the numerical analysis of the MHM method was extended to the asymptotic regime of $H \gg \varepsilon$ in [23] under the assumption that equation coefficients are $\varepsilon$- periodic. For that, [23] considered discontinuous polynomial interpolation for the Lagrange multiplier on the faces of $\mathcal{T}_H$ a partition of $\Omega$. This led to the following estimate in the simplest interpolation case (constant on faces)

$$ (1) \qquad \|u^{\varepsilon} - u_{H}^{\varepsilon}\|_{1, \mathcal{T}_{H}} = O\left(\left(\frac{\varepsilon}{H}\right)^{1/2} + H\right), $$

where $\|\cdot\|_{1, \mathcal{T}_H}$ is the $H^1$-broken norm, and the exact and approximate solutions are $u^\varepsilon$ and $u_H^\varepsilon$. Note that the estimate (1) differs from that claimed in [23] as it includes the resonance error $O((\varepsilon/H)^{1/2})$. So, the first objective of the present work is to revisit the analysis in [23] and demonstrate that the main error is, in fact, given by (1). This is done by correcting the proposed asymptotic expansion estimate in [23, Theorem 1]. In the process, we fill in the gap in [23, Theorem 1] by providing the details of the proof in the case of local Neumann-type problems. Here, we also assume that the coefficients are periodic for the numerical analysis of the method, but the proposed method is ready to approximate the solution of problems with more general coefficients. In the case of the stochastic coefficient, for example, we can use the MHM method to calculate the numerical approximation of each realization of the random coefficients.

Motivated by (1), the present work proposes a new multiscale finite element for the Poisson equation with oscillatory coefficients within the class of MHM methods. It extends the numerical analysis given in [23] to prove that the solution of the MHM method converges when we adopt piecewise continuous or discontinuous polynomial interpolation on faces to approximate the Lagrange multiplier. Notably, we prove that the MHM method, with
(dis)continuous polynomial interpolation on faces, yields an error estimate of the form

$$ (2) \qquad \|u^{\varepsilon} - u_{H}^{\varepsilon}\|_{1, \mathcal{T}_{H}} = O\left(\left(\frac{\varepsilon}{H}\right)^{1/2} + h\right), $$

where $h$ is the diameter of the face partitions. Interestingly, the estimate (2) indicates that the MHM method is indeed *without resonance* as claimed in [23] if the space-based refinement strategy is adopted (e.g. $h \to 0$ for a fixed $H$). The estimate was numerically anticipated in [23] in the discontinuous interpolation case. As for the continuous interpolation choice, a numerical test confirms the theoretical findings in Section 5.

It is also interesting to compare the MHM method with more classic approaches. For example, the present method computes the element-by-element multiscale basis as in the classic MsFEM [17] but differs in that MsFEM relies on local Dirichlet problems to drive the basis. The error estimates in the conforming version of MsFEM behave as in (1), and it is proven that the numerical resonance error (e.g. $(\varepsilon/H)^{1/2}$)) considerably diminishes when adopting an oversampling technique that extends the computational domain and makes the method nonconforming in $H^1(\Omega)$. The present method is also nonconforming and resonance-free but has the advantage that the multiscale basis functions are still computed element by element. On the other hand, MsFEM leads to a definite positive linear system, while the original MHM method has a saddle point structure. A coercive version of the MHM method was recently proposed in [20] by judiciously decomposing the Lagrange multiplier space used to enforce the weak continuity of the approximate solution. Other coercive and resonance-free multiscale methods exist, achieved by base functions computed on patches of elements, such as the LOD method.

The outline of this work is the following: Section 2 presents the model and notations and revisits the asymptotic results first presented in [23]. Section 3 defines the MHM method, and Section 4 provides the error analysis. The numerics are in Section 5 and conclusions in Section 6. The proof of the asymptotic result correcting [23][Theorem 1] is in the appendix.

## 2. THE MODEL AND ASYMPTOTIC EXPANSION

2.1. **The model.** Let $\Omega \subset \mathbb{R}^d$, $d \in \{2, 3\}$, be an open and bounded domain, with Lipschitz boundary $\partial\Omega := \partial\Omega_D \cup \partial\Omega_N$. The elliptic problem consists of finding $u_{\varepsilon}$ the solution of

$$ (3) \qquad \begin{aligned} -\nabla \cdot (A_{\varepsilon} \nabla u_{\varepsilon}) &:= -\frac{\partial}{\partial x_i} (a_{ij}(x/\varepsilon) \frac{\partial}{\partial x_j} u_{\varepsilon}) = f \quad \text{in } \Omega, \\ u_{\varepsilon} &= g \quad \text{on } \partial\Omega_D \quad \text{and} \quad A_{\varepsilon} \nabla u \cdot n = b \quad \text{on } \partial\Omega_N, \end{aligned} $$
where $\mathcal{A}_{\varepsilon}(\mathbf{x}) = \mathcal{A}(\mathbf{x}/\varepsilon) = (a_{ij}(\mathbf{x}/\varepsilon))$ is a symmetric positive definite matrix. Here $\varepsilon \in (0, 1)$ is the (small) parameter controlling the fine scale oscillations of the physical coefficient, $g \in H^{1/2}(\partial\Omega_D)$, $\mathbf{n}$ represents the unit outward normal vector on $\partial\Omega$, $b \in H^{-1/2}(\partial\Omega_N)$, and $f \in L^2(\Omega)$ (these spaces having their usual meaning), and $\mathbf{x} := (x_j)$ is a typical point in $\Omega$. Above, and throughout the paper, the indices $i, j$ run from $1, \dots, d$, even when not explicitly mentioned, and we employ the Einstein summation convention, i.e., repeated indices indicate summation.

We also assume that $a_{ij} \in L_{\text{per}}^{\infty}(Y)$, i.e., $a_{ij} \in L^{\infty}(\mathbb{R}^d)$ and it is $Y$-periodic, $Y := (0, 1)^d$, and there exist positive constants $\gamma_a$ and $\gamma_b$ such that

$$ (4) \qquad \gamma_a |\boldsymbol{\xi}|^2 \le a_{ij}(\boldsymbol{y}) \xi_i \xi_j \le \gamma_b |\boldsymbol{\xi}|^2 \quad \text{for all } \boldsymbol{\xi} := \{\xi_i\} \in \mathbb{R}^d \text{ and } \boldsymbol{y} \in Y, $$

where $|\cdot|$ represents the Euclidean norm. In the case $\partial\Omega_D = \emptyset$, we also assume that the following compatibility condition holds

$$ (5) \qquad \int_{\Omega} f \, d\mathbf{x} = \int_{\partial\Omega} b \, ds. $$

2.2. **Broken spaces and norms.** We note by $d_{\Omega}$ the side of the maximum square (or cube in 3D) contained in $\Omega$, and $L_{\Omega}$ is the length of $\partial\Omega$ in the 2D case (or area in 3D). Let $B \subset \mathbb{R}^d$ be an open set. We define

$$ \|v\|_{m,\infty,B} := \max_{|\alpha| \le m} \{\text{ess. sup}_{\mathbf{x} \in B} |\partial^\alpha v(\mathbf{x})|\} \quad \text{and} \quad |v|_{m,\infty,B} := \max_{|\alpha|=m} \{\text{ess. sup}_{\mathbf{x} \in B} |\partial^\alpha v(\mathbf{x})|\}, \\ \text{and for } 1 \le q < \infty $$

$$ \|v\|_{m,q,B} := \left( \int_B \sum_{|\alpha| \le m} |D^\alpha v|^q d\mathbf{x} \right)^{1/q} \quad \text{and} \quad |v|_{m,q,B} := \left( \int_B \sum_{|\alpha|=m} |D^\alpha v|^q d\mathbf{x} \right)^{1/q}. $$

We define the broken (semi-)norms related to a partition $\mathcal{T}_H$ of $B$ in elements $K$, with diameter $H > 0$, by

$$ \|v\|_{m,q,\mathcal{T}_H} := \left( \sum_{K \in \mathcal{T}_H} \|v\|_{m,q,K}^2 \right)^{1/2} \quad \text{and} \quad |v|_{m,q,\mathcal{T}_H} := \left( \sum_{K \in \mathcal{T}_H} |v|_{m,q,K}^2 \right)^{1/2}, $$

and the norm in the $H(\text{div};B)$ space, i.e., the space of functions belonging to $(L^2(B))^d$ with divergence also in $L^2(B)$, by

$$ \|\boldsymbol{\sigma}\|_{\text{div},B} := \left( \int_B |\boldsymbol{\sigma}|^2 d\mathbf{x} + \int_B |\nabla \cdot \boldsymbol{\sigma}|^2 d\mathbf{x} \right)^{1/2}. $$
Also, we simplify the notation with respect to the norms of a vector function $\chi$ with components $\chi^j$ by setting

$$ \|\chi\|_{s,p,Y} := \max_{1 \le j \le d} \|\chi^j\|_{s,p,Y}. $$

Hereafter, we do not make reference to the domain $B$, or to the coefficient $q$ when $B = \Omega$, or $q = 2$, respectively. In what follows, $c$ denotes a generic constant independent of $\varepsilon$ and $H$, although it may change in each occurrence.

2.3. **Asymptotics of $u_\varepsilon$**. The convergence analysis of the method is based on an asymptotic expansion of $u_\varepsilon$. For that, we assume throughout this work that

$$ (6) \qquad d_\Omega > c\varepsilon, $$

where the constant $c > 1$ is of order one. Here, we focus on the first order asymptotic expansion approximation, and then, we recall that $u_\varepsilon$ is approximated by

$$ (7) \qquad u_\varepsilon^1(\mathbf{x}) := u_0(\mathbf{x}) + \varepsilon \sum_{j=1}^n \chi^j(\mathbf{x}/\varepsilon) \partial_{x_j} u_0(\mathbf{x}), $$

where $\chi^j$ are the solutions of cell problems. More precisely, $\chi^j \in H_{\text{per}}^1(Y)$, i.e., $\chi^j \in H_{\text{loc}}^1(\mathbb{R}^d)$ and $Y$-periodic, is the weak solution with zero mean value on $Y$ of

$$ (8) \qquad \nabla_{\mathbf{y}} \cdot \mathcal{A}(\mathbf{y}) \nabla_{\mathbf{y}} \chi^j = \nabla_{\mathbf{y}} \cdot \mathcal{A}(\mathbf{y}) \nabla_{\mathbf{y}} y_j = \frac{\partial}{\partial y_i} a_{ij}(\mathbf{y}), $$

where $\mathbf{y} := \mathbf{x}/\varepsilon$ with coordinates $(y_j)$. Denote $\mathcal{A}_0$ the symmetric positive definite matrix

$$ (9) \qquad \mathcal{A}_0 := (a_{ij}^0), \quad a_{ij}^0 := \frac{1}{|Y|} \int_Y a_{lm}(\mathbf{y}) \frac{\partial}{\partial y_l} (y_i - \chi^i) \frac{\partial}{\partial y_m} (y_j - \chi^j) d\mathbf{y}. $$

Then, the function $u_0 \in H^1(\Omega)$ is the weak solution of the following homogenized problem (see [23, 19] for more details)

$$ (10) \qquad \begin{gathered} -\nabla \cdot (\mathcal{A}_0 \nabla u_0) = f \quad \text{in } \Omega, \\ u_0 = g \quad \text{on } \partial\Omega_D, \quad \mathcal{A}_0 \nabla u_0 \cdot \mathbf{n} = b \quad \text{on } \partial\Omega_N. \end{gathered} $$

The next theorem addresses the error estimate between $u_\varepsilon$ and $u_\varepsilon^1$ and the fundamental question of the dependence of the constants in terms of the domain $\Omega$. Its proof is postponed to the appendix.

**Theorem 1** (Asymptotic convergence). Let $u_\varepsilon$ be the solution of problem (3), and $u_0$ and $u_\varepsilon^1$ be defined by equations (10) and (7), respectively. Assume the coefficients $a_{ij}$ are $Y$-periodic, (6) holds and

$$ (11) \qquad u_0 \in H^2(\Omega) \quad \text{and} \quad \chi^j \in W_{\text{per}}^{1,q}(Y) \quad \text{with} \ q > d. $$
Then,

$$ (12) \qquad \|u_{\varepsilon} - u_{\varepsilon}^{1}\|_{1} \le c(p') \left[ \varepsilon + \left( \frac{L_{\Omega} \varepsilon}{|\Omega|} \right)^{1/2-1/p'} \right] (1+L_{\Omega}) \|\mathbf{x}\|_{1,q,Y} \|u_{0}\|_{2} $$

where $2 < p' < K_d$ with

$$ (13) \qquad \begin{cases} K_d = \infty, & \text{if } d=2, \\ K_d = 2d/(d-2) & \text{if } d > 2. \end{cases} $$

Also, the constant $c(p')$ depends on $p'$, and $c(p') \to \infty$ when $p' \to K_d$, and they may depend on the cone property of $\Omega$, but they do not depend on the size of $\Omega$ when $\Omega$ is a convex domain.

**Remark 2** (Revising Theorem 1 in [23]). The Theorem 1 above corrects [23][Theorem 1], in the sense that (12) contains a dependency of $\Omega$. Notably, we get

$$ c(p') \left( \frac{L_{\Omega} \varepsilon}{|\Omega|} \right)^{1/2-1/p'} \quad \text{rather than} \quad c(p') (L_{\Omega} \varepsilon)^{1/2-1/p'}. $$

Furthermore, the proof of the Theorem 1 (see appendix) details the estimation for Neumann problems that complements the results in [23]. □

### 3. THE MHM METHOD

The MHM method introduced in [2] is revisited in this section in the context of continuous or discontinuous polynomial interpolation on faces.

**3.1. Partitions and broken finite spaces.** Let $\mathcal{T}_H$ be a regular triangulation of $\Omega$ characterised by $H := \max_{K \in \mathcal{T}_H} h_K$, where given element $K \in \mathcal{T}_H$, $h_K = \text{diam}(K)$. 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 by $\mathcal{E}$. To each $E \in \mathcal{E}$, a normal $\mathbf{n}$ is associated, taking care to ensure this is directed outward on $\partial\Omega$. For each $K \in \mathcal{T}_H$, $\mathbf{n}^K$ indicates the outward normal on $\partial K$, and we adopt the notation $\mathbf{n}_E^K := \mathbf{n}^K|_E$ for each $E \subset \partial K$. We also introduce $\{\mathcal{E}_h\}_{h>0}$ a family of regular partition of $\mathcal{E}$ wherein each $E \subset \mathcal{E}$ is decomposed into sub-faces $F$ of length $h_F$, and we collect the faces of $F \subset \partial K$ in $\mathcal{E}_h^K$. We set $h := \max_{F \in \mathcal{E}_h} h_F$. We emphasise that each face $E \subset \mathcal{E}$ has its own family of partitions independent of one another.

Consider the decomposition

$$ V := H^1(\mathcal{T}_H) = V_0 \oplus V_0^{\perp}, $$
where $V_0$ is defined by

$$V_0 := \{v \in L^2(\Omega) : v_K \text{ is constant on } K \in \mathcal{T}_H\},$$

where $v_D := v|_D$ and $D$ is any measurable set, and $V_0^\perp$ is the orthogonal complement of $V_0$ with respect to the $L^2$ product, i.e., functions $v$ in $V$ such that $v_K \in L_0^2(K)$ (functions in $L^2(K)$ with mean value zero on $K$). Also, set

$$(14) \qquad \Lambda := \{ \boldsymbol{\sigma}_K \cdot \boldsymbol{n}^K |_{\partial K} : \boldsymbol{\sigma} \in H(\mathrm{div}; \Omega), \forall K \in \mathcal{T}_H \}.$$

With the notation

$$(w, v)_{\mathcal{T}_H} := \sum_{K \in \mathcal{T}_H} (w_K, v_K)_K,$$

for all $w, v \in L^2(D)$, where $(\cdot, \cdot)_D$ stands for the $L^2$ inner product over $D$, we equip $L^2(D)$ with the induced norm $\| \cdot \|_{0,D}$. Also, we define the broken gradient operator $\nabla_H : V \to L^2(\Omega)$ as such, for all $v \in V$,

$$(\nabla_H v)|_K := \nabla v_K \quad \text{for all } K \in \mathcal{T}_H.$$

With such a notation, $V$ space is equipped with the norm induced by the inner product

$$(15) \qquad (w, v)_V := d_{\Omega}^{-2}(w, v)_{\mathcal{T}_H} + (\nabla_H w, \nabla_H v)_{\mathcal{T}_H},$$

where $d_\Omega$ is the diameter of $\Omega$. We define the following norm in the broken space $V$

$$(16) \qquad \|v\|_V^2 := \sum_{K \in \mathcal{T}_H} d_{\Omega}^{-2} \|v_K\|_{0,K}^2 + \|\nabla v_K\|_{0,K}^2 \quad \text{for all } v \in V.$$

We equip the spaces $H(\mathrm{div}; \Omega)$ and $\Lambda$ with the norms,

$$(17) \qquad \|\boldsymbol{\sigma}\|_{\mathrm{div},\Omega}^2 := \sum_{K \in \mathcal{T}_H} (\|\boldsymbol{\sigma}_K\|_{0,K}^2 + d_{\Omega}^2 \|\nabla \cdot \boldsymbol{\sigma}_K\|_{0,K}^2),$$

$$(18) \qquad \|\mu\|_{\Lambda} := \inf_{\substack{\boldsymbol{\sigma} \in H(\mathrm{div};\Omega) \\ \boldsymbol{\sigma}_K \cdot \boldsymbol{n}_K = \mu \text{ on } \partial K, K \in \mathcal{T}_H}} \|\boldsymbol{\sigma}\|_{\mathrm{div},\Omega}.$$

The duality pairing between $H^{-1/2}(\partial K)^d$ and $H^{1/2}(\partial K)^d$ is $(\mu_K, v_K)_{\partial K}$, and we define

$$(\mu, v)_{\partial \mathcal{T}_H} := \sum_{K \in \mathcal{T}_H} (\mu_K, v_K)_{\partial K}.$$

Observe that the following equality holds (c.f. [14])

$$(19) \qquad \|\mu\|_{\Lambda} = \sup_{v \in V} \frac{(\mu, v)_{\partial \mathcal{T}_H}}{\|v\|_V} \quad \text{for all } \mu \in \Lambda.$$

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. We define
three local mappings $T_K^\epsilon : H^{-1/2}(\partial K) \to V_0^\perp(K)$, $\hat{T}_K^\epsilon : L^2(K) \to V_0^\perp(K)$ and $\bar{T}_K^\epsilon : H^{-1/2}(\partial K) \to V_0^\perp(K)$

$$ (20) \qquad (\mathcal{A}_\epsilon \nabla T_K^\epsilon \mu, \nabla v)_K = (\mu, v)_{\partial K} \quad \text{for all } v \in H^1(K) \cap L_0^2(K), $$

$$ (20) \qquad (\mathcal{A}_\epsilon \nabla \hat{T}_K^\epsilon q, \nabla v)_K = (q, v)_K \quad \text{for all } v \in H^1(K) \cap L_0^2(K). $$

$$ (\mathcal{A}_\epsilon \nabla \bar{T}_K^\epsilon \xi, \nabla v)_K = (\xi, v)_{\partial K \cap \partial \Omega_N} \quad \text{for all } v \in H^1(K) \cap L_0^2(K). $$

The corresponding global mappings $T_\epsilon \in \mathcal{L}(\Lambda, V_0^\perp)$, $\hat{T}_\epsilon \in \mathcal{L}(L^2(\Omega), V_0^\perp)$, and $\bar{T}_\epsilon \in \mathcal{L}(L^2(\Omega), V_0^\perp)$ are such that

$$ (21) \qquad T_\epsilon \mu|_K := T_K^\epsilon \mu_{\partial K} , \quad \hat{T}_\epsilon q|_K := \hat{T}_K^\epsilon q_K, \quad \text{and} \quad \bar{T}_\epsilon \xi|_K := \bar{T}_K^\epsilon \xi_{\partial K}, $$

for all $\mu, \xi \in \Lambda$ and $q \in L^2(\Omega)$, then they are well-defined and bounded (c.f. [2]). We also define the homogenized counterpart of the local mappings $T_K^\epsilon$ and $\hat{T}_K^\epsilon$, denoted by $T_K^0$, $\hat{T}_K^0$ and $\bar{T}_K^0$ wherein tensor $\mathcal{A}_\epsilon$ is replaced by $\mathcal{A}_0$, and their global counterpart by $T_0$, $\hat{T}_0$ and $\bar{T}_0$.

3.2. **The method.** The weak exact solution $u^\epsilon$ of (3) is equivalent to the solution of the hybrid formulation (c.f. [24]): Find $(u_\epsilon, \lambda_\epsilon) \in V \times \Lambda$ such that

$$ (22) \qquad \begin{aligned} (\mathcal{A}_\epsilon \nabla u_\epsilon, \nabla v)_{\mathcal{T}_H} - (\lambda_\epsilon, v)_{\partial \mathcal{T}_H} &= (f, v)_{\mathcal{T}_H} + (b, v)_{\partial \Omega_N} && \text{for all } v \in V, \\ (\mu, u_\epsilon)_{\partial \mathcal{T}_H} &= (\mu, g)_{\partial \Omega_D} && \text{for all } \mu \in \Lambda. \end{aligned} $$

In addition, we have that $u^\epsilon$ can be characterized as follows (c.f. [2])

$$ (23) \qquad u_\epsilon = u_\epsilon^0 + T_\epsilon \lambda_\epsilon + \hat{T}_\epsilon f + \bar{T}_\epsilon b, $$

where $u_\epsilon^0|_K = \frac{1}{|K|} \int_K u_\epsilon dx$ for all $K \in \mathcal{T}_H$ and $\lambda_\epsilon$ satisfies (22).

We set the finite element space $\Lambda_h^\ell \subset \Lambda \cap L^2(\mathcal{E}_h)$ as the space of (dis)continuous piecewise polynomials on each face $E$ of elements $K \in \mathcal{T}_H$, i.e.,

$$ (24) \qquad \Lambda_h^\ell := \{ \mu \in \Lambda \cap L^2(\mathcal{E}_h) : \mu_E \text{ is either continuous or discontinuous and} \\ \qquad \qquad \mu_F \in \mathbb{P}_\ell(F), \text{ for all } F \in \mathcal{E}_h \}, $$

where $\mathbb{P}_\ell(F)$ is the space of polynomials up degree $\ell \ge 1$ if $\Lambda_h^\ell$ is the continuous space and $\ell \ge 0$ otherwise, on $F$.

The MHM method corresponds to find $(u_\epsilon^{0,H}, \lambda_\epsilon^h) \in V_0 \times \Lambda_h^\ell$ such that

$$ (25) \qquad \begin{aligned} (\lambda_\epsilon^h, v^0)_{\partial \mathcal{T}_H} &= -(f, v^0)_{\mathcal{T}_H} - (b, v^0)_{\partial \Omega_N} && \text{for all } v^0 \in V_0, \\ (\mu^h, u_\epsilon^{0,H} + T_\epsilon \lambda_\epsilon^h)_{\partial \mathcal{T}_H} &= -(\mu^h, \hat{T}_\epsilon f)_{\partial \mathcal{T}_H} - (\mu^h, \bar{T}_\epsilon b)_{\partial \mathcal{T}_H} + (\mu^h, g)_{\partial \Omega_D} && \text{for all } \mu^h \in \Lambda_h^\ell, \end{aligned} $$
and construct the approximation $u_{\varepsilon}^h$ of $u_{\varepsilon}$ as follows

$$ (26) \qquad u_{\varepsilon} \approx u_{\varepsilon}^{h} := u_{\varepsilon}^{0,H} + T_{\varepsilon}\lambda_{\varepsilon}^{h} + \hat{T}_{\varepsilon}f + \overline{T}_{\varepsilon}b \notin H^{1}(\Omega), $$

where $T_{\varepsilon}$, $\hat{T}_{\varepsilon}$ and $\overline{T}_{\varepsilon}$ are given in (21) and (20), and $\sigma_{\varepsilon} := -\mathcal{A}_{\varepsilon}\nabla u_{\varepsilon}$ is approximated as follows

$$ (27) \qquad \sigma_{\varepsilon} \approx \sigma_{\varepsilon}^{h} := -\mathcal{A}_{\varepsilon}\nabla_{H}u_{\varepsilon}^{h} \in H(\operatorname{div}; \Omega). $$

Discretization decouples the local problems (20) from the global one (25). Thereby, a staggered algorithm can be adopted to solve the system. To see this more clearly, it is instructive to consider $T \lambda_{\varepsilon}^h$ in more detail. Suppose $\{\psi_i\}_{i=1}^{\dim \Lambda_h^\ell}$ is a basis for $\Lambda_h^\ell$, and define the set $\{\eta_i\}_{i=1}^{\dim \Lambda_h^\ell} \subset V_0^\perp$ such that, $\eta_i = T_{\varepsilon}\psi_i$. Now, giving $\lambda_{\varepsilon}^h = \sum_{i=1}^{\dim \Lambda_h^\ell} c_i \psi_i$ in $\Lambda_h^\ell$, $c_i \in \mathbb{R}$, the linearity of the operator $T_{\varepsilon}$ implies we may uniquely write

$$ T_{\varepsilon} \lambda^{h} = \sum_{i=1}^{\dim \Lambda_{h}^{\ell}} c_{i} T_{\varepsilon} \psi_{i} = \sum_{i=1}^{\dim \Lambda_{h}^{\ell}} c_{i} \eta_{i}. $$

Therefore, the degrees of freedom $c_i$'s of $\lambda_{\varepsilon}^h$ are “inherited” by $T_{\varepsilon} \lambda_{\varepsilon}^h$. It then follows from (26) that

$$ (28) \qquad u_{\varepsilon}^{h} = u_{\varepsilon}^{0,h} + \sum_{i=1}^{\dim \Lambda_{h}^{\ell}} c_{i} \eta_{i} + \hat{T}_{\varepsilon} f + \overline{T}_{\varepsilon} b. $$

As a result, the global formulation (25) is responsible for computing the degrees of freedom of $u_{\varepsilon}^{0,h}$ (one per element) and the $c_i$'s in (28), once the multiscale basis functions $\eta_i$'s, $\hat{T}_{\varepsilon}f$ and $\overline{T}_{\varepsilon}b$ are available from the local problems. Also, it is interesting to note that heterogeneous and/or high-contrast aspects of the media automatically impact the design of the basis functions $\eta_i$'s (see Figure 1 for an example) as well as $\hat{T}_{\varepsilon}f$ and $\overline{T}_{\varepsilon}b$.

# 4. NUMERICAL ANALYSIS

## 4.1. Existence and uniqueness, and best approximation.
The proof of well-posedness of (25) follows closely [2]. For completeness, we revisit it here assuming $\Lambda_h^\ell$ is given in (24), and we present a best approximation result on $\Lambda_h^\ell$ which drives convergence of $u_h^\varepsilon$ to $u^\varepsilon$. The proof differs from the one originally proposed in [2].
**Lemma 3** (Well-posedness and best approximation). There exist positive constants $\alpha_0$ and $\beta_0$, independent of mesh parameters, such that

$$ (29) \qquad \sup_{\mu^h \in \Lambda_h^\ell} \frac{(\mu^h, v_0)_{\partial \mathcal{T}_H}}{\|\mu^h\|_{\Lambda}} \geq \beta_0 \|v_0\|_V \quad \text{for all } v_0 \in V_0, $$

$$ (30) \qquad (\mu^h, T_\varepsilon \mu^h)_{\partial \mathcal{T}_H} \geq \alpha_0 \|\mu^h\|_{\Lambda}^2 \quad \text{for all } \mu^h \in \mathcal{N}_h^\ell, $$

where $\mathcal{N}_h^\ell$ is

$$ (31) \qquad \mathcal{N}_h^\ell := \{ \mu^h \in \Lambda_h^\ell : (\mu^h, v_0)_{\partial \mathcal{T}_H} = 0 \quad \text{for all } v_0 \in V_0 \}. $$

Hence, (25) is well-posed, and there exists $c$ such that

$$ (32) \qquad \|u_{\varepsilon} - u_{\varepsilon}^h\|_V \leq c \inf_{\mu^h \in \Lambda_h^\ell} \|\lambda_{\varepsilon} - \mu^h\|_{\Lambda}, $$

where $u_\varepsilon$ and $u_\varepsilon^h$ are given in (23) and (26), respectively.

*Proof.* Let us consider the space

$$ (33) \qquad X_H := \{ \tau \in H(\text{div}; \Omega) : \tau_K \in RT_0(K) \quad \text{for all } K \in \mathcal{T}_H \}, $$

where $RT_0(K)$ stands for the lowest order Raviart-Thomas space in $K$. Then, given $v_0 \in V_0$, there exists $\tau \in X_H$ such that $\nabla \cdot \tau = v_0$ in $\Omega$ and $\beta_0 \|\tau\|_{\text{div},\Omega} \le \|v_0\|_{0,\Omega} = \|v_0\|_V$, where $\beta_0$ does not depend on $\mathcal{T}_H$, $H$, or $h$. Using the definition of the norm in $\Lambda_h^\ell$, we pick $\mu_E^h := \tau \cdot n_E^K \in \Lambda_h^\ell$ on every $E \in \mathcal{E}$

$$ \beta_0 \|\mu^h\|_{\Lambda} \|v_0\|_V \le \beta_0 \|\tau\|_{\text{div},\Omega} \|v_0\|_V \le \int_{\Omega} \nabla \cdot \tau v_0 \,dx = (\mu^h, v_0)_{\partial \mathcal{T}_H}, $$

which proves (29). Next, given $\mu^h \in \mathcal{N}_h^\ell$, it holds

$$ \begin{aligned} \|\mu^h\|_{\Lambda} &= \sup_{v \in V} \frac{(\mu^h, v)_{\partial \mathcal{T}_H}}{\|v\|_V} \le \sup_{v_0^{\perp} \in V_0^{\perp}} \frac{(\mu^h, v_0^{\perp})_{\partial \mathcal{T}_H}}{\|v_0^{\perp}\|_V} = \sup_{v_0^{\perp} \in V_0^{\perp}} \frac{(A_{\varepsilon} \nabla_H T_{\varepsilon} \mu^h, \nabla_H v_0^{\perp})_{\partial \mathcal{T}_H}}{\|v_0^{\perp}\|_V} \\ &\le \gamma_b^{1/2} \|A_{\varepsilon}\|_{\partial \mathcal{T}_H}^{1/2} \|T_{\varepsilon} \mu^h\|_{0, \partial \mathcal{T}_H} \\ &= \gamma_b^{1/2} (\mu^h, T_{\varepsilon} \mu^h)_{\partial \mathcal{T}_H}^{1/2} \end{aligned} $$

and (30) follows. Well-posedness follows from standard saddle-point theory. As for the best approximation result (32), first observe from the second equations in (22) and (25) that

$$ (\mu^h, u_{\varepsilon} - u_{\varepsilon}^h)_{\partial \mathcal{T}_H} = 0 \quad \text{for all } \mu^h \in \Lambda_h^{\ell}. $$
Then, using that $(\lambda - \lambda_{\varepsilon}^{h}, v_{0})_{\partial\mathcal{T}_{H}} = 0$ for all $v_{0} \in V_{0}$ from the first equations in (22) and (25), it holds

$$
\begin{aligned}
|u_{\varepsilon} - u_{\varepsilon}^{h}|_{1, \mathcal{T}_{H}}^{2} &= |T_{\varepsilon}(\lambda - \lambda_{\varepsilon}^{h})|_{1, \mathcal{T}_{H}}^{2} \leq c \|\mathcal{A}_{\varepsilon}^{1/2} \nabla_{H}(T_{\varepsilon}(\lambda - \lambda_{\varepsilon}^{h}))\|_{0, \mathcal{T}_{H}}^{2} \\
&= c (\lambda - \lambda_{\varepsilon}^{h}, T_{\varepsilon}(\lambda - \lambda_{\varepsilon}^{h}))_{\partial \mathcal{T}_{H}} \\
&= c (\lambda - \lambda_{\varepsilon}^{h}, u^{\varepsilon} - u_{h}^{\varepsilon})_{\partial \mathcal{T}_{H}} \\
&= c (\lambda - \mu^{h}, u^{\varepsilon} - u_{h}^{\varepsilon})_{\partial \mathcal{T}_{H}} \\
&\leq c \|\lambda - \mu^{h}\|_{\Lambda} \|u_{\varepsilon} - u_{\varepsilon}^{h}\|_{V}.
\end{aligned}
$$

It remains to prove $\|u_{\varepsilon} - u_{\varepsilon}^{h}\|_{0, \mathcal{T}_{H}} \leq c|u_{\varepsilon} - u_{\varepsilon}^{h}|_{1, H}$. To see this, first notice that from the definition of $u_{\varepsilon}$ and $u_{\varepsilon}^{h}$ we have

$$
\begin{equation} \tag{34}
\begin{aligned}
\|u_{\varepsilon} - u_{\varepsilon}^{h}\|_{0, \mathcal{T}_{H}} &\le \|u_{\varepsilon}^{0} + T_{\varepsilon}\lambda_{\varepsilon} - u_{\varepsilon}^{0,H} - T_{\varepsilon}\lambda_{\varepsilon}^{h}\|_{0, \mathcal{T}_{H}} \\
&\le \|u_{\varepsilon}^{0} - u_{\varepsilon}^{0,H}\|_{0, \mathcal{T}_{H}} + c H |T_{\varepsilon}\lambda_{\varepsilon} - T_{\varepsilon}\lambda_{\varepsilon}^{h}|_{1, \mathcal{T}_{H}},
\end{aligned}
\end{equation}
$$

where we used the triangle inequality, the Poincaré inequality, and the assumption on the regularity of the mesh. To estimate $\|u_{\varepsilon}^{0} - u_{\varepsilon}^{0,H}\|_{0, \mathcal{T}_{H}}$ we follow [23], which we revisit here for completeness. Take $\boldsymbol{\sigma}^{*}$ be the vector-valued function belonging to the lowest order Raviart-Thomas space such that $\nabla \cdot \boldsymbol{\sigma}^{*} = u_{\varepsilon}^{0} - u_{\varepsilon}^{0,H}$. We recall that $\boldsymbol{\sigma}^{*} \cdot \mathbf{n}^{K}|_{\partial K}$ is piecewise constant for all $K$ in $\mathcal{T}_{H}$. Now, from the hybrid formulation of (22), the fact that $T_{K}^{\varepsilon}\lambda_{\varepsilon}$ and $T_{K}^{\varepsilon}\lambda_{\varepsilon}^{h}$ belong to $L_{0}^{2}(K)$ for all $K \in \mathcal{T}_{H}$, and Cauchy-Schwarz inequality

$$
\begin{aligned}
\|u_{\varepsilon}^{0} - u_{\varepsilon}^{0,H}\|_{0, \mathcal{T}_H}^2 &= (\nabla \cdot \boldsymbol{\sigma}^*, u_{\varepsilon}^0 - u_{\varepsilon}^{0,H})_{\mathcal{T}_H} \\
&= \sum_{K \in \mathcal{T}_H} (\boldsymbol{\sigma}_{\partial K}^* \cdot \mathbf{n}^K, u_{\varepsilon}^0 - u_{\varepsilon}^{0,H})_{\partial K} \\
&= - \sum_{K \in \mathcal{T}_H} (\boldsymbol{\sigma}_{\partial K}^* \cdot \mathbf{n}^K, T_K^\varepsilon \lambda_\varepsilon - T_K^\varepsilon \lambda_\varepsilon^h)_{\partial K} \\
&= -(\boldsymbol{\sigma}^*, \nabla_H (T_\varepsilon \lambda_\varepsilon - T_\varepsilon \lambda_\varepsilon^h))_{\mathcal{T}_H} + (\nabla \cdot \boldsymbol{\sigma}^*, T_\varepsilon \lambda_\varepsilon - T_\varepsilon \lambda_\varepsilon^h)_{\mathcal{T}_H} \\
&= -(\boldsymbol{\sigma}^*, \nabla_H (T_\varepsilon \lambda_\varepsilon - T_\varepsilon \lambda_\varepsilon^h))_{\mathcal{T}_H} \\
&\le \|\boldsymbol{\sigma}^*\|_{0, \mathcal{T}_H} \|\nabla_H (T_\varepsilon \lambda_\varepsilon - T_\varepsilon \lambda_\varepsilon^h)\|_{0, \mathcal{T}_H} \\
&\le c \|u_\varepsilon^0 - u_\varepsilon^{0,H}\|_{0, \mathcal{T}_H} \|\nabla_H (T_\varepsilon \lambda_\varepsilon - T_\varepsilon \lambda_\varepsilon^h)\|_{0, \mathcal{T}_H}
\end{aligned}
$$

where we used the regularity of the mesh. Collecting the previous results, we obtain from (34) the existence of $c$ such that

$$ \|u_{\varepsilon} - u_{\varepsilon}^{h}\|_{0, \mathcal{T}_{H}} \leq c |u_{\varepsilon} - u_{\varepsilon}^{h}|_{1, \mathcal{T}_{H}}, $$

and the result follows.

□
**Remark 4.** Optimality in the $L^2(\Omega)$ norm can be achieved using the classical duality argument (see [2] for details). As the additional regularity of the dual solution is mandatory, this is not an option when it comes to uniform error estimates with respect to the small parameter $\epsilon$. Another option is to assume the existence of $\boldsymbol{\sigma}^* \in H(\text{div}; \Omega)$ with $\nabla \cdot \boldsymbol{\sigma}^* = u_{\epsilon}^0 - u_{\epsilon}^{0,H}$ and $\boldsymbol{\sigma}^* \cdot \boldsymbol{n}_K|_{\partial K} \in \mathbb{P}_0(\partial K)$ such that $\|\boldsymbol{\sigma}^*\|_{0,K} \le cH\|\nabla \cdot \boldsymbol{\sigma}^*\|_{0,K}$ for all $K$ in $\mathcal{T}_H$. As a result, and following the same steps, we get $\|u_{\epsilon} - u_{\epsilon}^h\|_{0,\mathcal{T}_H} \le cH|u_{\epsilon} - u_{\epsilon}^h|_{1,H}$ without assuming extra regularity (c.f. [23]). Numerics in Section 5 verify the optimality of the error in the $L^2(\Omega)$ norm in terms of $H$, but the existence of $\boldsymbol{\sigma}^*$ remains an open question and deserves further investigation.

**4.2. Interpolation results.** Assume that, for each $K \in \mathcal{T}_H$, there exists $\{\mathcal{T}_h^K\}_{h>0}$ a shape-regular family of conforming simplicial partitions of $K$ associated to $\{\mathcal{E}_h^K\}_{h>0}$ in the sense that, for each $F \in \mathcal{E}_h^K$ there exists an element $\tau_F \in \mathcal{T}_h^K$ with $\tau_F \cap \partial K = F$. We denote the global partition

$$ (35) \qquad \mathcal{T}_h := \bigcup_{K \in \mathcal{T}_H} \mathcal{T}_h^K. $$

Also, we denote $\mathcal{F}_h^K$ the set of faces on $\mathcal{T}_h^K$, and $\mathcal{F}_0^K$ the set of internal faces. To each face $\gamma \in \mathcal{F}_h^K$ we associate a normal vector $\boldsymbol{n}_\gamma^\tau$, taking care to ensure this is facing outward on $\partial\tau$.

We need some standard interpolations and projections:

* Given $E \in \mathcal{E}$, let $\pi_E^1: L^2(E) \to \mathbb{P}_1(E) \cap C^0(E)$ be the local orthogonal projection $L^2$ into $\mathbb{P}_1(E)$ the space of of piecewise linear functions on each $F \in \mathcal{E}_h \cap E$. Also, we define the global $L^2$ projection $\pi^1(\cdot)$ such that $\pi^1(\mu)|_E := \pi_E^1(\mu)$, and then it holds

$$ (36) \qquad \|\pi^1(\mu)\|_{0,\partial\mathcal{T}_H} \le \|\mu\|_{0,\partial\mathcal{T}_H} \quad \text{for all } \mu \in L^2(\partial\mathcal{T}_H). $$

* Given $K \in \mathcal{T}_H$, let $\Pi_K^1: L^2(K) \to \mathbb{P}_1(\mathcal{T}_h^K) \cap C_0(K)$ be the $L^2$ orthogonal projection on $\mathbb{P}_1(\mathcal{T}_h^K)$ the space of piecewise linear functions with respect to the partition $\mathcal{T}_h^K$. Also, we define the global $L^2$ projection $\Pi^1(\cdot)$ such that $\Pi^1(v)|_K := \Pi_K^1(v)$, and then for all $v \in H^s(\mathcal{T}_H)$ and $0 \le m \le s$ and $s \in \{0, 1\}$, it holds

$$ (37) \qquad \begin{aligned} \|v - \Pi^1(v)\|_{m, \mathcal{T}_H} &\le c h^{s-m} |v|_{s, \mathcal{T}_H} \quad \text{for } m = 0, 1, \\ \|v - \Pi^1(v)\|_{0, \partial \mathcal{T}_H} &\le c h^{1/2} |v|_{1, \mathcal{T}_H}, \end{aligned} $$

with $h = \max_{K \in \mathcal{T}_H} \max_{\tau \in \mathcal{T}_h^K} h_\tau$, where $h_\tau$ is the diameter of element $\tau \in \mathcal{T}_h^K$.
• Let $\mathcal{P}_h: H^s(\mathcal{T}_h) \to \mathbb{P}_k(\mathcal{T}_h) \cap C_0(\bar{\Omega})$ be the Scott-Zhang projection $\mathcal{P}_h$ where $\mathbb{P}_k(\mathcal{T}_h)$ is the space of piecewise polynomial functions of degree up to $k \ge 1$ on $\mathcal{T}_h$ (c.f. [25], [11][Lemma 1.130]). Notably, given $K \in \mathcal{T}_H$ we have that $\mathcal{P}_h v|_K \in C_0(K)$ and $\mathcal{P}_h v|_{\tau} \in \mathbb{P}_k(\tau)$ for every $\tau \in \mathcal{T}_h^K$, and the following approximation property holds

$$ (38) \qquad \|v_{\tau} - \mathcal{P}_h(v)|_{\tau}\|_{m,\tau} \le C h_{\tau}^{s-m} |v_{\tau}|_{s,\Delta_{\tau}^K} \quad \text{for all } v \in H^s(K), $$

with $0 \le m \le s \le k + 1$, where $\Delta_{\tau}^K := \{\hat{\tau} \in \mathcal{T}_h^K : \hat{\tau} \cap \tau \neq 0\}$, and $h_{\tau}$ is the diameter of element $\tau \in \mathcal{T}_h^K$.

The next result is fundamental to establish rates of convergence of the MHM method using space $\Lambda_h^\ell$. It is an extension of [24] and [5] to deal with the case of continuous polynomial interpolation spaces on faces of elements $K \in \mathcal{T}_H$.

**Lemma 5 (Interpolation).** Suppose $w \in H^{s+2}(\mathcal{T}_H) \cap H_0^1(\Omega)$, $\mathcal{A}_\varepsilon \nabla_H w \in H^{s+1}(\mathcal{T}_H)$, with $0 \le s \le \ell$ and $\ell \ge 1$ if $\Lambda_h^\ell$ is the continuous space and $\ell \ge 0$ otherwise, and $\mathcal{A}_\varepsilon \nabla_H w \in H(\text{div};\Omega)$. Let $\mu \in \Lambda$ be defined by $\mu_E := (\mathcal{A}_\varepsilon \nabla_H w_K \cdot \mathbf{n}^K)|_E$ for each $E \in \mathcal{E}$. Then, there exist $\mu^h \in \Lambda_h^\ell$ and $c$ such that

$$ (39) \qquad \|\mu - \mu^h\|_{\Lambda} \le c h^{s+1} |\mathcal{A}_\varepsilon \nabla_H w|_{s+1, \mathcal{T}_H} \quad \text{and} \quad \mu - \mu^h \in \mathcal{N}, $$

where $\mathcal{N}$ is

$$ (40) \qquad \mathcal{N} := \{\mu \in \Lambda : (\mu, v_0)_{\partial \mathcal{T}_H} = 0 \text{ for all } v_0 \in V_0\}. $$

*Proof.* The estimate (39) follows from [5][Lemma 3] in the case $\Lambda_h^\ell$ is discontinuous. Next, we prove the result when $\Lambda_h^\ell$ is continuous, which builds on the strategy originally proposed in [24] for one-element mesh faces and extends to the submesh face case in [5]. Let $w \in H^{s+2}(\mathcal{T}_H)$ and $E \in \mathcal{E}$. We define $\chi_E := \mathcal{A}_\varepsilon \nabla_H w \cdot \mathbf{n}_E \in H^{s+1}(\mathcal{T}_H)$ where $\mathbf{n}_E$ is the trivial extension of the normal vector $\mathbf{n}^K|_E$ from $E$ to a constant function in the whole $K$, with $|\mathbf{n}_E| = 1$. Let $\mu \in L^2(\mathcal{E})$ be such that $\mu_E = \chi_E|_E$ on each $E \in \mathcal{E}$.

Define $\mu^h \in \Lambda_h^\ell$ such that $\mu_F^h := \mathcal{P}_h(\chi_E)|_F$, for all $F \in \mathcal{E}_h$, where $\mathcal{P}_h(\cdot)$ is the Scott-Zhang interpolation operator on the space of piecewise continuous polynomials of degree $\ell$ with respect to the partition $\mathcal{T}_h$. We recall that the Scott-Zhang operator is a projection (c.f. [11][Lemma 1.130]), i.e $\mathcal{P}_h(v) = v$ for all $v \in \mathbb{P}_\ell(\mathcal{T}_h) \cap C_0(\bar{\Omega})$, and then we can follow closely [5] to arrive at

$$ (41) \qquad \|\mu - \mu^h\|_{\Lambda} \le C h^{s+1} |\mathcal{A}_\varepsilon \nabla_H w|_{s+1, \mathcal{T}_H} \quad \text{for all } 0 \le s \le \ell. $$
Given any $K \in \mathcal{T}_H$ and $F \subset \partial K$, let $\tau_F$ be the unique element of $\mathcal{T}_h^K$ such that $\partial\tau_F \cap \partial K = F$. Using a trace inequality, the definition of $\mu$ and $\mu^h$, and from (38), we obtain

$$ \begin{align} \|\mu - \mu^h\|_{0,\partial\mathcal{T}_H}^2 &= \sum_{K \in \mathcal{T}_H} \sum_{F \in \mathcal{E}_h^K} \|\chi_E|_F - \mathcal{P}_h(\chi_E)|_F\|_{0,F}^2 \nonumber \\ &\le c \sum_{K \in \mathcal{T}_H} \sum_{F \in \mathcal{E}_h^K} \left( \frac{1}{h_{\tau_F}} \|\chi_E|_{\tau_F} - \mathcal{P}_h(\chi_E)|_{\tau_F}\|_{0,\tau_F}^2 \right. \nonumber \\ &\qquad \left. + h_{\tau_F} \|\chi_E|_{\tau_F} - \mathcal{P}_h(\chi_E)|_{\tau_F}\|_{1,\tau_F}^2 \right) \nonumber \\ &\le c h^{2(s+1)-1} |\chi_E|_{s+1,\mathcal{T}_H}^2. \tag{42} \end{align} $$

Hence, it holds

$$ (43) \qquad \|\mu - \mu^h\|_{0,\partial\mathcal{T}_H} \le C h^{s+1/2} |\mathcal{A}_\varepsilon \nabla_H w|_{s+1,\mathcal{T}_H}, $$

where we used $|\chi_E|_{s+1,K} = |\mathcal{A}_\varepsilon \nabla_H w \cdot n_E|_{s+1,K} \le |\mathcal{A}_\varepsilon \nabla_H w|_{s+1,K}$, for all $K \in \mathcal{T}_H$.

Unfortunately $\mu - \mu^h \notin \mathcal{N}$. Then, we propose a Fortin-type interpolation operator for $\mu$ relying on $\mu^h$ which keeps its approximation properties. Notably, let $\tilde{\mu}^h \in \Lambda_h$ be defined by

$$ \tilde{\mu}_E^h := \mu_E^h + \pi_E^1(\mu_E - \mu_E^h) \quad \text{for all } E \in \mathcal{E}. $$

Observe that

$$ \int_{\partial K} (\mu_{\partial K} - \tilde{\mu}_{\partial K}^h) d\boldsymbol{x} = \sum_{E \subset \partial K} \int_E (\mu_E - \mu_E^h - \pi_E^1(\mu_E - \mu_E^h)) d\boldsymbol{x} = 0 $$

and then $\mu - \tilde{\mu}^h \in \mathcal{N}$. It remains to prove that the function $\tilde{\mu}^h$ inherits the approximation properties of $\mu^h$. Next, using (19), (37), Cauchy-Schwartz inequality, (36), and (41), and
(43), we get

$$
\begin{align*}
\|\mu - \tilde{\mu}^h\|_{\Lambda} &= \|\mu - \mu^h - \pi_E^1(\mu - \mu^h)\|_{\Lambda} \\
&= \sup_{v \in V} \frac{(\mu - \mu^h - \pi_E^1(\mu - \mu^h), v)_{\partial \mathcal{T}_H}}{\|v\|_V} \\
&= \sup_{v \in V} \frac{(\mu - \mu^h - \pi_E^1(\mu - \mu^h), v - \Pi^1(v))_{\partial \mathcal{T}_H}}{\|v\|_V} \\
&\le \sup_{v \in V} \frac{(\mu - \mu^h, v - \Pi^1(v))_{\partial \mathcal{T}_H}}{\|v\|_V} + \sup_{v \in V} \frac{(\pi_E^1(\mu - \mu^h), v - \Pi^1(v))_{\partial \mathcal{T}_H}}{\|v\|_V} \\
&\le \sup_{v \in V} \frac{\|\mu - \mu^h\|_{\Lambda} \|v - \Pi_1(v)\|_V}{\|v\|_V} + \sup_{v \in V} \frac{\|\pi_E^1(\mu - \mu^h)\|_{0, \partial \mathcal{T}_H} \|v - \Pi^1(v)\|_{0, \partial \mathcal{T}_H}}{\|v\|_V} \\
&\le C \left( \|\mu - \mu^h\|_{\Lambda} + h^{1/2} \|\pi_E^1(\mu - \mu^h)\|_{0, \partial \mathcal{T}_H} \right) \\
&\le C \left( \|\mu - \mu^h\|_{\Lambda} + h^{1/2} \|\mu - \mu^h\|_{0, \partial \mathcal{T}_H} \right) \\
&\le C h^{s+1} |\mathcal{A}_{\varepsilon} \nabla_H w|_{s+1, \mathcal{T}_H},
\end{align*}
$$

and the result follows.

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

4.3. **Error estimates.** Now, we address the main convergence result. For that, we first set up the asymptotic regime in which the future convergence results will be proved. Let $d_K$ be the side of the largest square (or cube in 3D) contained in $K \in \mathcal{T}_h$. Hereafter, we shall assume that

$$
(44) \qquad \inf_{K \in \mathcal{T}_h} d_K > c \varepsilon,
$$

where $c > 1$ is of order one. Also, we recall that from the regularity theory for elliptic equations and the convexity of $K \in \mathcal{T}_H$, the following conditions hold

$$
(45) \qquad \hat{T}_K^0 f_K \in H^2(K) \quad \text{and} \quad \|\hat{T}_K^0 f_K\|_{2,K} \le c \|f_K\|_{0,K}.
$$

Next, we consider the analysis for the homogeneous Neumann boundary case ($b=0$) and address the general case in Remark 12. For that, we first revisit an important Lemma from [23]. Its proof follows the same strategy proposed in [23, Lemma 4] and Theorem 1.

**Lemma 6.** Assume (44) holds, the coefficients $a_{ij}$ are $Y$-periodic $u_0 \in H^2(\Omega)$ and $\chi^j$
$\in W_{per}^{1,q}(Y)$ satisfies (11) Also, assume $\hat{T}_0 f|_K \in H^2(K)$ for all $K \in \mathcal{T}_h$. Then,

$$
\left\| T_{\varepsilon} \lambda_{\varepsilon} - T_0 \lambda_0 - \varepsilon \chi_{\varepsilon}^j \frac{\partial T_0 \lambda_0}{\partial x_j} \right\|_{1, \mathcal{T}_H} \le c(p') \left(\frac{\varepsilon}{H}\right)^{1/2-1/p'} \|\chi\|_{1,q,Y} (\|u_0\|_2 + \|\hat{T}_0 f\|_{2, \mathcal{T}_H})
$$

where $2 < p' < K_d$, and $K_d$ satisfies (13). Furthermore, the constant $c(p')$ is such that
$c(p') \to \infty$ when $p' \to K_d$.
The next lemma improves [23, Lemma 7], allowing us to study the convergence of the MHM method with piecewise discontinuous polynomials $\Lambda_h^\ell$. Unlike [23, Lemma 7] the strategy adopted here avoids estimating the broken semi-norm $H^1$ of $T_\varepsilon\lambda_0 - T_\varepsilon\mu^h$. Consequently, we do not work with the asymptotic expansion approximation for $T_\varepsilon\mu^h$, which requires less regularity than $T_0\mu^h$.

**Lemma 7.** Let $u_\varepsilon$ be the exact solution of (3) with $b=0$ and $u_\varepsilon^h$ the approximate solution given by (26). Assume the coefficients $a_{ij}$ are Y-periodic, $u_0 \in H^2(\Omega)$ and $\chi^j \in W_{per}^{1,q}(Y)$, with $q > d$, $f \in L^2(\Omega)$, and that (44) and (45) hold. Then,

$$ (46) \qquad \begin{aligned} |u_\varepsilon - u_\varepsilon^h|_{1, \tau_H} &\le c(p') \left(\frac{\varepsilon}{H}\right)^{1/2-1/p'} \|\mathbf{x}\|_{1,q,Y} (\|u_0\|_2 + \|f\|_0) \\ &\quad + \inf_{\mu^h \in \tilde{\Lambda}_h^\ell} \|\lambda_0 - \mu^h\|_{\Lambda}, \end{aligned} $$

where

$$ \tilde{\Lambda}_h^\ell := \left\{ \mu^h \in \Lambda_h^\ell : \int_{\partial K} \mu^h ds = \int_K f x d\mathbf{x} \quad \forall K \in \mathcal{T}_H \right\}, $$

and we can choose $p' > 2$ satisfying (13), and the constant $c(p')$ depends on $p'$ ($c(p') \to \infty$ when $p' \to K_d$).

*Proof.* We recall from [2, Lemma 3.5] that

$$ (47) \qquad |u_\varepsilon - u_\varepsilon^h|_{1, \tau_H} = |T_\varepsilon \lambda_\varepsilon - T_\varepsilon \lambda_\varepsilon^h|_{1, \tau_H} \le c \inf_{\mu^h \in \tilde{\Lambda}_h^\ell} \|\lambda_\varepsilon - \mu^h\|_{\Lambda}. $$

Next, we restrict the infimum in (47) to $\tilde{\Lambda}_h^\ell$ and use that $\lambda_\varepsilon - \mu^h$ and $\lambda_\varepsilon - \lambda_0$ belong to the space defined by (40) to get

$$ (48) \qquad \begin{aligned} |u_\varepsilon - u_\varepsilon^h|_{1, \tau_H} &\le c \inf_{\mu^h \in \tilde{\Lambda}_h^\ell} \|\lambda_\varepsilon - \mu^h\|_{\Lambda} \\ &\le c \left( \|\lambda_\varepsilon - \lambda_0\|_{\Lambda} + \inf_{\mu^h \in \tilde{\Lambda}_h^\ell} \|\lambda_0 - \mu^h\|_{\Lambda} \right) \\ &\le c \left( |T_\varepsilon \lambda_\varepsilon - T_\varepsilon \lambda_0|_{1, \tau_H} + \inf_{\mu^h \in \tilde{\Lambda}_h^\ell} \|\lambda_0 - \mu^h\|_{\Lambda} \right), \end{aligned} $$

where we used that if $\mu$ belongs to the space defined by (40) then

$$ (49) \qquad \|\mu\|_{\Lambda} \le c |T_{\varepsilon} \mu|_{1, \tau_H}. $$

Indeed, from the definitions of $\|\cdot\|_{\Lambda}$ and the operator $T_\varepsilon$, the ellipticity condition (4), and the fact that $\nabla \cdot (A_\varepsilon \nabla T_\varepsilon \mu) = 0$ in each $K \in \mathcal{T}_h$, we obtain (49). The term $|T_\varepsilon \lambda_\varepsilon - T_\varepsilon \lambda_0|_{1, \tau_H}$
in (48) is bounded as follows

$$ (50) \qquad |T_{\varepsilon}\lambda_{\varepsilon} - T_{\varepsilon}\lambda_{0}|_{1,\mathcal{T}_H} \le \left|T_{\varepsilon}\lambda_{\varepsilon} - T_{0}\lambda_{0} - \varepsilon \chi_{\varepsilon}^{j} \frac{\partial T_{0}\lambda_{0}}{\partial x_{j}}\right|_{1,\mathcal{T}_H} + \left|T_{0}\lambda_{0} + \varepsilon \chi_{\varepsilon}^{j} \frac{\partial T_{0}\lambda_{0}}{\partial x_{j}} - T_{\varepsilon}\lambda_{0}\right|_{1,\mathcal{T}_H}. $$

Next, Lemma 6 provides an estimate for the first term in the right-hand side of (50), while the second one is estimated using Theorem 1 as $T_0\lambda_0|_K \in H^2(K)$. The final result follows summing up all contributions. □

We are ready to present the main result of this section.

**Theorem 8 (Convergence $H \gg \varepsilon$).** Let $u_\varepsilon$ be the solution of (3) (or (22)) with $b=0$, $\sigma_\varepsilon := -A_\varepsilon \nabla u_\varepsilon$, and $u_\varepsilon^h$ given in (26) and $\sigma_\varepsilon^h$ in (27). Assume the coefficients $a_{ij}$ are Y-periodic, $u_0 \in H^2(\Omega)$, $\chi^j \in W_{per}^{1,\infty}(Y)$, $f \in L^2(\Omega)$, and (44) holds. Then,

$$ (51) \qquad \|u_\varepsilon - u_\varepsilon^h\|_V \le \left[ c_1(p') \left(\frac{\varepsilon}{H}\right)^{1/2-1/p'} + ch \right] \|\chi\|_{1,\infty,Y} (\|u_0\|_2 + \|f\|_0), $$

$$ (52) \qquad \|\sigma_{\varepsilon} - \sigma_{\varepsilon}^{h}\|_{\text{div}} \le \left[ c_{1}(p') \left(\frac{\varepsilon}{H}\right)^{1/2-1/p'} + ch \right] \|\chi\|_{1,\infty,Y} (\|u_{0}\|_{2} + \|f\|_{0}), $$

where $p'$ satisfies (13).

*Proof.* The proof follows the strategy proposed in [23, Theorem 8], using Lemma 3, Lemma 5, Lemma 7 and Theorem 1. □

We present a result of convergence assuming less regularity of the functions $\chi^j$. However, unlike the Theorem 1, the Sobolev embedding here uses fractional Sobolev spaces. To get the result, we make an additional assumption about the triangulation family $\mathcal{T}_H$ and use a scaling argument to get the correct dependence of the Sobolev embedding constant on the size of $K$.

**Theorem 9 (Convergence $H \gg \varepsilon$ under low regularity).** Let $u_\varepsilon$ be the solution of (22) with $b=0$ and let $u_\varepsilon^h$ be its numerical approximation defined by (26). Assume the coefficients $a_{ij}$ are Y-periodic, $u_0 \in H^2(\Omega)$, $\chi^j \in W_{per}^{1,q}(\Omega)$, $q > d$, $f \in L^2(\Omega)$, and (44) holds. Also, assume that every element $K \in \mathcal{T}_H$ is of the form $H\tilde{K}$, where $\tilde{K}$ belongs to a finite collection of elements of size one. Then,

$$ \|u_\varepsilon - u_\varepsilon^h\|_V \le c(p') \left( \left(\frac{\varepsilon}{H}\right)^{1/2-1/p'} + \frac{h^{1-d/q}}{H^{d/q}} \right) \|\chi\|_{1,q,Y} (\|u_0\|_2 + \|f\|_0), $$

$$ \|\sigma_{\varepsilon} - \sigma_{\varepsilon}^{h}\|_{\text{div}} \le c(p') \left( \left(\frac{\varepsilon}{H}\right)^{1/2-1/p'} + \frac{h^{1-d/q}}{H^{d/q}} \right) \|\chi\|_{1,q,Y} (\|u_{0}\|_{2} + \|f\|_{0}), $$

where $p'$ satisfies (13).
Proof. The proof follows the same strategy from [23][Theorem 9]. Next, we present the main estimates needed to obtain our result considering the precise dependence of the constants on $\mathcal{T}_H$. We start recalling the definition of a fractional Sobolev norm

$$ \|f\|_{s,2,\tilde{K}}^2 := \|f\|_{0,2,\tilde{K}}^2 + \int_{\Omega} \int_{\Omega} \frac{|f(x) - f(y)|^2}{|x - y|^{d+2s}} dx dy. $$

It is well known that the following Sobolev embedding

$$ \|f\|_{0,p,\tilde{K}} \le c \|f\|_{s,2,\tilde{K}} \quad \text{with} \quad \frac{2d}{d-2s} = p \text{ and hence } s = \frac{d}{q} $$

holds. Using this result on a mesh element $K$, applying a scaling argument and the hypothesis that every element $K \in \mathcal{T}_H$ is of the form $H\tilde{K}$, where $\tilde{K}$ belongs to a finite collection of elements of size one, we obtain that

$$ \|f\|_{0,2,K} \le c H^{-d/2} \|f\|_{s,2,K}. $$

Now, let $\frac{1}{p} + \frac{1}{q} = \frac{1}{2}$, then it holds

$$ \begin{aligned} \left\| \varepsilon \frac{\partial \chi_{\varepsilon}^j}{\partial x_k} \frac{\partial (T_0 \lambda_0 - T_0 \mu^h)}{\partial x_j} \right\|_0 &\le \|\chi\|_{1,q,Y} \|T_0 \lambda_0 - T_0 \mu^h\|_{1,p,\mathcal{T}_H} \\ &\le c H^{d/p} H^{-d/2} \|\chi\|_{1,q,Y} \|T_0 \lambda_0 - T_0 \mu^h\|_{1+s,\mathcal{T}_H} \\ &\le c \frac{h^{1-d/q}}{H^{d/q}} \|\chi\|_{1,q,Y} \|u_0\|_2. \end{aligned} $$

The other terms are bounded following the proof of Theorem 1 with similar changes. $\square$

**Remark 10** (Convergence under low regularity). *From the Theorem 9, we conclude that, under the assumption of mild regularity at the local level, the MHM method converges in the space-based approach ($h \to 0$ and $H$ fixed). $\square$*

**Remark 11** (Convergence $H \le \varepsilon$). *The (optimal) error estimates in the asymptotic regime $H \le \varepsilon$ turn out to be a direct consequence of (32) and the Lemma 5. In fact, assuming that $u_\varepsilon$ satisfies the regularity conditions of the Lemma 5, then*

$$ \|u_{\varepsilon} - u_{\varepsilon}^{h}\|_{V} \le c h^{s+1} |u_{\varepsilon}|_{s+2, \mathcal{T}_{H}}, $$

with $0 \le s \le l$ and $l \ge 1$ if $\Lambda_h^l$ is the continuous space and $l \ge 0$ otherwise. If $u_\varepsilon$ is more regular, it follows from (32) and [7][Theorem 3.3] that the error above super-converges to zero as $O(h^{s+3/2})$. $\square$

**Remark 12** (Non-homogeneous Neumann boundary condition). *Our analysis extends to treat the case $b \ne 0$ for which the additional term $\|\overline{T}_0 b\|_{2,H}$ should be included on the right*
side of the estimates in Theorems 8 and 9. Typically, one would expect from the theory of elliptic equations that

$$ \|\bar{T}_0 b\|_{2,K} \le c \|b\|_{H^{3/2}(K)}. $$

However, in the case of a non-homogeneous Neumann boundary condition, it is not clear the dependence of the constant c on the domain K, and then we suggest keeping $\|\bar{T}_0 b\|_{2,H}$ in the error estimates.

## 5. NUMERICAL RESULTS

5.1. **Assessing convergence.** The domain is a unit square with prescribed homogeneous Dirichlet boundary conditions, $f(x) = \sin(x_1) \sin(x_2)$, and coefficient given by

$$ \mathcal{A}_{\varepsilon}(x) = \left[ 1 + 100 \cos^2 \left( \frac{\pi x_1}{\varepsilon} \right) \sin^2 \left( \frac{\pi x_2}{\varepsilon} \right) \right] \mathcal{I}, $$

where $\varepsilon = \frac{\pi}{150}$ is the small parameter defining the periodicity, and $\mathcal{I}$ is the identity matrix. The reference solution is constructed using a mesh of 16,777,216 quadrilateral bilinear elements. The MHM method is validated using quadrilateral elements with continuous linear interpolation on faces to approximate the Lagrange multipliers ($\ell = 1$). The multiscale basis functions and $\hat{T}_{\varepsilon}f$ (see Figure 1 for an illustrative example) are approximated at the local level by the standard Galerkin method over the bilinear continuous polynomial space defined on structured sub-meshes. The sub-meshes are selected such that the multiscale base and $\hat{T}_{\varepsilon}f$ are accurately approximated so that the underlying errors do not impact the MHM method.

<figure><img src="image_4.png" /><figcaption>FIGURE 1. A typical multiscale basis function (left) and $\hat{T}_{\varepsilon}f$ (right) on a quadrangular element of $\mathcal{T}_H$. Here $\varepsilon = \frac{1}{64}$ and $\ell = 1$.</figcaption></figure>
Specifically, we define a structured quadrilateral mesh with $H = \frac{1}{8}$ and progressively increase the degrees of freedom on each edge. We calculate the relative error in the broken $H^1$-norm and in the $L^2$ norm using this strategy (space-based) and compare it with that obtained from successive mesh refinements (mesh-based). The results are shown in Figure 2, wherein we analyze the result from the perspective of the diameter $h$ corresponding to the edge partition's diameter. Remember that in the space-based case, the mesh is fixed (with diameter $\frac{1}{8}$), and we define $h = H$ in the mesh-based strategy. In Figure 3, we perform the same analysis, but now for the number of degrees of freedom $N_{DOF}$.

We observe that in the space-based enhancement approach, the underlying error is drastically reduced with the result that the error divergence region (of order $h^{-1}$) is no longer displayed. Also, considerably fewer degrees of freedom are needed to reach a given error threshold. It is worth mentioning that such behavior is predicted by current theory, and it is achieved without any oversampling technique.

<figure><img src="image_5.png" /><figcaption>FIGURE 2. Convergence history with respect to $h$ in the relative $L^2(\Omega)$ (left) and $H^1(\mathcal{T}_H)$ (right) norms. Comparison between the one piecewise linear interpolation per edge case with mesh refinement (mesh-based) and the multiple piecewise linear interpolation per edge case on a fixed mesh (space-based).</figcaption></figure>

Figure 4 illustrates a comparison between the isolines of the reference and MHM solutions with a one-element mesh on faces and $\ell = 1$.
<figure><img src="image_6.png" /><figcaption>FIGURE 3. Convergence history with respect to $N_{DOF}$ in the relative $L^2(\Omega)$ (left) and $H^1(\mathcal{T}_H)$ (right) norms. Comparison between the one piecewise linear interpolation per edge case with mesh refinement (mesh-based) and the multiple piecewise linear interpolation per edge case on a fixed mesh (space-based). The latter induces a tremendous improvement in the quality of the numerical results with fewer degrees of freedom.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 4. Isolines of the MHM solution overlaid with the rectangular mesh with 64 squares (left) alongside the reference solution (right). Here $\epsilon = \frac{1}{64}$, $\ell = 1$ and $h = H$.</figcaption></figure>

5.2. **A porous media fluid flow problem.** We use the permeability defined in layer 36$^{th}$ from Model 2 of the 10th Comparative Solution Project of the Society of Petroleum Engineers [22] (see Figure 5) to validate the MHM method with continuous linear interpolation
(ℓ = 1) on the faces of a quadrilateral mesh with 66 squares. The model problem is defined by an entry pressure imposed equal to 1 at the bottom side and 0 at the top. We set the homogeneous Neumann conditions on the other two boundaries, and f = 0. We consider a reference solution (see Figure 5) obtained with continuous piecewise cubic polynomials on a quadrilateral mesh with 1,081,344 elements representing 9,738,625 degrees of freedom.

<figure><img src="image_8.png" /><figcaption>FIGURE 5. Logarithm of the permeability proposed in [22] (left) and the isolvalues of the reference solution (right).</figcaption></figure>

In Figure 6 we compare the isolines of the reference solution and the MHM solution on a coarse mesh with 66 squares, accounting for 4,257 degrees of freedom. We observe qualitatively equivalent results from the MHM method even using much fewer degrees of freedom (Figure 6). To see this more clearly, we depict in Figure 7 the profiles of the solutions obtained from three refinements of *h* and the reference solution at *x*₁ = 199, and zoom those profiles in Figure 8.
<figure><img src="image_9.png" /><figcaption>FIGURE 6. Isolines of pressure from the MHM method with the overlaid quadrilateral mesh composed by 66 squares (left) and the reference solution (right). Here $\ell = 1$ and $h = \frac{1}{32} H$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 7. Comparison of the profiles at the $x_1 = 199$ between three MHM solutions obtained different levels of refinement on edges and the reference solution. Here $\ell = 1$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 8. Zoomed solutions in the area around $x_2 = 1,650$ in Figure 7.</figcaption></figure>

## 6. CONCLUSIONS

In previous works (c.f. [5, 7]), the MHM method proved to be convergent with respect to $h$ when the space-based strategy is used (i.e. $h \to 0$ with fixed $H$) assuming the asymptotic regime $H \le \varepsilon$. The present work proved that convergence also occurs in the asymptotic regime $H \gg \varepsilon$, without resonance, when we adopt piecewise continuous or discontinuous polynomial interpolation on faces. However, numerical pollution is still present and represented by the factor $(\varepsilon/H)^{1/2}$. Possibilities to reduce such error exist by using nonlocal strategies to select the multiscale basis or by involving physics in the construction of the basis functions for the Lagrange multipliers. Despite the pollution error, the numerical results point to a drastic decrease in the number of degrees of freedom to reach a given error threshold compared to the standard convergence way of refining the one-level mesh ($H \to 0$) which do have the presence of resonance error. Discontinuous interpolation on faces seems to avoid resonance errors in the same way and produce results qualitatively similar to the case of continuous interpolation (c.f. [23] for numerics).

## 7. APPENDIX

*Proof of Theorem 1.* This section fills in the gap in [23][Theorem 1], where only the Dirichlet boundary condition case was detailed. Here, we consider (3) with the prescribed Neumann boundary condition and follow [19, Section 1.4]. This case was not dealt with in the proof presented in [23].
First, we establish the existence of a function $\mathbf{z}_{\varepsilon}^1 \in H^1(\Omega)$ such that

$$ (53) \qquad \int_{\Omega} \mathbf{z}_{\varepsilon}^1 \cdot \nabla \phi \, dx = \int_{\Omega} \mathcal{A}_{\varepsilon} \nabla u_{\varepsilon} \cdot \nabla \phi \, d\mathbf{x} \quad \text{for all } \phi \in H^1(\Omega), $$

which is “close” to $\mathcal{A}_{\varepsilon} \nabla u_{\varepsilon}^1$ in the sense that it satisfies

$$ (54) \qquad \|\mathbf{z}_{\varepsilon}^1 - \mathcal{A}_{\varepsilon} \nabla u_{\varepsilon}^1\|_0 \le \left[ c\varepsilon + c(p') \left( \frac{L_{\Omega}\varepsilon}{|\Omega|} \right)^{1/2-1/p'} \right] \|\mathbf{\chi}\|_{1,q,Y} \|u_0\|_2. $$

Observe that, owing to (53) and (54) we get

$$ \begin{aligned} \int_{\Omega} \mathcal{A}_{\varepsilon} \nabla (u_{\varepsilon} - u_{\varepsilon}^{1}) \cdot \nabla (u_{\varepsilon} - u_{\varepsilon}^{1}) \, d\mathbf{x} &= \int_{\Omega} (\mathbf{z}_{\varepsilon}^{1} - \mathcal{A}_{\varepsilon} \nabla u_{\varepsilon}^{1}) \cdot \nabla (u_{\varepsilon} - u_{\varepsilon}^{1}) \, d\mathbf{x} \\ &\le \left[ c\varepsilon + c(p') \left( \frac{L_{\Omega}\varepsilon}{|\Omega|} \right)^{1/2-1/p'} \right] \|\mathbf{\chi}\|_{1,q,Y} \|u_0\|_2 |u_{\varepsilon} - u_{\varepsilon}^{1}|_1, \end{aligned} $$

and, then, the ellipticity of $\mathcal{A}_{\varepsilon}$ yields

$$ (55) \qquad |u_{\varepsilon} - u_{\varepsilon}^{1}|_{1} \le \left[ c\varepsilon + c(p') \left( \frac{L_{\Omega}\varepsilon}{|\Omega|} \right)^{1/2-1/p'} \right] \|\mathbf{\chi}\|_{1,q,Y} \|u_0\|_2. $$

The proof of $\mathbf{z}_{\varepsilon}^1$ fulfilling (53) and (54) is constructive. Let $\mathbf{y} := \mathbf{x}/\varepsilon$ be the slow variable, and observe that

$$ \begin{aligned} (\mathcal{A}_{\varepsilon} \nabla u_{\varepsilon}^{1})_{i} &= a_{ij}(\mathbf{y}) \frac{\partial u_{\varepsilon}^{1}}{\partial x_{j}} = \left( a_{ij}(\mathbf{y}) + a_{ik}(\mathbf{y}) \frac{\partial \chi^{j}(\mathbf{y})}{\partial y_{k}} \right) \frac{\partial u_{0}}{\partial x_{j}} + \varepsilon a_{ij}(\mathbf{y}) \chi^{k}(\mathbf{y}) \frac{\partial^{2} u_{0}}{\partial x_{j} \partial x_{k}} \\ (56) \qquad &= a_{ij}^{0} \frac{\partial u_{0}}{\partial x_{j}} + g_{i}^{j}(\mathbf{y}) \frac{\partial u_{0}}{\partial x_{j}} + \varepsilon a_{ij}(\mathbf{y}) \chi^{k}(\mathbf{y}) \frac{\partial^{2} u_{0}}{\partial x_{j} \partial x_{k}}, \end{aligned} $$

where $g_i^j(\mathbf{y}) = a_{ij}(\mathbf{y}) + a_{ik}(\mathbf{y}) \frac{\partial \chi^j(\mathbf{y})}{\partial y_k} - a_{ij}^0$. We have from (8) that the vector fields $\mathbf{g}^j$ are solenoidal, i.e. $\nabla_y \cdot \mathbf{g}^j = 0$, as their i-th component is $g_i^j$. Hence, by Theorem 3.4 and Remark 3.11 from [13] there exists $\boldsymbol{\alpha}^j \in W_{per}^{1,q}(Y)^d$, $\nabla \cdot \boldsymbol{\alpha}^j = 0$ such that

$$ (57) \qquad \mathbf{g}^k = \mathrm{curl}_y \boldsymbol{\alpha}^k \text{ with } \|\boldsymbol{\alpha}^k\|_{1,q,Y} \le c \|\chi^k\|_{1,q,Y} \text{ with } q > d. $$

Equation (56) yields

$$ (58) \qquad \mathcal{A}_{\varepsilon} \nabla u_{\varepsilon}^{1} - \mathcal{A}_{0} \nabla u_{0} = \mathrm{curl}_{y} \boldsymbol{\alpha}^{k}(\mathbf{y}) \frac{\partial u_{0}}{\partial x_{j}} + \varepsilon \left( a_{ij}(\mathbf{y}) \chi^{k}(\mathbf{y}) \frac{\partial^{2} u_{0}}{\partial x_{j} \partial x_{k}} \right), $$

where the last term on the right-hand side of (58) is the vector whose i-th component is $a_{ij}(\mathbf{y}) \chi^k(\mathbf{y}) \frac{\partial^2 u_0}{\partial x_j \partial x_k}$. Next, we observe that

$$ \mathrm{curl} \left( \boldsymbol{\alpha}^k(\mathbf{y}) \frac{\partial u_0}{\partial x_k} \right) = \boldsymbol{\alpha}^k(\mathbf{y}) \times \nabla \frac{\partial u_0}{\partial x_k} + \mathrm{curl} \boldsymbol{\alpha}^k(\mathbf{y}) \frac{\partial u_0}{\partial x_k}. $$
Since $\mathrm{curl}_{\mathbf{y}}\boldsymbol{\alpha}^k(\mathbf{y})\frac{\partial u_0}{\partial x_j} = \varepsilon \mathrm{curl}\,\boldsymbol{\alpha}^k(\mathbf{x}/\varepsilon)\frac{\partial u_0}{\partial x_j}$ it follows from (58)

$$ (59) \qquad \begin{aligned} \mathcal{A}_{\varepsilon} \nabla u_{\varepsilon}^1 - \mathcal{A}_0 \nabla u_0 &= \varepsilon \mathrm{curl} \left( \boldsymbol{\alpha}^k(\mathbf{x}/\varepsilon) \frac{\partial u_0}{\partial x_k} \right) \\ &\quad - \varepsilon \boldsymbol{\alpha}^k(\mathbf{x}/\varepsilon) \times \nabla \frac{\partial u_0}{\partial x_k} + \varepsilon a_{ij}(\mathbf{x}/\varepsilon) \chi^k(\mathbf{x}/\varepsilon) \frac{\partial^2 u_0}{\partial x_j \partial x_k}. \end{aligned} $$

Next, from (57) and Sobolev inequalities we obtain

$$ (60) \qquad \begin{aligned} \left\| \varepsilon \boldsymbol{\alpha}^k (\mathbf{x}/\varepsilon) \times \nabla \frac{\partial u_0}{\partial x_k} + \varepsilon a_{ij} (\mathbf{x}/\varepsilon) \chi^k (\mathbf{x}/\varepsilon) \frac{\partial^2 u_0}{\partial x_j \partial x_k} \right\|_0 &\le \varepsilon \left( \|\boldsymbol{\alpha}^k\|_{0,\infty,Y} + \|\boldsymbol{\chi}\|_{0,\infty,Y} \right) \|u_0\|_2 \\ &\le c \varepsilon \|\boldsymbol{\chi}\|_{1,q,Y} \|u_0\|_2, \end{aligned} $$

where the constant $c$ depends only on $Y$. Thereby, induced by (59) and the estimate (60), we define the ansatz

$$ (61) \qquad \mathbf{z}_{\varepsilon}^{1} := \mathcal{A}_{0} \nabla u_{0} + \varepsilon \mathrm{curl} \left( \boldsymbol{\alpha}^{k}(\mathbf{x}/\varepsilon) \frac{\partial u_{0}}{\partial x_{k}} \tau_{\varepsilon} \right), $$

where $\tau_{\varepsilon} \in C_0^{\infty}(\Omega)$ is a cut-off function satisfying $\|\nabla \tau_{\varepsilon}\|_{\infty} \le \frac{c}{\varepsilon}$ and $\tau_{\varepsilon}(x) = 1$ if $\mathrm{dist}(\mathbf{x}, \partial\Omega) > \varepsilon$. Set

$$ \Omega_{\varepsilon} := \{ \mathbf{x} \in \Omega : \mathrm{dist}(\mathbf{x}, \partial\Omega) \le \varepsilon \}, $$

and observe that, from (59) and (60), the leading error between $\mathbf{z}_{\varepsilon}^1$ and $\mathcal{A}_{\varepsilon} \nabla u_{\varepsilon}^1$ can be bounded as follows

$$ (62) \qquad \left\| \varepsilon \mathrm{curl} \left( \boldsymbol{\alpha}^k(\mathbf{x}/\varepsilon) \frac{\partial u_0}{\partial x_j} (1-\tau_\varepsilon) \right) \right\|_0 \le \left\| \varepsilon (1-\tau_\varepsilon) \boldsymbol{\alpha}^k(\mathbf{x}/\varepsilon) \times \nabla \frac{\partial u_0}{\partial x_j} \right\|_{0,\Omega_\varepsilon} + \\ \qquad \left\| \varepsilon (1-\tau_\varepsilon) \frac{\partial u_0}{\partial x_j} \mathrm{curl} \boldsymbol{\alpha}^k(\mathbf{x}/\varepsilon) \right\|_{0,\Omega_\varepsilon} + \left\| \varepsilon \boldsymbol{\alpha}^k(\mathbf{x}/\varepsilon) \times \frac{\partial u_0}{\partial x_j} \nabla \tau_\varepsilon \right\|_{0,\Omega_\varepsilon}. $$

The first term on the right-hand side of (62) is bounded using Sobolev embedding Theorem as follows

$$ \left\| \varepsilon (1 - \tau_{\varepsilon}) \boldsymbol{\alpha}^k (\mathbf{x}/\varepsilon) \nabla \frac{\partial u_0}{\partial x_j} \right\|_{0,\Omega_{\varepsilon}} \le \varepsilon \|\boldsymbol{\chi}\|_{0,\infty,\Omega_{\varepsilon}} \|u_0\|_{2,\Omega_{\varepsilon}} \le \varepsilon \|\boldsymbol{\chi}\|_{1,q,Y} \|u_0\|_2. $$

Let $\alpha_j^k$ denotes the j-th element of $\boldsymbol{\alpha}^k$. To estimate the remaining terms in (62), we use the Y-periodicity of the function $\alpha_j^k$ to obtain

$$ \left( \int_{\Omega_{\varepsilon}} \alpha_j^k(\mathbf{x}/\varepsilon)^q d\mathbf{x} \right)^{1/q} \le \left( \frac{L_{\Omega}}{\varepsilon^{d-1}} \int_{\varepsilon Y} \alpha_j^k(\mathbf{x}/\varepsilon)^q d\mathbf{x} \right)^{1/q} \le \left( \varepsilon L_{\Omega} \int_{Y} \alpha_j^k(\mathbf{y})^q d\mathbf{y} \right)^{1/q}. $$
The second term on the right-hand side of (62) is estimated as follows

$$
\begin{aligned}
\left\| \varepsilon (1 - \tau^{\varepsilon}) \frac{\partial \alpha_j^k}{\partial x_l} \frac{\partial u_0}{\partial x_j} \right\|_{0, \Omega_{\varepsilon}} &\leq c(p') \|(1 - \tau^{\varepsilon})\|_{0,s, \Omega_{\varepsilon}} \left\| \varepsilon \frac{\partial \alpha_j^k}{\partial x_l} \right\|_{0,q, \Omega_{\varepsilon}} \left\| \frac{\partial u_0}{\partial x_j} \right\|_{0,p', \Omega_{\varepsilon}} \\
&\leq c(p') |\Omega_{\varepsilon}|^{\frac{1}{s}} (L_{\Omega} \varepsilon)^{\frac{1}{q}} |\Omega|^{1/p' - 1/2} \|\mathbf{x}\|_{1,q,Y} \|u_0\|_{2,\Omega} \\
&\leq c(p') (L_{\Omega} \varepsilon)^{1/2 - 1/p'} |\Omega|^{1/p' - 1/2} \|\mathbf{x}\|_{1,q,Y} \|u_0\|_{2,\Omega}
\end{aligned}
$$

where

$$ (63) \qquad \frac{1}{s} + \frac{1}{p'} + \frac{1}{q} = \frac{1}{2}. $$

The constant $c(p')$ depends on $p'$ and its dependence on $\Omega_{\varepsilon}$ relies only on the cone property of $\partial\Omega_{\varepsilon}$, and the term $|\Omega|^{1/p'-1/2}$ appears from the Sobolev embedding constant; [26]. As for the third term on the right-hand side of (62), we observe that

$$
\begin{aligned}
\left\| \varepsilon \frac{\partial \tau^{\varepsilon}}{\partial x_l} \alpha_l^k \frac{\partial u_0}{\partial x_j} \right\|_{0, \Omega_{\varepsilon}} &\leq \left\| \varepsilon \frac{\partial \tau^{\varepsilon}}{\partial x_l} \right\|_{0,s, \Omega_{\varepsilon}} \|\boldsymbol{\alpha}^k\|_{0,q, \Omega_{\varepsilon}} \|u_0\|_{1,p', \Omega_{\varepsilon}} \\
&\leq \left\| \varepsilon \frac{\partial \tau^{\varepsilon}}{\partial x_l} \right\|_{0, \infty, \Omega_{\varepsilon}} |\Omega_{\varepsilon}|^{\frac{1}{s}} (L_{\Omega} \varepsilon)^{\frac{1}{q}} \|\mathbf{x}\|_{0,q,Y} \left\| \frac{\partial u_0}{\partial x_j} \right\|_{0,p', \Omega_{\varepsilon}} \\
(64) \qquad &\leq c(p') (L_{\Omega} \varepsilon)^{1/2 - 1/p'} |\Omega|^{1/p' - 1/2} \|\mathbf{x}\|_{1,q,Y} \|u_0\|_{2,\Omega},
\end{aligned}
$$

where we used (7), and that $p'$ and $s$ satisfy (63). We now estimate the last term on the right-hand side of (62) as follows

$$ \left\| \varepsilon \tau^{\varepsilon} \chi_{\varepsilon}^{j} \frac{\partial^2 u_0}{\partial x_k \partial x_j} \right\|_{0, \Omega_{\varepsilon}} \leq \varepsilon \|\chi_{\varepsilon}^j\|_{0, \infty, \Omega_{\varepsilon}} \|u_0\|_{2, \Omega_{\varepsilon}} \leq \varepsilon \|\mathbf{x}\|_{1, q, Y} \|u_0\|_{2}. $$

Finally, gathering previous contributions, we conclude that

$$ (65) \qquad \varepsilon \left\| \operatorname{curl} \left( \boldsymbol{\alpha}^k(\boldsymbol{x}/\varepsilon) \frac{\partial u_0}{\partial x_k} (1-\tau_{\varepsilon}) \right) \right\|_0 \leq [c\varepsilon + c(p')(L_{\Omega}\varepsilon)^{1/2-1/p'} |\Omega|^{1/p'-1/2}] \|\boldsymbol{\chi}\|_{1,q,Y} \|u_0\|_2, $$

where $p'$ satisfy (63), and we used estimate (57). Hence, from the definition of $z_{\varepsilon}^1$ in (61), identity (59) and inequality (60), we prove (54) holds. Now, since $\tau_{\varepsilon} \in C_0^{\infty}(\Omega)$, we get

$$ \int_{\Omega} \varepsilon \operatorname{curl} \left( \boldsymbol{\alpha}^k(\boldsymbol{x}/\varepsilon) \frac{\partial u_0}{\partial x_k} \tau_{\varepsilon} \right) \cdot \nabla \phi \, d\boldsymbol{x} = 0 \quad \text{for all } \phi \in H^1(\Omega), $$

and as a result, from the definition of $z_{\varepsilon}^1$ in (61) and the weak formulations satisfied by $u_0$ and $u_{\varepsilon}$, it follows

$$ \int_{\Omega} (\mathbf{z}_{\varepsilon}^{1} - \mathcal{A}_{\varepsilon} \nabla u_{\varepsilon}) \cdot \nabla \phi \, d\mathbf{x} = \int_{\Omega} (\mathcal{A}_{0} \nabla u_{0} - \mathcal{A}_{\varepsilon} \nabla u_{\varepsilon}) \cdot \nabla \phi \, d\mathbf{x} = 0 \quad \text{for all } \phi \in H^{1}(\Omega), $$
and then (55) holds. It remains to proof that the estimate (55) also holds on the $L^2(\Omega)$ norm. Nevertheless, in the pure Neumann case, $u_{\epsilon}^1$ does not have mean value zero in $\Omega$ and neither vanishes at part of the boundary $\partial\Omega$. Consequently, we cannot directly apply Poincaré inequality to obtain the desired result. Thereby, we define

$$M_{u_1} := \frac{1}{|\Omega|} \int_{\Omega} \chi^j \partial_j u_0 \, d\boldsymbol{x},$$

and, using $u_{\epsilon}, u_0 \in L_0^2(\Omega)$ we have that $u_{\epsilon} - u_{\epsilon}^1 - \epsilon M_{u_1} \in L_0^2(\Omega)$, and applying Poincaré inequality, it holds

$$ (66) \quad \begin{aligned} \|u_{\epsilon} - u_{\epsilon}^1 - \epsilon M_{u_1}\|_0 &\le c d_{\Omega} |u_{\epsilon} - u_{\epsilon}^1|_1 \\ &\le c(p') d_{\Omega} \left[ c\epsilon + c(p') \left( \frac{L_{\Omega}\epsilon}{|\Omega|} \right)^{1/2-1/p'} \right] \|\boldsymbol{x}\|_{1,q,Y} \|u_0\|_2. \end{aligned} $$

To estimate $M_{u_1}$, we use Sobolev inequality which yields

$$ |M_{u_1}| \le |\Omega|^{-1} \|\chi^j\|_0 \|u_0\|_1 \le |\Omega|^{-1} |\Omega|^{1/2} \|\chi^j\|_{\infty,Y} \|u_0\|_1 \le c |\Omega|^{-1/2} \|\chi^j\|_{1,q,Y} \|u_0\|_1, $$

and, as a result, we arrive at

$$ \|M_{u_1}\|_0 \le c \|\chi^j\|_{1,q,Y} \|u_0\|_1. $$

From triangle inequality, we conclude

$$ \|u_{\epsilon} - u_{\epsilon}^{1}\|_{0} \le c(p') \left[ c\epsilon + c(p') \left( \frac{L_{\Omega}\epsilon}{|\Omega|} \right)^{1/2-1/p'} \right] (1 + d_{\Omega}) \|\boldsymbol{x}\|_{1,q,Y} \|u_{0}\|_{2}, $$

and the result follows.

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

## REFERENCES

[1] R. Araya, G. R. Barrenechea, L. P. Franca, and F. Valentin. Stabilization arising from PGEM: a review and further developments. *Appl. Numer. Math.*, 59(9):2065–2081, 2009.

[2] R. Araya, C. Harder, D. Paredes, and F. Valentin. Multiscale hybrid-mixed method. *SIAM J. Numer. Anal.*, 51(6):3505–3531, 2013.

[3] T. Arbogast, G. Pencheva, M. F. Wheeler, and I. Yotov. A multiscale mortar mixed finite element method. *SIAM Multiscale Modeling and Simulation*, 6:319–346, 2007.

[4] I. Babuška and E. Osborn. Generalized finite element methods: Their performance and their relation to mixed methods. *SIAM J. Num. Anal.*, 20(3):510–536, 1983.

[5] G.R. Barrenechea, F. Jaillet, D. Paredes, and F. Valentin. The multiscale hybrid mixed method in general polygonal meshes. *Numer. Math.*, 145:197–237, 2020.

[6] T. Chaumont-Frélet, A. Ern, S. Lemaire, and F. Valentin. Bridging the multiscale hybrid-mixed and multiscale hybrid high-order methods. *ESAIM: M2AN*, 56(1):261–285, 2022.
[7] T. Chaumont-Frelet, D. Paredes, and F. Valentin. Flux approximation on unfitted meshes and application to multiscale hybrid-mixed methods. Technical report, hal-03834748, 2022.

[8] M. Cicuttin, A. Ern, and S. Lemaire. A multiscale hybrid high-order method. *Comput. Methods Appl. Math*, 19:723-748, 2019.

[9] O. Duran, P. R. B Devloo, S. M. Gomes, and F. Valentin. A multiscale hybrid method for darcy's problems using mixed finite element local solvers. *Comp. Meth. Appl. Mech. Eng.*, 354:213-244, 2019.

[10] Y. R. Efendiev, T. Y. Hou, and Xiao-Hui Wu. Convergence of a nonconforming multiscale finite element method. *SIAM J. Numer. Anal.*, 37(3):888-910, 2000.

[11] A. Ern and J.-L. Guermond. *Theory and practice of finite elements*. Springer-Verlag, Berlin, New York, 2004.

[12] H. Fernando, C. Harder, D. Paredes, and F. Valentin. Numerical multiscale methods for a reaction-dominated model. *Computer Methods in Applied Mechanics and Engineering*, 201:228-244, 2012.

[13] V. Girault and P.-A. Raviart. *Finite element methods for Navier-Stokes equations*, volume 5 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1986. Theory and algorithms.

[14] A.T.A. Gomes, W. Pereira, and F. Valentin. The mhm method for linear elasticity on polytopal meshes. *to appear in IMA Journal of Numerical Analysis*, 2022.

[15] 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:107-130, 2013.

[16] C. Harder and F. Valentin. Foundations of the MHM method. In G. R. Barrenechea, F. Brezzi, A. Cangiani, and E. H. Georgoulis, editors, *Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations*, Lecture Notes in Computational Science and Engineering, Edinburgh, 2016. Springer.

[17] 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(227):913-943, 1999.

[18] T. J. R. Hughes, G. R. Feijó, L. Mazzei, and J. Quincy. The variational multiscale method - a paradigm for computational mechanics. *Computer Methods in Applied Mechanics and Engineering*, 166(1-2):3-24, 1998.

[19] V. V. Jikov, S. M. Kozlov, and O. A. Olešnik. *Homogenization of differential operators and integral functionals*. Springer-Verlag, Berlin, 1994. Translated from the Russian by G. A. Yosifian [G. A. Iosif'yan].

[20] A. L. Madureira and M. Sarkis. Hybrid localized spectral decomposition for multiscale problems. *SIAM J. Numer. Anal.*, 59(2):829-863, 2021.

[21] A. Malqvist and D. Peterseim. Localization of elliptic multiscale problems. *Math. Comp.*, 83(290):2583-2603, 2014.

[22] Society of Petroleum Engineers. Spe comparative solution project, 2000.

[23] D. Paredes, F. Valentin, and H. M. Versieux. On the robustness of multiscale hybrid-mixed methods. *Math. Comp.*, 86(304):525-548, 2017.

[24] P.-A. Raviart and J. M. Thomas. Primal hybrid finite element methods for 2nd order elliptic equations. *Math. Comp.*, 31(138):391-413, 1977.
[25] L. Scott and S. Zhang. Finite element interpolation of nonsmooth functions satisfying boundary conditions. *Mathematics of Computation*, 54:483–493, 1990.

[26] K. Tanaka, K. Sekine, M. Mizuguchi, and et al. Estimation of sobolev-type embedding constant on domains with minimally smooth boundary using extension operator. *J. Inequal. Appl.*, 389, 2015.

[27] E Weinan and B. Engquist. The heterogeneous multiscale methods. *Comm. Math. Sci.*, 1(1):87–132, 2003.

(*) DEPARTAMENTO DE INGENIERÍA MATEMÁTICA AND CI$^2$MA, UNIVERSIDAD DE CONCEPCIÓN, CHILE

E-mail address: dparedes@udec.cl

(†) LNCC - NATIONAL LABORATORY FOR SCIENTIFIC COMPUTING, AV. GETÚLIO VARGAS, 333, 25651-070 PETRÓPOLIS - RJ, BRAZIL AND CENTRE INRIA DE L'UNIVERSITÉ CÔTE D'AZUR, FRANCE

E-mail address: valentin@lncc.br; frederic.valentin@inria.fr

(‡) DEPARTAMENTO DE MATEMÁTICA, UNIVERSIDADE FEDERAL DE MINAS GERAIS - UFMG, BELO HORIZONTE - MG, BRAZIL

E-mail address: versieux@mat.ufmg.br