**1. Introduction**

The Darcy-Forchheimer equation is commonly used for describing the high velocity flow near oil and gas wellbores and fractures, which is a correction formula of the well known Darcy's law by supplementing a nonlinear velocity quantity as follows:

$$
\mu k^{-1} \mathfrak{u} + \beta \rho \left| \mathfrak{u} \right| \mathfrak{u} + \nabla p = 0,\tag{1}
$$

where *μ*, *k*, *ρ* and *β* represent the viscosity, the permeability, the density, and the dynamic viscosity coefficient of the fluid, respectively. *β* is also mentioned as the Forchheimer coefficient, whose values stand for the nonlinear intensity. In contrast, Darcy's law, which is valid for the extremely small velocity case, is usually used to show the linear relationship between the velocity vector *u* and the

pressure gradient ∇*p*. The Darcy-Forchheimer model can be obtained by coupling Equation (1) with the following conservation law equation:

$$
\nabla \cdot \mathfrak{w} = f.\tag{2}
$$

In recent years, the Darcy-Forchheimer model has been studied by many researchers within numerical discretized methods. Girault et al. in [1] proved the existence and uniqueness of the solution of the Darcy-Forchheimer model. Then, they considered mixed finite element methods by piecewise constant and nonconforming Crouzeix-Raviart elements to approximate the velocity and the pressure, respectively. Park in [2] gave a mixed finite element method (MFEM) for generalized Darcy-Forchheimer flow. Pan et al. in [3] presented an MFEM to approximate the velocity based on the Raviart–Thomas element or the Brezzi–Douglas–Marini element and piecewise constant for the pressure of the Darcy-Forchheimer model. Rui et al. in [4] published a block-centered finite difference method (BCFDM) for the Darcy-Forchheimer model. The authors in [5] established a BCFDM for the Darcy-Forchheimer model with variable parameter *β*(*x*). Rui and Liu in [6] introduced a two level BCFDM for the Darcy-Forchheimer model. Huang, Chen, and Rui in [7] designed a nonlinear multigrid method for the two-dimensional Darcy-Forchheimer model with a Peaceman–Rachford-type iteration as a smoother and proposed a better choice of the parameter used in the splitting. We point out that the constructed multigrid method with an almost linear computational cost is convergen<sup>t</sup> independent of the critical parameters.

To solve the problem in heterogeneous media, a very fine mesh should be used for solving a heterogeneity scale problem, which results to a large number of degrees of freedom. For dimension reduction and fast solvers of such problems, some model reduction techniques are needed. In [8], the authors introduced the mixed multiscale finite element methods and presented the main convergence results for the solution of second order elliptic equations with heterogeneous coefficients, which oscillate rapidly. Mixed multiscale finite element methods are widely used for solving the reservoir simulation problems [9–11]. In [12–14], we developed a mixed generalized multiscale finite element method (GMsFEM), where we enriched a multiscale space by new degrees of freedom, which was obtained by solving local spectral problems on the snapshot space.

In this paper, we consider a solution of the Darcy-Forchheimer model in heterogeneous media. We construct an efficient algorithm of the mixed generalized multiscale finite element method to make an approximation of the problem on the coarse grid. The construction is based on solving the local problems for calculating multiscale basis functions of the velocity. Meanwhile, we use piecewise constant elements to approximate the pressure. In the mixed formulation, we firstly define a snapshot space, which provides a solution space in each local area. Then we solve the local spectral problem in the snapshot space to find out multiscale basis functions. We use the Picard iteration to address the non-linearity when solving the problems on the multiscale spaces [15,16]. Note that when constructing multiscale basis functions, we do not take into account the nonlinear part of the Darcy-Forchheimer equation.

The remainder of this article is organized as follows: The model problem and its weak formulation in mixed form and the discrete weak formulation are demonstrated in Section 2. The fine grid approximation and coarse gird approximation are presented in Section 3 and Section 4, respectively. Some numerical experiments using our mixed generalized multiscale finite element method are carried out in Section 5 to verify the efficiency of the presented method. Finally, conclusions and further ideas are presented in Section 6.

## **2. Mathematical Model**

We consider the steady state Darcy-Forchheimer model to describe a single phase fluid flow in a heterogeneous porous medium.

$$\begin{aligned} \mu k^{-1}(\mathbf{x})\mathfrak{u} + \beta(\mathbf{x})\rho|\mathfrak{u}|\mathfrak{u} + \nabla p &= 0, \quad \mathbf{x} \in \Omega, \\ \nabla \cdot \mathfrak{u} &= f, \quad \mathbf{x} \in \Omega, \end{aligned} \tag{3}$$

where Ω ∈ R<sup>2</sup> is a bounded domain and *∂*Ω is Lipschitz continuous; |*v*| = (*<sup>v</sup>*, *v*) 12 , (·, ·) denotes the *L*<sup>2</sup> inner product; *k* and *β* are the heterogeneous permeability and the heterogeneous non-Darcy coefficient, respectively.

Complemented with Dirichlet boundary condition,

$$p = f\_{D, \prime} \quad \mathfrak{x} \in \partial \Omega \tag{4}$$

Without loss of generality, we suppose that *fD* = 0, namely the homogeneous Dirichlet boundary condition. Then, we ge<sup>t</sup> the following problem:

$$\begin{aligned} \mu k^{-1}(\mathbf{x})\boldsymbol{\mu} + \beta(\mathbf{x})\rho|\boldsymbol{u}|\boldsymbol{u} + \nabla p &= 0, \quad \mathbf{x} \in \Omega, \\ \nabla \cdot \boldsymbol{\mu} &= \boldsymbol{f}, \quad \mathbf{x} \in \Omega, \\ \boldsymbol{p} &= 0, \quad \mathbf{x} \in \partial\Omega. \end{aligned} \tag{5}$$

**Remark 1.** *The boundary condition* (4) *can be replaced by the Neumann boundary condition,*

$$
\mathfrak{u} \cdot \mathfrak{n} = f\_{\mathcal{N}}, \quad \mathfrak{x} \in \partial \Omega,\tag{6}
$$

*f and fN should satisfy the compatibility condition by the Gauss theorem,*

$$\int\_{\Omega} f\left(\mathbf{x}\right) \, \mathrm{d}\mathbf{x} = \int\_{\partial\Omega} f\_N\left(s\right) \, \mathrm{d}s.\tag{7}$$

#### **3. Fine Grid Approximation**

In this section, mixed finite element methods have been borrowed to handle Problem (5) in the fine grid. Here, we use the standard Sobolev spaces notation to define the function spaces and their norms as follows:

$$V = \left\{ \mathbf{v} \in L^3(\Omega)^2; \nabla \cdot \mathbf{v} \in L^2(\Omega) \right\}, \quad \|\mathbf{u}\|\_V = \|\mathbf{v}\|\_{0,3\Omega} + \|\nabla \cdot \mathbf{v}\|\_{0,2,\Omega}, \forall \mathbf{v} \in V.$$

$$Q = L^2(\Omega), \quad \|q\|\_{Q} = \|q\|\_{0,2,\Omega}, \forall q \in Q.$$

Then, we obtain the following variational formulation: find (*<sup>u</sup>*, *p*) ∈ *V* × *Q* such that:

$$\begin{split} \int\_{\Omega} \mu k^{-1}(\mathbf{x}) \boldsymbol{\mu} \, \mathbf{v} \, \mathbf{d} \mathbf{x} + \int\_{\Omega} \beta(\mathbf{x}) \rho |\boldsymbol{\mu}| \boldsymbol{\mu} \, \mathbf{v} \, \mathbf{d} \mathbf{x} - \int\_{\Omega} p \, \nabla \cdot \boldsymbol{\mathbf{v}} \, \mathbf{d} \mathbf{x} = \mathbf{0}, \quad \forall \mathbf{v} \in V, \\ -\int\_{\Omega} q \, \nabla \cdot \mathbf{u} \, \mathbf{d} \mathbf{x} = -\int\_{\Omega} f \, q \, \mathbf{d} \mathbf{x}, \quad \forall q \in Q. \end{split} \tag{8}$$

The variational formulation (8) and the problem (5) are equivalent by the following Green's formula:

$$
\int\_{\Omega} \nabla p \, \mathbf{v} \, \mathbf{d} \mathbf{x} = -\int\_{\Omega} p \, \nabla \cdot \mathbf{v} \, \mathbf{d} \mathbf{x} + \int\_{\partial \Omega} p \, \mathbf{v} \cdot \mathbf{n} \, \mathbf{d} s, \quad \forall \mathbf{v} \in V. \tag{9}
$$

Let Ω be a polygon in two dimensions, which can be entirely covered by a shape regular decomposition T*<sup>h</sup>*, in the sense of Ciarlet [17], into triangles, with *h* being the maximum diameter of the elements of the triangles. Therefore, T*h* is a family of conforming triangulations of Ω,

$$\overline{\Omega} = \bigcup\_{\mathbb{T} \in \mathfrak{T}\_h} \mathcal{T}.$$

Here, we discretize the velocity *u* in the lowest order Raviart–Thomas space: Given a simplex *T* ∈ R2, the local Raviart–Thomas space of order *k* ≥ 0 is defined by:

$$\mathcal{R}T\_k(T) = \mathcal{P}\_k(T)^2 + \mathbf{x}\,\mathcal{P}\_k(T)\_{\prime\prime}$$

where P*k* denotes the set of all polynomials in two variables of degree ≤ *k*. Then, we can ge<sup>t</sup> the lowest order global Raviart–Thomas space, which is the conforming finite element space of the velocity *u*:

$$V\_{\hbar} = RT\_0 \left( \Omega \right) = \left\{ \boldsymbol{\sigma} \in V \, : \, \boldsymbol{\sigma} \vert\_{T} \in RT\_0(T) \quad \forall \, \boldsymbol{T} \in \mathcal{T}\_{\hbar} \right\}.\tag{10}$$

The pressure *p* is approximated in the following piecewise constant space:

$$Q\_{\rm li} = \left\{ q \in L^2(\Omega) \, : \, q|\_T \in \mathcal{P}\_0 \quad \forall \, T \in \mathcal{T}\_{\rm h} \right\}. \tag{11}$$

Then, we can obtain the discrete weak formulation of (8): find a pair (*<sup>u</sup>h*, *ph*) ∈ *Vh* × *Qh*:

$$\begin{split} \int\_{\Omega} \mu \mathbf{k}^{-1} (\mathbf{x}) \boldsymbol{\mu}\_{\hbar} \boldsymbol{\nu}\_{\hbar} \, \mathbf{dx} + \int\_{\Omega} \beta (\mathbf{x}) \rho |\boldsymbol{\mu}\_{\hbar}| \boldsymbol{\mu}\_{\hbar} \boldsymbol{\nu}\_{\hbar} \, \mathbf{dx} - \int\_{\Omega} p\_{\hbar} \nabla \cdot \boldsymbol{\boldsymbol{\nu}}\_{\hbar} \, \mathbf{dx} &= \mathbf{0}, \quad \forall \boldsymbol{\nu}\_{\hbar} \in V\_{\hbar} \\ - \int\_{\Omega} q\_{\hbar} \nabla \cdot \boldsymbol{\boldsymbol{\nu}}\_{\hbar} \, \mathbf{dx} &= - \int\_{\Omega} f \, q\_{\hbar} \, \mathbf{dx}, \quad \forall q\_{\hbar} \in Q\_{\hbar}. \end{split} \tag{12}$$

In [3], the authors demonstrated the existence and uniqueness of the continuous problem and the discrete problem, respectively. Moreover, if T*h* is quasi-uniform and (*<sup>u</sup>*, *p*) ∈ *Ws*,<sup>3</sup>(Ω)<sup>2</sup> × *Ws*, 3 2 (Ω), then the following error estimates can be proven; see ([3], Theorem 4.4) for details:

$$\|\|\boldsymbol{u} - \boldsymbol{u}\_{h}\|\|\_{0,2}^{2} + \|\|\boldsymbol{u} - \boldsymbol{u}\_{h}\|\|\_{0,3}^{3} \le Ch^{2s}, \quad 1 \le s \le k+1,\tag{13}$$

$$||p - p\_h||\_{0,2} \le Ch^s, \quad 1 \le s \le k+1. \tag{14}$$

where *Ws*,*<sup>p</sup>*(Ω) is the standard Sobolev space of index (*s*, *p*), where *s* is a nonnegative integer and *p* ≥ 1.

$$\mathcal{W}^{s,p}(\Omega) = \{ v \in L^p(\Omega) : D^a v \in L^p(\Omega) \text{ for all } |a| \le s \},$$

where *D<sup>α</sup>v* is the weak derivative of *v*; see the details in [18].

Let:

$$\mathfrak{u}\_{\hbar} = \sum\_{i=1}^{m\_1} \mathfrak{u}\_i \mathfrak{f}\_{i\prime} \quad p\_{\hbar} = \sum\_{i=1}^{m\_2} p\_i \theta\_{i\prime}$$

where *u* = [*<sup>u</sup>*1, *u*2, ... , *um*1 ] *T*, *p* = [*p*1, *p*2, ... , *pm*2 ] *T* are the coefficients of the finite element approximations with *m*1, *m*2 dimensions, in terms of a basis *ξ* of the velocity and a basis *θ* of the pressure, respectively.

We apply the Picard iteration for solving the resulting discrete nonlinear system. Find *un*+<sup>1</sup> *h* ∈ *Vh*, *p<sup>n</sup>*+<sup>1</sup> *h* ∈ *Qh* with an arbitrary initial guess *u*0 *h* ∈ *Vh*, such that:

$$\begin{split} \int\_{\Omega} \mu \mathbf{k}^{-1}(\mathbf{x}) \boldsymbol{\upmu}\_{h}^{n+1} \boldsymbol{\upmu}\_{h} \, \mathrm{d}\mathbf{x} + \int\_{\Omega} \beta(\mathbf{x}) \rho |\mathbf{u}\_{h}^{n}| \boldsymbol{\upmu}\_{h}^{n+1} \boldsymbol{\upnu}\_{h} \, \mathrm{d}\mathbf{x} - \int\_{\Omega} p\_{h}^{n+1} \, \nabla \cdot \boldsymbol{\upnu}\_{h} \, \mathrm{d}\mathbf{x} = \boldsymbol{0}, \quad \forall \boldsymbol{\upnu}\_{h} \in V\_{\mathrm{h}}. \\ - \int\_{\Omega} q\_{h} \, \nabla \cdot \mathbf{u}\_{h}^{n+1} \, \mathrm{d}\mathbf{x} = - \int\_{\Omega} f \, q\_{h} \, \mathrm{d}\mathbf{x}, \quad \forall q\_{h} \in Q\_{\mathrm{h}}. \end{split} \tag{15}$$

We rewrite the iteration (15) into the following matrix form.

$$
\begin{pmatrix} A\left(\mathfrak{u}\_h^n\right) & B^T \\ B & 0 \end{pmatrix} \begin{pmatrix} \mathfrak{u}\_h^{n+1} \\ p\_h^{n+1} \end{pmatrix} = \begin{pmatrix} 0 \\ \mathcal{F} \end{pmatrix},\tag{16}
$$

where *<sup>A</sup>*(*<sup>u</sup>h*) is the matrix associated with the term:

$$\int\_{\Omega} \mu k^{-1}(\mathbf{x}) \mu\_{\hbar} \boldsymbol{\sigma}\_{\hbar} \, \mathrm{d}\mathbf{x} + \int\_{\Omega} \beta(\mathbf{x}) \rho |\boldsymbol{\mu}\_{\hbar}| \boldsymbol{\mu}\_{\hbar} \, \boldsymbol{\sigma}\_{\hbar} \, \mathrm{d}\mathbf{x},$$

where *B* is the matrix corresponding to − 'Ω *qh* ∇ · *uh* d*<sup>x</sup>*, and *F* is the vector associated with the linear functional − 'Ω *f qh* d*<sup>x</sup>*.

In the practical implementation, we use the following termination criterion to control the iteration,

$$\max(r\_{\mathfrak{u}\_{\prime}}, r\_{\mathcal{P}}) < tol\_{\prime}$$

where:

$$\begin{array}{rcl} \mathbf{r}\_{\mathsf{u}} &=& \left\| \mu \mathbf{k}^{-1} (\mathbf{x}) \boldsymbol{\mu}\_{\mathsf{h}}^{\mathsf{u}+1} + \beta(\mathbf{x}) \rho |\boldsymbol{\mu}\_{\mathsf{h}}^{\mathsf{u}+1}| \boldsymbol{\mu}\_{\mathsf{h}}^{\mathsf{u}+1} + \nabla p\_{\mathsf{h}}^{\mathsf{u}+1} \right\|\_{0^{\prime}} \\\\ \mathbf{r}\_{\mathsf{p}} &=& \left\| \begin{array}{c} \left\| f - \nabla \cdot \boldsymbol{\mu}\_{\mathsf{h}}^{\mathsf{u}+1} \right\|\_{0^{\prime}} / \left\| f \right\|\_{0^{\prime}} & \text{ when } \left\| f \right\|\_{0} \neq 0, \\\left\| f - \nabla \cdot \boldsymbol{\mu}\_{\mathsf{h}}^{\mathsf{u}+1} \right\|\_{0^{\prime}} & \text{ when } \left\| f \right\|\_{0} = 0. \end{array} \right. \end{array}$$

#### **4. Coarse Grid Approximation**

In this section, we describe the construction of the approximation on a coarse grid using the mixed generalized multiscale finite element method (mixed GMsFEM). We construct a square uniform coarse grid T*H* for the computational domain Ω with a coarse grid size *H*; E*H* = ∪*NE <sup>i</sup>*=<sup>1</sup>*Ei* is the set of all facets of a coarse mesh, and *NE* is the number of facets of a coarse mesh. To construct the multiscale basis function, we build a uniform triangular fine grid, which is obtained by refinement of a coarse grid. In the mixed GMsFEM, we compute the multiscale basis functions for the velocity in the local domains *ωi* that correspond to the coarse edge *Ei* ∈ E*H* (Figure 1).

$$
\omega\_i = \cup\_j \{ \mathcal{K}\_j \in \mathcal{T}\_H | E\_i \subset \partial \mathcal{K}\_j \},
$$

where *Kj* is the coarse grid cell, which is equal to a square element of T*<sup>H</sup>*.

We build a multiscale space for the velocity *ums* ∈ *Vms*:

$$V\_{ms} = \text{span}\{\Psi\_1, \dots, \Psi\_{M\_{\omega\_l} \cdot N\_E}\}\_\prime \tag{17}$$

where *ψi* are the multiscale basis functions, which are calculated in the local domain *ωi* and *M<sup>ω</sup>i* is the number of basis functions in each local domain. For pressure, we use the space *Qms* of piecewise constant functions on the coarse cells.

We start with constructing a snapshot space in the local domain *ωi* and then make an approximation on the coarse grid using the solution of the local spectral problem on the snapshot space. The snapshot space is obtained by solving the next local problem in *ωi*: find (*φj*, *η*) ∈ *V<sup>ω</sup><sup>i</sup> h* × *Q<sup>ω</sup>i h* such that:

$$\begin{split} \int\_{\omega\_{i}} k^{-1} \boldsymbol{\Phi}\_{\rangle}^{i} \boldsymbol{\nu} \, \mathrm{d}\mathbf{x} - \int\_{\omega\_{i}} \eta \, \nabla \cdot \boldsymbol{\nu} \, \mathrm{d}\mathbf{x} &= \mathbf{0}, \quad \boldsymbol{\sigma} \in V\_{h}^{\omega\_{i}} \, \mathrm{} \\ \int\_{\omega\_{i}} r \, \nabla \cdot \boldsymbol{\Phi}\_{\rangle}^{i} \mathrm{d}\mathbf{x} &= \int\_{\omega\_{i}} c\_{j} \, r \, \mathrm{d}\mathbf{x}, \quad r \in \bigotimes\_{h}^{\omega\_{i}} \, \mathrm{} \end{split} \tag{18}$$

with boundary condition:

$$
\boldsymbol{\phi}\_{\dot{j}}^{\dot{i}} \cdot \boldsymbol{\mathfrak{n}} = 0, \quad \boldsymbol{\mathfrak{x}} \in \partial \boldsymbol{\omega}\_{\dot{i}\prime} \tag{19}
$$

where *n* is the unit exterior normal vector to the boundary *∂ωi*.

On the fine facet *ej*, we apply an additional boundary condition:

$$
\boldsymbol{\Phi}^l\_j \cdot \boldsymbol{\mathfrak{n}} = \boldsymbol{\delta}\_{\dot{\boldsymbol{i}}j}, \quad \boldsymbol{\mathfrak{x}} \in \boldsymbol{\mathfrak{e}}\_{\dot{\boldsymbol{\lambda}}} \tag{20}
$$

where *j* = 1, ...,*J<sup>ω</sup><sup>i</sup>* is the number of the fine facets *ej* on which we define the boundary condition (20), *cj* = |*ej*| *<sup>S</sup><sup>ω</sup>i* is chosen by the compatibility condition '*∂ωi φij* · *n* = ' ∇ · *φij*, |*ej*| presents the length of the fine grid facet *ej*, and *S<sup>ω</sup>i* denotes the volume of the local domain *ωi*. Here, *J<sup>ω</sup>i* is the number of fine grid edges *ej* on *Ei*, *Ei* = ∪*J<sup>ω</sup><sup>i</sup> j*=<sup>1</sup>*ej*, and *δij* is a piecewise constant function defined on *Ei*, which takes the value of one on *ej* and zero on each other edge of the fine grid.

**Figure 1.** Illustration of the heterogeneous property, coarse grid, and multiscale basis functions in the local domain.

Therefore, we can define the snapshot space in each local domain *ωi*:

$$V\_{suap}^{\dot{i}} = \text{span}\{\phi\_1^{\dot{i}}, \dots, \phi\_{\|i\|}^{\dot{i}}\}.$$

To compute a multiscale basis function, we solve a local spectral problem in each local domain *ωi*. The spectral problem allows us to find the most important characteristics of the problem. Then, we use multiscale basis functions to define the multiscale space to ge<sup>t</sup> an approximation on the coarse grid. In each *ωi*, we solve the spectral problem on *<sup>V</sup>isnap*:

$$
\tilde{A}^{\omega\_i} \tilde{\Psi}\_l^{\omega\_i} = \lambda\_l \tilde{S}^{\omega\_i} \tilde{\Psi}\_l^{\omega\_i}, \quad \tilde{A}^{\omega\_i} = R^{\omega\_i} A^{\omega\_i} (R^{\omega\_i})^T, \quad \tilde{S}^{\omega\_i} = R^{\omega\_i} S^{\omega\_i} (R^{\omega\_i})^T,\tag{21}
$$

where:

$$(\mathcal{R}^{\omega\_i})^T = [\Phi^i\_1, \dots, \Phi^i\_{\parallel^{\omega\_i}}]\_{\prime\prime}$$

and:

$$A^{\omega\_i} = [a^{\omega\_i}\_{mn}], \quad a^{\omega\_i}\_{mn} = a\_{\omega\_i} (\Phi\_m \cdot \Phi\_n) = \int\_{E\_i} k^{-1} (\Phi\_m \cdot \mathfrak{n}) (\Phi\_n \cdot \mathfrak{n}) \, ds,\tag{22}$$

$$S^{\omega\_i} = [s^{\omega\_l}\_{nm}], \quad s^{\omega\_l}\_{mn} = s\_{\omega\_i} (\boldsymbol{\Phi}\_{\mathcal{W}}, \boldsymbol{\Phi}\_{\mathcal{n}}) = \int\_{\omega\_i} k^{-1} \boldsymbol{\Phi}\_{\mathcal{W}} \boldsymbol{\Phi}\_{\mathcal{n}} \, d\mathbf{x} + \int\_{\omega\_i} \nabla \cdot \boldsymbol{\Phi}\_{\mathcal{n}} \nabla \cdot \boldsymbol{\Phi}\_{\mathcal{n}} \, d\mathbf{x}.\tag{23}$$

For construction of the multiscale space, we select the first *M<sup>ω</sup>i* smallest eigenvalues and take the corresponding eigenvectors *ψ<sup>ω</sup><sup>i</sup> l* = (*R<sup>ω</sup>i*)*<sup>T</sup>ψ*˜ *<sup>ω</sup>i l* as basis functions (*l* = 1, 2, ... , *<sup>M</sup><sup>ω</sup>i*). We note that the presented spectral decomposition in the snapshot space is motivated by theoretical analysis (see Theorem 4.3 in [12]). Moreover, the oversampling techniques can be used for the construction of the multiscale basis functions to enhance the accuracy of mixed GMsFEM. The main idea of oversampling techniques is to introduce a small dimensional snapshot space using the POD (proper orthogonal decomposition) approach, where snapshot vectors are constructed in larger regions that contain the interfaces of two adjacent coarse blocks [12].

**Remark 2.** *For the construction of the multiscale basis functions, we can use other types of spectral problems (see [12,19]). For example, we can solve the following spectral problem:*

$$\vec{S}^{\omega\_{\overline{\imath}}} \vec{\psi}\_{l}^{\omega\_{\overline{\imath}}} = \mu\_{l} \vec{\Psi}\_{l}^{\omega\_{\overline{\imath}}}, \quad \vec{S}^{\omega\_{\overline{\imath}}} = R^{\omega\_{\overline{\imath}}} S^{\omega\_{\overline{\imath}}} (R^{\omega\_{\overline{\imath}}})^T,\tag{24}$$

*with:*

$$S^{\omega\_i} = [s\_{mn}^{\omega\_i}], \quad s\_{mn}^{\omega\_i} = s\_{\omega\_i}(\Phi\_{m\nu}\Phi\_n) = \int\_{\omega\_i} k^{-1} \Phi\_m \Phi\_n \,d\mathbf{x},\tag{25}$$

*on the snapshot space and taking eigenvectors corresponding to the largest eigenvalues as a multiscale basis functions (λl* = 1/*μl).*

Additionally, we calculate the first basis function by the solution of the following local problem in the domain *ωi*: find (*χ<sup>ω</sup><sup>i</sup>* , *η*) ∈ *V<sup>ω</sup><sup>i</sup> h*× *Q<sup>ω</sup>i h*such that:

$$\begin{split} \int\_{\omega\_{i}} k^{-1} \mathfrak{X}^{\omega\_{i}} \mathfrak{v} \, d\mathfrak{x} - \int\_{\omega\_{i}} \eta \, \nabla \cdot \mathfrak{v} \, d\mathfrak{x} &= 0, \quad \mathfrak{v} \in V\_{\hbar}^{\omega\_{i}}, \\ \int\_{\omega\_{i}} r \, \nabla \cdot \mathfrak{X}^{\omega\_{i}} \, d\mathfrak{x} &= \int\_{\omega\_{i}} \mathfrak{c} \, r \, d\mathfrak{x}, \quad r \in Q\_{\hbar}^{\omega\_{i}}. \end{split} \tag{26}$$

with the following boundary conditions:

$$\mathbb{X}^{\omega\_{i}} \cdot \mathfrak{n} = 0, \quad \mathfrak{x} \in \partial \omega\_{i}, \quad \mathbb{X}^{\omega\_{i}} \cdot \mathfrak{n} = 1, \quad \mathfrak{x} \in E\_{i}, \tag{27}$$

where *c* = |*Ei*| *<sup>S</sup><sup>ω</sup>i*, |*Ei*| is the length of coarse facet *Ei*.

 Finally, we obtain the following multiscale space for velocity:

$$V\_{\rm ms} = \text{span}\{\chi^{\omega\_i}, \Psi\_l^{\omega\_i}, 1 \le l \le M\_{\omega\_l}, 1 \le i \le N\_E\}. \tag{28}$$

The illustration of the local multiscale basis functions with an additional basis are presented in Figure 1.

For the construction of the coarse grid system, we define a projection matrix:

$$R = \begin{bmatrix} R\_{\mu} & 0 \\ 0 & R\_{p} \end{bmatrix}, \quad R\_{\mu} = \begin{bmatrix} R\_{\mu,1}, \dots, R\_{\mu,N\_{E}} \end{bmatrix}^{T}, \tag{29}$$

where (*Ru*,*<sup>i</sup>*)*<sup>T</sup>* = [*χ<sup>ω</sup><sup>i</sup>* , *ψ<sup>ω</sup><sup>i</sup>* 1 , ..., *ψ<sup>ω</sup><sup>i</sup> <sup>M</sup><sup>ω</sup>i* ] and *Rp* is the projection matrix for pressure, where we set one for each fine grid cell in the current coarse grid cell. Here, *NE* is the number of facets of the coarse grid, and *M<sup>ω</sup>i* is the number of multiscale basis functions in local domain *ωi*.

We use the constructed multiscale space and write the approximation in the coarse grid in matrix form: 

$$
\begin{pmatrix} A\_c^k & B\_c^T \\ B\_c & 0 \end{pmatrix} \begin{pmatrix} \mathfrak{u}\_c^{k+1} \\ p\_c^{k+1} \end{pmatrix} = \begin{pmatrix} 0 \\ F\_c \end{pmatrix},\tag{30}
$$

where:

$$A\_c^k = R\_{\text{ll}} A^k R\_{\text{ll}}^T, \quad B\_{\text{\textell}} = R\_{\text{ll}} B R\_{p\text{\textell}}^T, \quad F\_{\text{\textell}} = R\_{p} F\_{\text{\textell}} \tag{31}$$

and finally, we reconstruct the solution on the fine grid *uk*+<sup>1</sup> *ms* = *RTu uk*+<sup>1</sup> *c* .

## **5. Numerical Results**

In this section, we present several numerical results of the Darcy-Forchheimer model in the domain Ω = [0, 1]2. We used a structured 160 × 160 fine grid with 77,120 edges and 51,200 cells (triangles) and a 10 × 10 coarse grid with 220 edges and 100 cells (see Figure 2).

**Figure 2.** Coarse grid (**left**) and fine grid (**right**).

In the numerical study, we set *μ* = 1, *ρ* = 1, and *f* = 1. For the Picard iteration, we set *tol* = 10−8. The permeability tensor *k* is heterogeneous and presented in Figure 3, where we considered two test cases. We set the Darcy-Forchheimer coefficient *β* = *C* · *k*−<sup>1</sup> from [20,21], where the parameter *C* controls the influence of the nonlinear part of the equation. We studied the proposed multiscale solver for *C* = 10.24, *C* = 34.93, *C* = 1584.14, and *C* = 71,554.17. With such values of *C*, we could investigate the behavior of a method with the various influences of the nonlinear part. By increasing the parameter *C*, we obtained a Darcy-Forchheimer equation with the dominant nonlinear part.

**Figure 3.** Heterogeneous coefficient for Test 1 (**left**) and Test 2 (**right**).

At first, we considered a test case with *β* = 0. Fine scale and multiscale solutions using eight multiscale basis functions are presented in Figures 4 and 5 for coefficient *k* from Test 1 and Test 2. In Table 1, the relative errors in the *L*2 norm for different numbers of multiscale basis functions are presented for *β* = 0. Here, *DOFc* and *DOFf* denote the size of the multiscale and fine grid solutions (*DOFf* = 128,320), and *M* is the number of multiscale basis functions. By #iter, we denote the number of Picard iterations, and *eu* and *ep* are the relative errors in *L*2 norm for velocity and pressure, respectively. To calculate the error of the pressure, the average values over the coarse grid cells are used.

**Figure 4.** Fine grid solution (**top**) and multiscale solution using eight multiscale basis functions (**bottom**). Coefficient *k* from Test 1 with *β* = 0.

**Figure 5.** Fine grid solution (**top**) and multiscale solution using eight multiscale basis functions (**bottom**). Coefficient *k* from Test 2 with *C* = 34.93.

The numerical results for *C* = 34.93 with coefficient *k* from Test 1 and Test 2 using eight multiscale basis functions are presented in Figures 6 and 7. In Tables 2–5, we show the relative error in the *L*2 norm between the multiscale solution and the fine grid solution for different numbers of multiscale basis functions. The errors are presented for different values of *C* to see the influence of the nonlinear part on the accuracy of mixed GMsFEM. According to the obtained results, we observed that the method worked well with the presented problem. When we increased the number of multiscale bases, we observed that the error decreased. For large values of *C* the error was greater than for smaller values. This was due to the fact that for large *C* values, the influence of the nonlinear part of the equation increased, and our multiscale bases did not take into account the nonlinear part. From the

tables, we can observe that the number of Picard iterations of mixed GMsFEM was significantly smaller than the number of iterations of the fine grid solution for large *C*. For Test 1 and Test 2, we obtained similar results for the presented multiscale solver.

**Figure 6.** Fine grid solution (**top**) and multiscale solution using eight multiscale basis functions (**bottom**). Coefficient *k* from Test 1 with *C* = 34.93.

**Figure 7.** Fine grid solution (**top**) and multiscale solution using eight multiscale basis functions (**bottom**). Coefficient *k* from Test 2 with *β* = 0.

**Table 1.** Relative error in the *L*2 norm for different numbers of multiscale basis functions. Coefficient *k* from Test 1 (left) and Test 2 (right) with *β* = 0.



**Table 2.** Relative error in the *L*2 norm for different numbers of multiscale basis functions. Coefficient *k* from Test 1 (left) and Test 2 (right) with *C* = 10.24.

**Table 3.** Relative error in the *L*2 norm for different numbers of multiscale basis functions. Coefficient *k* from Test 1 (left) and Test 2 (right) with *C* = 34.93.


**Table 4.** Relative error in the *L*2 norm for different numbers of multiscale basis functions. Coefficient *k* from Test 1 (left) and Test 2 (right) with *C* = 1581.14.


**Table 5.** Relative error in the *L*2 norm for different numbers of multiscale basis functions. Coefficient *k* from Test 1 (left) and Test 2 (right) with *C* = 71,554.17.

