A hybrid high-order method for a coupled Stokes-Darcy problem on general meshes

Journal Pre-proof A Hybrid High-Order Method for a Coupled Stokes-Darcy Problem on General Meshes

Yongchao Zhang, Liquan Mei, Rui Li

PII:

S0021-9991(19)30769-7

DOI:

https://doi.org/10.1016/j.jcp.2019.109064

Reference:

YJCPH 109064

To appear in:

Journal of Computational Physics

Received date:

4 April 2018

Revised date:

28 September 2019

Accepted date:

23 October 2019

Please cite this article as: Y. Zhang et al., A Hybrid High-Order Method for a Coupled Stokes-Darcy Problem on General Meshes, J. Comput. Phys. (2019), 109064, doi: https://doi.org/10.1016/j.jcp.2019.109064.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Highlights • • • • •

A hybrid high-order (HHO) method for a coupled Stokes-Darcy problem is presented. The all primal variables are approximated by polynomials of arbitrary degree k >= 0. We obtain the optimal (high-order) error estimates in energy-norm. HHO is convenient to implement, mass conserving, robust for physical parameters. Numerical simulation of flow in the free fluid region and fractured porous media.

A Hybrid High-Order Method for a Coupled Stokes-Darcy Problem on General Meshes $ Yongchao Zhanga , Liquan Meia,∗, Rui Lib a School

b School

of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, P. R. China of Mathematics and Information Science, Shaanxi Normal University, Xi’an, Shaanxi 710062, P. R. China

Abstract In this work, a hybrid high-order (HHO) method on general meshes is presented to solve a coupled StokesDarcy problem with the Beavers-Joseph-Saffman interface condition. Constructed on polynomials of arbitrary degree k ≥ 0, the numerical method is established in terms of discrete unknowns attached to mesh faces and cells (or elements). The unified discrete scheme for Stokes equation and Darcy equation is given by the continuity condition of the interface. The unique solvability of the discrete scheme is proved. Moreover, the energy error estimate for the velocity and L2 -error estimate for pressure of order (k + 1) are derived. Finally, a series of numerical experiments are reported to illustrate the accuracy, mass conservation and robustness of our method. Keywords: Stokes-Darcy; Hybrid high-order method; Beavers-Joseph-Saffman condition; Error estimates; General meshes

1. Introduction Over the past decades, a model of coupled fluid flow and porous media flow has received increasing attention, this type of coupled model arises from many applications, for instance the interaction of ground and surface water [1, 2], blood motion in the vessels [3] and engineering filtration problems [4, 5]. In this paper, we are mainly interested in the Stokes-Darcy problem, which describes the coupling of the fluid flow (governed by the Stokes equations) and the porous media flow (governed by the Darcy’s law) through three transmission conditions at the interface. A wide variety of numerical methods have been developed to solve the Stokes-Darcy problem, and a nonexhaustive list of these methods includes finite element method [6, 7, 8], domain decomposition method [9, 10], Lagrange multiplier method [11, 12], multigrid method [13, 14], discontinuous Galerkin method [15, 16], discontinuous finite volume element method [17, 18], mortar finite element method [19, 20], least square method [21, 22], boundary integral method [23, 24] and spectral method [25]. However, with the fact that given the complexity of geometric resulting from compaction, erosion, and the onset of fractures or fault, the meshes cannot meet the needs of the standard methods. In recent years, numerical methods on general meshes, possibly with nonconforming and polytopal (polygonal/polyhedral) elements, have witnessed a great development. Such novel approaches include hybrid high-order (HHO) method [26, 27, 28], virtual element method (VEM) [29, 30], discontinuous Galerkin method [31, 32], hybridizable discontinuous Galerkin (HDG) method [33, 34], weak Galerkin (WG) method [35, 36], mixed and hybrid finite volume (MHFV) method [37, 38], and the family of mimetic finite difference (MFD) methods [39, 40], and so on. Among them, MHFV and MFD methods, which are mainly concerned lowest-order schemes, have been extensively studied and widely applied to different types of problems. Indeed, the algebraic equivalence between the MHFV and the MFD methods have been demonstrated in [41]. While, the high-order polytopal discretization methods, e.g., HHO, VEM, DG, HDG and WG methods, have received significant attention over the last few years. Nevertheless, as far as we know, there are few study of the Stokes-Darcy problem by using these methods. In fact, we always hope that the proposed numerical methods can be applied to practical problems. Considering the complexity of the Stokes-Darcy problem and in order to obtain good numerical results, it is required that the numerical methods can handle the interface conditions simply and $ This work is supported by Science Challenge Project of China (TZ2016002) and a scholarship from the China Scholarship Council (CSC, 201806280141). ∗ Corresponding author. Email addresses: [email protected] (Yongchao Zhang), [email protected] (Liquan Mei ), [email protected] (Rui Li)

conveniently, adapt to complex geometries, maintain mass conservation, and be robust to physical parameters. Recently, the VEM [42], WG [43, 44] and DG [45] finite element discretizations, on general meshes, have been proposed to solve the Stokes-Darcy problem. Our focus is here on the HHO method, which is originally derived for the scalar and mixed diffusion problems [26, 28] as well as the linear elasticity problems [27]. Nowadays, this framework has been successfully extended to the discretization of advection-dominated transport problem [46], Stokes problem [47], Brinkman model [48] non-linear Leray-Lions problems [49, 50] and Cahn-Hilliard equation [51]. It is defined hybrid because both element- and face-based discrete unknowns are used. Meanwhile, the term high-order referring to the fact that the order of the polynomial functions can be an arbitrary integer k ≥ 0, which allows arbitrary orders of convergence of the error estimates. Compared to DG method, the HHO method does not need the penalty constant. What’s more, thanks to the compactness techniques, HHO method is also appealing in terms of computational cost, this is similar to HDG method. Links between HDG, VEM and HHO methods have been pointed out in [52, 53, 54]. For smooth solutions, another remarkable feature of HHO method is that the corresponding discrete schemes lead to energy error estimates of order (k + 1) and super-closeness of order (k + 2) in L2 -error estimates. In the following, we consider a HHO method for the Stokes-Darcy problem. Benefiting from the property of supporting general meshes, we illustrate an easy treatment of non-matching meshes on interface (c.f., Remark 2.1). The prior work [48] has studied the HHO discretisation of Brinkman problem, which is also a form of combination of the Stokes and Darcy models. While, in our paper, by using the Stokes equations in primal velocity-pressure formulation, and the mixed formulation for Darcy’s law, we derive the existence and uniqueness of the discrete solutions. We analyze error estimates of order (k+1) in the energy-norm for velocity and L2 -norm for pressure. Besides, in our numerical implementations, we show the strong inclusiveness of the HHO method for polynomial degrees (c.f., Example 5.2). Unlike the construction of arbitrary higher-order polynomials in Example 5.2, in Example 5.3, we will see that for the lowest order case (k = 0), which shares strong links with the MHFV method [26], the method still maintains a good convergence order. This is different from the recent works. It is worth mentioned that the HHO method for Stokes-Darcy problem holds a natural treatment of physical parameters [55] and keeps robust over the entire range of variation of hydraulic conductivity in Darcy region (c.f., Example 5.3−5.5). All the above indicate that the HHO method is a simple and robust finite element method that has flexible mesh generation, good accuracy, robustness properties for the Stokes-Darcy problem. The rest of the paper is organized as follows. In Section 2, we recall and state the coupled Stokes-Darcy model and its weak formulation. In Section 3, the fundamental definitions and HHO method discrete scheme are proposed. The solvability of the discrete problem in Section 4 and the optimal error estimates in suitable norms are proved. A series of numerical experiments are presented to validate the theoretical results in Section 5. Finally, we end the paper with some conclusions in Section 6. 2. Preliminaries 2.1. Model problem We consider the model problem in an open, bounded and convex polygonal or polyhedral domain Ω ⊂ Rd (d = 2 or 3), consisting of a fluid subregion ΩS and a porous medium subregion ΩD . Assume that both ΩS and ΩD have Lipschitz continuous boundaries. For any subset X ⊂ ΩS or ΩD , the space L2 (X) is equipped with the standard L2 -scalar product (·, ·)X and L2 -norm ·X := ·L2 (X) . Similarly, the Sobolev space H k (X) = W k,2 (X) is defined in the usual way with the norm ·k,X and seminorm | · |k,X . In the following, we often use bold fonts to express vector variables and spaces, such as, L2 (X) := L2 (X)d , H01 (X) := H01 (X)d , H k (X) := H k (X)d . For further use, we allow X to denote the discrete element T or face F (see Section 2.2 for details), then (·, ·)X , ·X , ·k,X and | · |k,X have the same definitions as above. Denote by Γ = ∂ΩS ∩ ∂ΩD the sharp interface between the Stokes and Darcy regions, and by Γi = ∂Ωi \Γ for i = {S, D} the outer boundary, as illustrated in Figure 1. Let nS and nD be the unit outward normal vectors on ∂ΩS and ∂ΩD , respectively. Especially, we have nS = −nD on the interface Γ. To this end, we define the unit normal vector of Γ as nS , as shown in Figure 1. For convenience, we denote by u = (uS , uD ) the velocity and p = (pS , pD ) the pressure, where, ui = u|Ωi

2

Figure 1: A sketch of the free fluid domain ΩS , porous media domain ΩD and interface Γ.

and pi = p|Ωi , i = {S, D}. In ΩS , the free fluid flow is assumed to be governed by the Stokes equations: −νΔuS + ∇pS = fS ∇ · uS = 0 S uS = gD

in ΩS ,

(2.1a)

in ΩS ,

(2.1b)

on ΓS ,

(2.1c)

where ν is the kinematic viscosity, uS is the fluid velocity, pS is the kinematic pressure, and fS denotes a general body force term that includes gravitational acceleration. In ΩD , the porous medium flow motion is governed by Darcy’s law: K−1 uD = −∇pD ∇ · uD = fD uD · n =

D gN

in ΩD ,

(2.2a)

in ΩD ,

(2.2b)

on ΓD ,

(2.2c)

where K is the hydraulic conductivity tensor, uD is the specific discharge rate in the porous medium, pD is the hydraulic head, and fD is a sink/source term. Furthermore, we assume that K is piecewise constant and takes symmetric positive definite values with eigenvalues λ in the interval 0 < λmin ≤ λ ≤ λmax < +∞. For S = 0 on ΓS and simplicity, we set zero Dirichlet and Neumann boundary conditions for uS and uD , i.e., gD D gN = 0 on ΓD , respectively. For the coupled Stokes-Darcy model, the key part is the interface conditions that describe how different types of flow interact at the fluid/porous medium interface Γ. Significantly, the interface conditions of mass conservation (2.3a), forces balance (2.3b) and the classical Beavers-Joseph-Saffman condition [56, 57] are imposed as follows: uS · nf = uD · nf , ∂uS = pD , pS − νnf · ∂nf √ ∂uS α ν τj · uS , j = 1, ..., d − 1, = −ντj · ∂nf trace(K)

(2.3a) (2.3b) (2.3c)

where α is the Beavers-Joseph constant depending on the properties of the porous medium, τj (j = 1, ..., d − 1) is the unit tangent vectors on Γ, as depicted in Figure 1. We introduce the following two function spaces for the velocity and the pressure, V := {v ∈ L2 (Ω) : ∇ · v ∈ L2 (Ω), v|ΩS ∈ H 1 (ΩS ), v = 0 on ΓS , v · nD = 0 on ΓD },  2 2 Q := L0 (Ω) = {q ∈ L (Ω) : q = 0}. Ω

By using the interface conditions (2.3) and integration by parts, the weak formulation of coupled Stokes-Darcy

3

model reads as: Find u ∈ V and p ∈ Q such that a(u, v) − b(v, p) = (fS , v)ΩS b(u, q) = (fD , q)ΩD

∀v ∈ V ,

(2.4a)

∀q ∈ Q,

(2.4b)

with the bilinear forms are defined as a(u, v) = aS (u, v) + aD (u, v) + aΓ (u, v), b(v, q) = (∇ · v, q)Ω , where aS (u, v) = ν(∇u, ∇v)ΩS , aD (u, v) = (K−1 u, v)ΩD , √ α ν (Pτ (u), Pτ (v))Γ . aΓ (u, v) =  trace(K) d−1 Here, Pτ is the projection onto the tangent space on Γ defined as Pτ (v) = j=1 (v · τj )τj , and (·, ·)Γ denotes the L2 inner product on the interface Γ. The well-posedness of the problem (2.4) can be found in [6]. 2.2. Meshes and basic results In this subsection, we refer to the settings of meshes by [47, 58]. Let H ⊂ R+ ∗ be a countable set of meshsizes having 0 as its unique accumulation point. We consider an h-refined mesh sequence (Th )h∈H where, for all h ∈ H, Th = {T } is a matching simplicial mesh characterized by the scalar h := maxT ∈Th hT (hT is the diameter of the element T ). We assume that (Th )h∈H is regular in the sense of [58]. The set of internal faces is denoted by Fhi , the set of boundary faces is denoted by Fhb , and it holds that Fh := Fhi ∪ Fhb . For all F ∈ Fh , TF := {T ∈ Th : ∂T ∩ F = F }, the diameter of F is denoted by hF . For all T ∈ Th , FT denotes the set of faces of T , for all F ∈ FT , nT F is the unit normal to F pointing out of T , and nF is a prescribed normal direction for F . Inheriting the above properties, we let Th be the decomposition of Ω into polygons or polyhedrons which are aligned with the interface Γ, and Fh be the collection of all the faces. Since any T ∈ Th cannot be cut by Γ, the partition Th can be naturally divided into two sets denoted by ThS = Th ∩ ΩS and ThD = Th ∩ ΩD . Moreover, let FhS and FhD be the set of all faces of ThS and ThD , respectively, and denote FhS,b := FhS ∩ (∂ΩS \Γ), FhD,b := FhD ∩ (∂ΩD \Γ) and FhΓ := FhS ∩ FhD . Furthermore, all the meshes in ThD are assumed to be compatible with the known partition on which the diffusion tensor is piecewise constant. The restriction of K to a mesh cell T ∈ ThD is denoted KT . The lowest and largest eigenvalue of KT are denoted by λmin,T and λmax,T , respectively, and the local anisotropy ratio ρKT := λmax,T /λmin,T . Throughout the paper, we use the abbreviation A  B for the inequality A ≤ cB, where the letter c denotes a positive constant independent of the mesh size h, but may depend on the domain Ω (or subdomains ΩS , ΩD ) and mesh regular parameter  (see [58, Definition 1.38] for details). Remark 2.1. We point out that the non-matching discrete meshes on Γ are allowed in practice, such as the situation illustrated in Figure 2. In this case, no special treatment is required, each face containing hanging nodes is treated as multiple coplanar faces. We use the same trick to treat the nonconforming meshes in ΩS and ΩD . What’s more, for the curved interface, we replace the original interface with short straight line segments (or flat planes), which coincide with FhΓ . Therefore, in each discrete element T ∈ Th , the faces are straight lines or flat planes. The above considered cases are examined by the numerical Example 5.1−Example 5.4. In the following, we recall some basic results [58, Section 1.4]. Let l ≥ 0 be a nonnegative integer. For an n-dimensional subset X of Ω (n ≤ d), Pl (X) is the space spanned by the restrictions to X of n-variate polynomials of total degree ≤ l, the set X will represent a mesh element or face. There exists a real number ctr depending on  and l, but independent of h, the following discrete trace inequality holds for all T ∈ Th , [58, Lemma 1.46] −1/2

vF ≤ ctr hF

vT

∀v ∈ Pl (T ), F ∈ FT .

4

(2.5)

Figure 2: Treatment of non-matching interface discretizations, the interface Γ is colored in red. Gray elements are squares and blue elements are pentagons.

And for all T ∈ Th , there exists a constant cinv depending on  and l, but independent of h ([58, Lemma 1.44]), such that ∇vT ≤ cinv h−1 T vT

∀v ∈ Pl (T ).

(2.6)

l : L1 (X) → Pl (X) the L2 -orthogonal projector such that, for all v ∈ L1 (X), Furthermore, we denote πX  l (v − πX v)w = 0 ∀w ∈ Pl (X). X

Then, using [58, Lemma 1.40], there exists a real number capp depending on  and l, but independent of h, such that, for all T ∈ Th , s ∈ {1, ..., l + 1} and v ∈ H s (T ), 1/2

|v − πTl v|m,T + hT |v − πTl v|m,∂T ≤ capp hs−m |v|s,T T

∀m ∈ {0, ...(s − 1)}.

(2.7)

l , are obtained by Remark 2.2. The vector- and matrix-valued L2 -orthogonal projectors, both denoted by πX l applying πX component-wise, and we have the similar approximation properties as above.

3. Hybrid high-order discretization In this section, we first introduce the local discrete spaces and formulate reconstructions of different operators appearing in (2.4). Subsequently, the global discrete spaces and problems are established. 3.1. Local reconstructions Let the polynomials of degree k ≥ 0 be fixed. For each T ∈ Th , we define the local hybrid discrete space for the velocity variable, V kT :={v T := (vT , (vF )F ∈FT ) such that vT ∈ Pk (T )d , vF ∈ Pk (F )d for T ∈ ThS , F ∈ FT ; vT ∈ KT ∇Pk (T ), vF = vF nF where vF ∈ Pk (F ) for T ∈ ThD , F ∈ FT , especially, for T ∈ ThD and F ∈ FT ∩ FhΓ , we define vF ∈ Pk (F )d }.

(3.1)

For the pressure variable, we define the local discrete space as QkT := Pk (T ).

(3.2)

Now we present the local reconstructions based on [26, 28, 59]. Definition 3.1. (Local discrete divergence reconstruction). For each element T ∈ Th , the local divergence operator DTk := V kT → QkT is such that, for all v T ∈ V kT , q ∈ QkT ,   k  DT v T , q T = − (vT , ∇q)T + (vF , qnT F )F , F ∈FT

= (∇ · vT , q)T +



F ∈FT

5

(vF − vT , qnT F )F .

(3.3)

Definition 3.2. (Local discrete gradient reconstruction). Let T ∈ ThS be fixed, define the local gradient reconstruction operator rTk+1 : V kT → Pk+1 (T )d such that for all v T , ∈ V kT , w ∈ Pk+1 (T )d ,  (vF − vT , ∇w nT F )F , (3.4) (∇rTk+1 v T , ∇w)T = (∇vT , ∇w)T + F ∈FT

which further satisfies (rTk+1 v T − vT , 1)T = 0. Definition 3.3. (Local discrete flux reconstruction). Let T ∈ ThD be fixed, the consistent flux reconstruction operator FTk+1 : V kT → KT ∇Pk+1 (T ) is such that, for all v T ∈ V kT , w ∈ Pk+1 (T ), FTk+1 v T solves  (vF · nT F , w)F . (3.5) (FTk+1 v T , ∇w)T = −(DTk v T , w)T + F ∈FT

We define the local interpolation operator I kT : H 1 (T ) → V kT such that, for all v ∈ H 1 (T ), I kT v :={(ITk v, (IFk v)F ∈FT ) such that ITk v := πTk v, IFk v := πFk v for T ∈ ThS , F ∈ FT ; k ITk v := πK v, IFk v := πFk (v · nF )nF for T ∈ ThD , F ∈ FT , T

especially, for T ∈ ThD and F ∈ FT ∩ FhΓ , we define IFk v := πFk v},

(3.6)

k : H 1 (T ) → KT ∇Pk (T ) is the projector such that where πK T k v − v, ∇w)T = 0 (πK T

∀w ∈ Pk (T ).

In the following, we present some technical lemmas for the operators DTk , rTk+1 , FTk+1 and I kT , which have been reported in [26, 47, 59]. Lemma 3.4. (Commuting property). For all T ∈ Th , v ∈ H 1 (T ), the following holds: DTk (I kT v) = πTk (∇ · v).

(3.7)

Proof. For all T ∈ ThS , the validity of the conclusion has been verified in [47], so here we only focus on T ∈ ThD . For all T ∈ ThD and v ∈ H 1 (T ), by taking w ∈ Pk (T ), we obtain  (v · nT F , w)F . (πTk (∇ · v), w)T = (∇ · v, w)T = −(v, ∇w)T + F ∈FT k k With the definition of πK and I kT , we get −(v, ∇w)T = −(πK v, ∇w)T = −(ITk v, ∇w)T . Then, for the second T T term in the right-hand side of above equation, we have the following cases:

(v · nT F , w)F = (πFk v, w nT F )F = (IFk v · nT F , w)F for F ∈ FT ∩ FhΓ , (v · nT F , w)F = (πFk (v · nT F ), w)F = (πFk (v · nT F )nF · nF , w)F = (πFk (v · nF )nF · nT F , w)F = (IFk v · nT F , w)F for F ∈ FT but F ∈ / FT ∩ FhΓ , where we have use the definition of IFk on T ∈ ThD and the fact that nF = nT F or nF = −nT F . Therefore, by the Definition 3.1, it holds that  (πTk (∇ · v), w)T = −(ITk v, ∇w)T + (IFk v · nT F , w)F = (DTK (I kT v), w)T , F ∈FT

which completes the proof. Lemma 3.5. (Approximation property [47]). For all T ∈ ThS , v ∈ H k+2 (T ), the following optimal approximation property holds:  1/2 v − rTk+1 I kT vT + hT v − rTk+1 I kT vF + hT ∇(v − rTk+1 I kT v)T +

 F ∈FT

F ∈FT

3/2

hT ∇(v − rTk+1 I kT v)F  hk+2 T vk+2,T .

6

(3.8)

Lemma 3.6. (Polynomial preservation). For all T ∈ ThD , v ∈ KT ∇Pk+1 (T ), the following holds: FTk+1 (I kT v) = v.

(3.9)

Proof. For all T ∈ ThD , let v ∈ KT ∇Pk+1 (T ) and use the definition (3.5) of FTk+1 , we have that, for all w ∈ Pk+1 (T ),  (IFk v · nT F , w)F . (FTk+1 (I kT v), ∇w)T = −(DTk (I kT v), w)T + F ∈FT

Using Lemma 3.4 and v ∈ KT ∇Pk+1 (T ) ⊂ Pk (T )d , we infer that (DTk (I kT v), w)T = (πTk (∇ · v), w)T = (∇ · v, w)T , (IFk v · nT F , w)F = (πFk v · nT F , w)F = (v · nT F , w)F for F ∈ FT ∩ FhΓ , (IFk v · nT F , w)F = ((πFk (v · nF )nF ) · nT F , w)F = ((v · nT F )nF · nF , w)F = (v · nT F , w)F for F ∈ FT but F ∈ / FT ∩ FhΓ . Then (FTk+1 (I kT v), ∇w)T = −(∇ · v, w)T +



(v · nT F , w)F = (v, ∇w)T ,

F ∈FT

which gives the conclusion. 3.2. Local discretization To define local discrete bilinear forms for the interface term, we first introduce the notations: S D TΓF := ThS ∩ TF , TΓF := ThD ∩ TF for all F ∈ FhΓ .

aΓF

k k Then, using the Definition 3.1−3.3, the local discrete bilinear forms aST : V kT ×V kT → R, aD T : V T ×V T → R, k k k k : V T S × V T S → R, bT : V T × QT → R are defined as follows ΓF

ΓF

aST (wT , v T ) = ν(∇rTk+1 wT , ∇rTk+1 v T )T + νsST (wT , v T ) aD T (w T , v T )

=

K+1 (K−1 wT , FTk+1 v T )T T FT

aΓF (wT , v T ) = 

√

α ν D ) trace(KTΓF

bT (v T , qT ) = (DTk v T , qT )T

+

sD T (w T , v T )

(Pτ (wF ), Pτ (vF ))F

∀T ∈ ThS ,

(3.10)

ThD ,

(3.11)

∀T ∈

S ∀F ∈ FhΓ , T ∈ TΓF ,

∀T ∈ Th ,

(3.12) (3.13)

where the stabilization terms sST , sD T are defined by  k Tk+1 wT ), πFk (vF − r Tk+1 v T ))F , h−1 sST (wT , v T ) = F (πF (wF − r

(3.14)

F ∈FT

sD T (w T , v T ) =



F ∈FT

k+1 hF κ−1 wT − wF ) · nT F , (FTk+1 v T − vF ) · nT F )F , T F ((FT

(3.15)

with r Tk+1 v T := vT + (rTk+1 v T − πTk rTk+1 v T ), and κT F := KT nT F · nT F . Note that the stabilization terms have the following properties (c.f., e.g., [26] and [28]): 1

sST (I kT v, I kT v) /2  hr+1 T vr+2,T k sD T (I T w, v T )

= 0 ∀T ∈

ThD ,

∀T ∈ ThS , v ∈ H k+2 (T ), r ∈ {0, ..., k},

w ∈ KT ∇Pk+1 (T ), v T ∈ V

7

k T.

(3.16) (3.17)

3.3. Global discretization Here we give the following global discrete velocity space V kh and pressure space Qkh ,

V kh := v h = ((vT )T ∈Th , (vF )F ∈Fh ) such that v h |T = v T ∈ V kT for T ∈ Th ,  vF = 0 for F ∈ Fhb ,

(3.18)

Qkh :=   QkT .

(3.19)

T ∈Th

What’s more, we define the global interpolation operator I kh : H 1 (Ω) → V kh such that (I kh v)|T = I kT (v|T )

∀T ∈ Th , v ∈ H 1 (Ω).

Then from (3.10)−(3.19), the HHO discrete scheme for the coupled Stokes-Darcy flow can be written as follows: Find uh ∈ V kh and ph ∈ Qkh such that ah (uh , v h ) − bh (v h , ph ) = lhS (v h ) bh (uh , qh ) = where

lhD (qh )

∀v h ∈ V kh ,

(3.20a)

∀qh ∈

(3.20b)



 aST (uT , v T ) + aD aΓF (uTΓF , v TΓF ), T (uT , v T ) +

ah (uh , v h ) =

T ∈Th



bh (v h , qh ) =

Qkh ,

F ∈FhΓ

bT (v T , qT ),

T ∈Th

lhS (v h ) =



(fS , vT )T ,

T ∈ThS

lhD (qh ) =



(fD , qT )T .

T ∈ThD

4. Solvability and error analysis In this section, we study the solvability and stability of HHO discrete scheme for the Stokes-Darcy problem. We first define some useful norms for all v h ∈ V kh , qh ∈ Qkh and then recall several lemmas in subsection 4.1. Next, in subsections 4.3−4.2, the consistency error and the error estimates in the energy- and L2 -norm are proved. In the first place, we introduce the following discrete energy-norm on V kh ,

1/2 v h  := v h 2S,h + v h 2D,h + v h 2Γ,h ,

(4.1)

where v h 2S,h =



νv T 2S,T ,

v T 2S,T := ∇vT 2T +

F ∈FT

T ∈ThS

v h 2D,h =



v T 2D,T ,

v T 2D,T := vT 2T +

=



F ∈FhΓ



2 h−1 F vF − vT F ,

(4.2)

hF vF · nT F 2F ,

(4.3)

Pτ (vF )2F .

(4.4)

F ∈FT

T ∈ThD

v h 2Γ,h



2 v TΓF S Γ,F ,

2 v TΓF S Γ,F

The L2 -norm for the qh on Qkh is defined by qh 2 :=



:= 

√ α ν

D ) trace(KTΓF

qT 2T = qh 2Ω .

T ∈Th

8

(4.5)

4.1. Solvability To start with the proof of the solvability, we introduce two lemmas for norms equivalence, and then present that ah (·, ·) is continuous and coercive, and bh (·, ·) satisfies the discrete inf-sup condition. As a consequence, from the standard framework for saddle point problem [60], the discrete problem (3.20) is solvable. Lemma 4.1. ([26, Lemma 4]). There exists ζ > 0 such that: For all T ∈ ThS , v T ∈ V kT , ζνv T 2S,T ≤ aST (v T , v T ) ≤ ζ −1 νv T 2S,T .

(4.6)

Lemma 4.2. ([28, Lemma 4, Lemma 8]). There is η > 0, uniform with respect to the mesh size and the diffusion tensor, such that 2 D −1 −1 ηλ−1 λmin,T v T 2D,T , max,T v T D,T ≤ aT (v T , v T ) ≤ η

hT DTk v T T

1/2

− 1/2

+ λmin,T KT

FTk+1 v T T

+

1/2

1/2 λmin,T sD T (v T , v T )

(4.7)  v T D,T ,

(4.8)

for all T ∈ ThD and all v T ∈ V kT . In the rest of this subsection, we shall prove the solvability of (3.20). First, we have the following lemmas for ah (·, ·) and bh (·, ·). Lemma 4.3. For all v h , wh ∈ V kh , the bilinear form ah (·, ·) has following properties |ah (v h , wh )| ≤ acon  v h  wh , 2

ah (v h , v h ) ≥ acoe  v h  ,

(4.9) (4.10)

−1 where acon = ccon max{ζ −1 , maxT ∈ThD {η −1 λ−1 min,T }}, acoe = ccoe min{ζ, minT ∈ThD {ηλmax,T }}, and ccon , ccoe are constants independent of h.

Proof. In combination with Lemma 4.1 and Lemma 4.2 and the definition (4.1) of the norm  · , we can directly get (4.10). And then using the Cauchy-Schwarz inequality we have (4.9). Lemma 4.4. There exists a positive constant β, independent of h, such that: For all qh ∈ Qkh , βqh  ≤

sup v h ∈V k h ,v h ≤1

bh (v h , qh ).

(4.11)

Proof. Using the inf-sup condition for the continuous problem, we infer that there has a constant β0 > 0 such that, for each qh ∈ Qkh ⊂ L2 (Ω), there exists v ∈ H01 (Ω) with ∇vΩ ≤ 1 and  β0 qh  ≤ (∇ · v, qh )T . T ∈Th

Multiplying the above inequality by a constant c0 > 0, and choosing v h := I kh (c0 v), we get  (DTk (I kT (c0 v)), qh )T = bh (v h , qh ), c0 β0 qh  ≤

(4.12)

T ∈Th

where we have used the commuting property (3.7) and c0 > 0 will be chosen such that v h  ≤ 1. In the following of the proof, we only need to prove the existence of such an h-independent constant c0 . Recalling the definition of  · , we derive that I kh v  I kh vS,h + I kh vD,h + I kh vΓ,h := T1 + T2 + T3 . From [47, Lemma 3], we have T1 ≤ c1 ∇vΩS , where c1 is an h-independent constant. To estimate T2 , by the definition (3.6) of I kT for T ∈ ThD , we infer that    k (T2 )2  πK v2T + hF πFk v2F := T2,1 + T2,2 . T T ∈ThD

T ∈ThD F ∈FT

9

k v = KT ∇mT for some mT ∈ Pk (T ) and have Using the definition (3.6) of I kT , we write πK T −1 −1 −1 k k 2 k k λ−1 max,T πKT vT ≤ |(KT KT ∇mT , πKT v)T | = |(KT πKT v, v)T | ≤ λmin,T πKT vT vT , k vT ≤ ρKT vT , and T2,1 ≤ maxT ∈ThD {ρ2KT }v2ΩD . so we get πK T Then, using (2.7) and trace inequality (2.5), we derive

T2,2 ≤





T ∈ThD F ∈FT



  F ∈FT

T ∈ThD





1/2

1/2

hF πFk v − πTk vF + hF πTk vF 1/2

hF v − πTk vF +

 F ∈FT

capp hT ∇vT + ctr πTk vT

1/2

2

hF πTk vF

2

2

T ∈ThD

c2app h2 ∇v2ΩD + c2tr c2stab v2ΩD , where cstab is a constant only independent of h. So there holds (T2 )2  (capp h∇vΩD + (ctr cstab + maxT ∈ThD {ρKT })vΩD )2 . From v ∈ H01 (Ω), we have D vΩD ≤ cD p ∇vΩD , here, cp is a constant that only depends on ΩD . Consequently, we obtain T2  c2 ∇vΩD with c2 := (capp h + (ctr cstab + max {ρKT })cD p ), T ∈ThD

note that, c2 contains one high order term capp h which can be bounded when h → 0. For the T3 part in the v h  norm, we use the trace inequality to get √ √  α ν α ν 2 k 2  Pτ (πF v)F  max {  }v2Γ  c3 ∇v2ΩS , (T3 ) = F ∈FhΓ D D trace(K ) trace(K ) Γ TΓF TΓF F ∈F h

√ α ν } trace(KT D )

where c3 := cStr maxF ∈FhΓ { 

and cStr is the constant that only depends on ΩS .

ΓF

Summing over the above inequalities, we arrive v h   c4 ∇vΩ with c4 = max{c1 , c2 , Therefore, by choosing c0 = the conclusion of (4.11).

1 c4 ,

√

c3 }.

we have I kh (c0 v)  ∇vΩ ≤ 1. Finally, we let β = β0 c0 and this gives

From Lemma 4.3 and Lemma 4.4, the HHO scheme (3.20) has a unique solution. 4.2. Consistency error In this subsection, we present some useful lemmas which will be used in the error analysis of the HHO discretization. Inspired by [61], we introduce the coupled global space Zh := V kh × Qkh and the coupled global bilinear form Ah : Zh × Zh → R such that: For all (wh , rh ), (v h , qh ) ∈ Zh , Ah ((wh , rh ), (v h , qh )) := ah (wh , v h ) − bh (v h , rh ) + bh (wh , qh ).

(4.13)

Then, the discrete scheme of (3.20) can be written as follows: Find (uh , ph ) ∈ Zh such that Ah ((uh , ph ), (v h , qh )) = lhS (v h ) + lhD (qh )

∀(v h , qh ) ∈ Zh .

(4.14)

To get the a priori error estimate, we need the stability like [61, Definition 4]. Firstly, we define the norm on Zh : For all (v h , qh ) ∈ Zh , (v h , qh )2Zh := v h 2 +qh 2 . Subsequently, we have the following stability of bilinear form Ah : 10

(4.15)



− 1/2 −2 2 Lemma 4.5. (Stability of Ah ). There exists γ = a−2 acon )2 +4β −2 such that: For all (wh , rh ) ∈ coe (1+2β Zh , γ(wh , rh )Zh ≤

Ah ((wh , rh ), (v h , qh )) , (v h , qh )Zh (v h ,qh )∈Zh \{0} sup

(4.16)

where acon , acoe and β are defined in Lemma 4.3 and Lemma 4.4. Proof. From the definition (4.13) of Ah and the coercivity (4.10) of ah , we have −1 −1 ∗ wh 2 ≤ a−1 coe ah (w h , w h ) = acoe Ah ((w h , rh ), (w h , rh )) ≤ acoe A (w h , rh )Zh ,

(4.17)

where A∗ is the supremum of the right-hand side of (4.16). Using the definition (4.13) of Ah , the continuity (4.9) of ah and the inf-sup condition (4.11) of bh , we infer that rh  ≤ β −1 = β −1

sup v h ∈V k h ,v h ≤1

sup v h ∈V k h ,v h ≤1

≤β

−1

sup v h ∈V k h ,v h ≤1

bh (v h , rh )

ah (wh , v h ) − Ah ((wh , rh ), (v h , 0))

|ah (wh , v h ) − Ah ((wh , rh ), (v h , 0))|

≤ β −1 (acon  wh  +A∗ ). Hence, combing above two inequalities and using the Cauchy-Schwarz inequality, it holds (wh , rh )2Zh =  wh 2 +rh 2 ∗ −2 2 ≤a−1 acon  wh 2 +2β −2 A∗2 coe A (w h , rh )Zh + 2β −2 2 ≤a−1 acon )A∗ (wh , rh )Zh + 2β −2 A∗2 coe (1 + 2β 1 1 ≤ a−2 (1 + 2β −2 a2con )2 A∗2 + (wh , rh )2Zh + 2β −2 A∗2 , 2 coe 2

thus, we have



−2 2 2 −2 (1 + 2β a ) + 4β A∗2 . (wh , rh )2Zh ≤ a−2 coe con



− 1/2 −2 2 By taking γ = a−2 acon )2 + 4β −2 , we complete the proof. coe (1 + 2β We now introduce the consistency error Eh,u,p (v h ) for the solution u, p of problem (2.4) is defined as Eh,u,p (v h , qh ) := lhS (v h ) + lhD (qh ) − Ah ((I kh u, πhk p), (v h , qh ))

∀(v h , qh ) ∈ V kh × Qkh ,

(4.18)

where πhk is the L2 -orthogonal projector onto Qkh such that (πhk p)|T = πTk (p|T ). Afterwards, we have the following results: Lemma 4.6. Let u, p be the exact solution of problem (2.4) with the additional regularity ((u)|ΩS , p) ∈ ˇ T := I kT (KT ∇ˇ pT ). Then, H 2 (ΩS ) × H 1 (Ω). Furthermore, for all T ∈ ThD , denote pˇT := πTk+1 (p|T ) and u the consistency error Eh,u,p (v h , qh ) satisfies S D (v h ) + Eh,u,p (v h ) Eh,u,p (v h , qh ) = Eh,u,p

∀(v h , qh ) ∈ V kh × Qkh ,

(4.19)

S D where Eh,u,p (v h ), Eh,u,p (v h ) are defined as follows:    S (v h ) = ν(vF − vT , ∇(u − rTk+1 I kT u)nT F )F − νsST (I kT u, v T ), Eh,u,p T ∈ThS F ∈FT

−

 

(vF − vT , (p − πTk p)nT F )F +



√ α ν

(Pτ (u − πFk u), Pτ (vF ))F ,

D ) trace(KTΓF

  k ˇ T , v T ) − (DTk v T , πTk p − pˇT )T − aD (vF , (p − pˇT )nT F )F . T (I T u − u

T ∈ThS F ∈FT D Eh,u,p (v h ) = −



T ∈ThS

F ∈FhΓ



T ∈ThD F ∈FT

T ∈ThD

11

Proof. Using the equations of (2.1a), (2.2a) and the definition of lhS in (3.20): For all v h ∈ V kh , it holds  (ff , vT )T lhS (v h ) = T ∈ThS

=





 ν (∇u, ∇vT )T + (vF − vT , ∇unT F )F F ∈FT

T ∈ThS



+



− (∇ · vT , p)T −

(vF − vT , pnT F )F

F ∈FT

T ∈ThS −1

+ (K u, vT )ΩD + (∇p, vT )ΩD

  − ν(vF , ∇unT F )F − (vF , pnT F )F , T ∈ThS F ∈FT

where (K−1 u, vT )ΩD + (∇p, vT )ΩD = 0 will vanish due to (2.2a). What’s more, using the fact that ν∇unT F , pnT F are continuous at the interior faces, then for all v h ∈ V kh , we have

  ν(vF , ∇unT F )F − (vF , pnT F )F − T ∈ThS F ∈FT



=

F ∈(FhS \FhΓ )

=



(vF , [pnT F ])F +

F ∈(FhS \FhΓ )



(vF , ν∇unT F )F

F ∈(FhS ∩FhΓ )



−



(vF , ν[∇unT F ])F −



(vF , pD nS )F +

F ∈FhΓ

(vF , pnT F )F

F ∈(FhS ∩FhΓ )



F ∈FhΓ

√ α ν

D ) trace(KTΓF

(Pτ (u), Pτ (vF )),

(4.20)

where we have used the interface conditions (2.3), and the notation [·], borrowed from the discontinuous Galerkin literature (c.f., e.g., [16]), denotes the jump. From the commuting property (3.7), and the L2 -projection, the following equality holds    lhD (qh ) = (∇ · u, qT )T = (πTk ∇ · u, qT )T = (DTk I kT u, qT )T , T ∈Th

T ∈Th

T ∈Th

which means that lhD (qh ) − bh (I kh u, qh ) = 0. Next, let us compute ah (I kh u, v h ) − bh (v h , πhk p). From (3.20), we have ah (I kh u, v h ) − bh (v h , πhk p)

  k k aST (I kT u, v T ) − bT (v T , πTk p) + aD = (I u, v ) − b (v , π p) T T T T T T T ∈ThS

+



F ∈FhΓ

T ∈ThD

√

α ν  (Pτ (πFk u), Pτ (vF ))F . D ) trace(KTΓF

On the one hand, for all T ∈ ThS , using the fact that rTk+1 I kT is the L2 -elliptic projector ([49, Equation (4.6)]), it holds  (∇rTk+1 I kT u, ∇rTk+1 v T )T = (∇vT , ∇rTk+1 I kT u)T + (vF − vT , ∇rTk+1 I kT unT F )F 

= (∇u, ∇vT )T +

F ∈FT

F ∈FT

(vF − vT , ∇rTk+1 I kT unT F )F ,

so we obtain aST (I kT u, v T ) = ν(∇u, ∇vT )T +

 F ∈FT

ν(vF − vT , ∇rTk+1 I kT unT F )F + νsST (I kT u, v T )T .

12

Furthermore, −bT (v T , πhk p) = −(∇ · vT , πTk p)T −



(vF − vT , (πTk p)nT F )F

F ∈FT



= −(∇ · vT , p)T −

(vF − vT , (πTk p)nT F )F .

F ∈FT

On the other hand, for all T ∈ ThD , from the polynomial preservation property (3.9), the bilinear forms − bT (v T , πTk p) can be expressed as follows

k aD T (I T u, v T )

k k aD T (I T u, v T ) − bT (v T , πT p) k ˇ T , v T ) − (DTk v T , πTk p − pˇT )T + aD =aD uT , v T ) − (DTk v T , pˇT )T T (I T u − u T (ˇ    k ˇ T , v T ) − (DTk v T , πTk p − pˇT )T + =aD (vF , (p − pˇT )nT F )F − (vF , pD nD )F , T (I T u − u T ∈ThD F ∈FT

F ∈FhΓ

where we have used (3.17), the homogeneous Dirichlet boundary condition of v h and the fact that p is single valued at interior faces. We complete the proof of the lemma by collecting the above equalities. Lemma 4.7. Let (u, p) and (uh , ph ) be the solutions to (2.4) and (3.20), respectively. Moreover, assume that (uS )|Γ ∈ H k+1 (Γ), uS ∈ H k+2 (ΩS ) and pD ∈ H k+2 (ΩD ). Then, for all v h ∈ V kh , qh ∈ Qkh , the following estimates hold S |Eh,u,p (v h )|  c5 hk+1 (ν D |Eh,u,p (v h )|

 c6 h

k+1

1/2

uk+2,ΩS + ν

ΓF

defined in Lemma 4.6.

pk+1,ΩS + ν

1/4

uk+1,Γ )  v h ,

pk+2,ΩD  v h  .

1 α ) /2 }}, trace(KT D )

where c5 = max{1, maxF ∈FhΓ {( 

− 1/2

S D c6 = max{1, maxT ∈ThD ρKT }, and Eh,u,p (v h ), Eh,u,p (v h ) are

Proof. Using the definition (4.1) of norm  · , Lemma 3.5, approximation property (3.8), and the CauchySchwarz inequality, we have   ν(vF − vT , ∇(u − rTk+1 I kT u)nT F )F T ∈ThS F ∈FT



 

 ν

T ∈ThS F ∈FT

ν

1/2

2 h−1 T vF − vT F

1/2

 

ν

T ∈ThS F ∈FT

hT ∇(u − rTk+1 I kT u)nT F 2F

1/2

hk+1 uk+2,ΩS  v h  .

Using the stability property (3.16), it yields    1 1 1 ν sST (I kT u, v T )  |ν sST (I kT u, I kT u)| /2 |ν sST (v T , v T )| /2  ν /2 hk+1 uk+2,ΩS  v h  . T ∈ThS

T ∈ThS

T ∈ThS

According to the trace inequality,   (vF − vT , (p − πTk p)nT F )F T ∈ThS F ∈FT



 ν

 

T ∈ThS F ∈FT

ν

− 1/2

2 h−1 T vF − vT F

1/2

ν −1

  T ∈ThS F ∈FT

hk+1 pk+1,ΩS  v h  .

13

hT (p − πTk p)nT F 2F

1/2

For the interface term, we have √  α ν  (Pτ (u − πFk u), Pτ (vF ))F D trace(K ) Γ TΓF F ∈Fh √ √ 

1/2 

1/2 α ν α ν k 2   u − πF uF Pτ (vF )2F  D ) D ) trace(KTΓF trace(KTΓF F ∈F Γ F ∈F Γ h

 max (  F ∈FhΓ

h

α

)

1/2

ν

1/4

D ) trace(KTΓF

By setting c5 = max{1, maxF ∈FhΓ {( 

hk+1 uk+1,Γ  v h  .

1 α ) /2 }} trace(KT D )

and combining the above inequalities, we get the estima-

ΓF

S tion of Eh,u,p (v h ). D As for the Eh,u,p (v h ) term, the following estimations have been shown in [28, Section 4.3], such that for all D T ∈ Th , there holds k ˇ T , v T )  ρKT hk+1 aD T (I T u − u T pk+2,T v T D,T ,

(DTk v T , πTk p − pˇT )T  hk+1 T pk+2,T v T D,T , ˇ T and pˇT have been defined in Lemma 4.6. Subsequently, using the trace inequality, we have where u

1/2 

1/2   2 2 (vF , (p − pˇT )nT F )F ≤ h−1 p − p ˇ  h v · n   hk+1 T T F T F F F T T pk+2,T v T D,T . F ∈FT

F ∈FT

F ∈FT

and setting c6 = max{1, maxT ∈ThD ρKT } to obtain the estimation of Summing over the element T ∈ this completes the proof of the lemma. ThD

D (v h ), Eh,u,p

4.3. Error estimates With the consistency error and related estimations established in Lemma 4.6 and Lemma 4.7, we now present the energy-error estimate of discrete scheme (3.20). Theorem 4.8. Under the assumptions of Lemma 4.7, let (uh , ph ) ∈ V kh × Qkh solve the discrete equation (3.20). Then the following estimate holds

1 −1 1 uh − I kh u  +ph − πhk p  γ −1 c7 hk+1 ν /2 uk+2,ΩS + ν /2 pk+1,ΩS + pk+2,ΩD + ν /4 uk+1,Γ , (4.21) where γ is defined in Lemma 4.5, c7 = max{c5 , c6 } with c5 and c6 are defined in Lemma 4.7. Proof. By invoking [61, Theorem 10] and Lemma 4.5, we have (uh − I kh u, ph − πhk p)Zh ≤ γ −1 Eh,u,p (v h , qh )Zh∗ , where ·Zh∗ denotes the norm dual to ·Zh . Using Lemma 4.6 and Lemma 4.7, we obtain Eh,u,p (v h , qh )Zh∗ =

|Eh,u,p (v h , qh )| E h,u,p (v h , qh )Zh (v h ,qh )∈Zh \{0} sup

|Eh,u,p (v h , qh )| v h  (v h ,qh )∈Zh \{0}

1 −1 1  max{c5 , c6 }hk+1 ν /2 uk+2,ΩS + ν /2 pk+1,ΩS + pk+2,ΩD + ν /4 uk+1,Γ . (4.22) ≤

sup

Collecting the above estimates and setting c7 = max{c5 , c6 } yields the assertion. 5. Numerical experiments In this section, we present some two-dimensional numerical experiments to confirm our theoretical results. The polynomial degrees range from k = 0 to k = 5 are employed in different examples to get numerical discretization. The convergence of errors uh − I kh u and ph − πhk p are presented in the first three examples. Furthermore, the more realistic simulation of complex regions are illustrated in the last two examples. 14

5.1. Example 5.1 Let the computational domain be Ω = (0, 1) × (0, 2), where ΩS = (0, 1) × (1, 2), ΩD = (0, 1) × (0, 1), and √ α ν √ Γ = (0, 1) × {1}. Choose = 1, ν = 1 and K = I, where I is the unit matrix, the exact solution is given trace(K)

by ⎧

 πx 1 πx ⎪ 2 πy ⎪ ⎪uS = − cos ( ) sin( ), cos( )(sin(πy) + πy) , ⎪ 2 2 4 2 ⎪ ⎪ ⎪ πx π ⎪ 2 πy ⎪ ⎨pS = − cos( )(y − 2 cos ( )), 4 2 2 π2 ⎪ πx π πx  ⎪ ⎪ u sin( )y, cos( ) , = − ⎪ D ⎪ 8 2 4 2 ⎪ ⎪ ⎪ ⎪ ⎩pD = − π cos( πx )y. 4 2

(5.1)

The boundary conditions and the source terms follow from the exact solution. In this numerical example we test different type meshes, as shown in Figure 3, i.e., uniform triangulation mesh Th1 , nonconforming rectangle mesh Th2 , quadrilateral mesh Th3 by perturbing the interior nodes of uniform rectangle mesh with a parameter 0.25 (and see the package [62] for details), polygonal mesh Th4 generated by the dual of the triangle mesh Th1 , distorted polygonal mesh Th5 , hybrid distorted and rectilinear mesh Th6 , centroid Voronoi Tessellation (CVT) mesh Th7 generated by PolyMesher package [63], and the non-convex mesh Th8 . Here we take the polynomial degree k = 2 for all variables. The errors and the corresponding estimated convergence rates(ECR) are presented in Table 1, where Dof refers to the number of (cell-based and face-based) degree of freedoms (DOFs) of the discrete velocity uh . The conclusions of Theorem 4.8 are confirmed by the numerical results Table 1, with the energy-error of the velocity and the L2 -error of the pressure converging as order O(h3 ) as expected.

(a) Th1

(b) Th2

(c) Th3

(d) Th4

(e) Th5

(f) Th6

(g) Th7

(h) Th8

Figure 3: Illustrations of meshes.

15

Table 1: The errors for a series of the meshes Th1 −Th8 for Example 5.1. (a) Th1 −Th4

Th Th1

Th2

Th3

Th4

Dof uh − I kh u 1240 1.2560e-02 4816 1.5841e-03 18976 1.9839e-04 75328 2.4806e-05 1802 1.3598e-02 6884 1.8294e-03 26888 2.3362e-04 106256 2.9470e-05 728 2.5103e-02 2768 3.1266e-03 10784 3.7615e-04 42560 4.7607e-05 1367 1.1825e-02 4311 1.7542e-03 14999 2.3608e-04 55575 3.0540e-05

ECR – 2.99 3.00 3.00 – 3.00 3.03 3.02 – 3.01 3.06 2.98 – 3.14 3.15 3.10

(b) Th5 −Th8

ph − πhk p 2.1987e-03 2.9324e-04 3.7583e-05 4.5474e-06 1.0069e-03 1.1888e-04 1.1040e-05 9.7044e-07 1.6549e-03 1.9705e-04 2.3431e-05 2.8863e-06 6.8164e-04 7.7121e-05 9.0983e-06 1.1051e-06

Th

ECR – 2.91 2.96 3.05 – 3.20 3.50 3.55 – 3.07 3.07 3.02 – 3.58 3.36 3.20

Th5

Th6

Th7

Th8

Dof uh − I kh u 1367 1.3411e-02 4311 2.0514e-03 14999 2.7566e-04 55575 3.5587e-05 728 4.4967e-02 2768 5.0459e-03 10784 6.1620e-04 42560 7.6316e-05 809 1.7442e-02 3009 2.3459e-03 14509 2.0363e-04 51209 3.1145e-05 1088 1.2855e-02 4064 1.6205e-03 15680 2.0137e-04 61568 2.5046e-05

ECR – 3.09 3.15 3.10 – 3.16 3.03 3.01 – 2.74 2.93 2.89 – 3.03 3.03 3.02

ph − πhk p 9.7912e-04 1.3150e-04 1.5823e-05 1.9458e-06 3.3150e-03 3.4028e-04 3.9491e-05 3.6438e-06 1.0442e-03 1.1773e-04 7.9397e-06 1.2491e-06 1.5398e-03 2.0089e-04 2.5417e-05 3.1956e-06

ECR – 3.30 3.33 3.18 – 3.28 3.11 3.44 – 2.98 3.23 2.85 – 2.98 3.00 3.00

5.2. Example 5.2 By taking polynomial degrees k = 0, ..., 5, we investigate on the numerical errors and convergence rates of the proposed method. Considering the coupled Stokes-Darcy problem (2.1)−(2.3) in Ω = (0, 1) × (0, 2) with ⎧

 2 2 2 3 ⎪ u x(y − 1) = x (y − 1) + y, − + 2 − π sin(πx) , ⎪ S ⎪ ⎪ 3 ⎪ ⎪ ⎪ πy ⎨ pS = (2 − π sin(πx)) sin( ), 2 (5.2)

 ⎪ ⎪ 2 ⎪ uD = − π cos(πx)(1 − y − cos(πy)), −(2 − π sin(πx)(−1 + π sin(πy))) , ⎪ ⎪ ⎪ ⎪ ⎩ pD = (2 − π sin(πx))(1 − y − cos(πy)). where ΩS = (0, 1) × (1, 2), ΩD = (0, 1) × (0, 1), and Γ = (0, 1) × {1}. All the parameters remain unchanged as set in Example 5.1. The boundary conditions and the source terms follow from the exact solution. For various polynomial degrees k, the different errors for the HHO method on polygon meshes Th7 are displayed in Figure 4. This example shows the strong inclusiveness of the HHO method for polynomial degrees. Particularly, we show that the HHO scheme (3.20a)−(3.20a) achieves the order of convergence in accordance with Theorems 4.8. 5.3. Example 5.3 In this case we focus on the polynomial degree k = 0. We show that even with the lowest degree the proposed method still works well with large hydraulic conductivity tensor K. Let the computational domain be Ω = (0, 1) × (0, 2), where ΩS = (0, 1) × (1, 2), ΩD = (0, 1) × (0, 1) and Γ = (0, 1) × {1}. The exact solution is given by ⎧

 ⎪ ⎪ uS = y 2 − 2y + 1, x2 − x , ⎪ ⎪ ⎪ ⎪ ⎪ 1 ⎪ ⎪ ⎨pS = 2ν(x + y − 1) + , 3K (5.3)

 ⎪ 2 ⎪ u = − 2νK − ((1 − 2x)(y − 1)), −(x(1 − x) + y − 2y + 1) , ⎪ D ⎪ ⎪ ⎪ ⎪ ⎪ 1 1 ⎪ ⎩pD = 2νx + (x(1 − x)(y − 1) + y 3 − y 2 + y). K 3 Set K = Kj I, j = 1, ..., 4, where K1 = 10−1 , K2 = 10−3 , K3 = 10−6 , K4 = 10−12 . All the remaining parameters still retain unchanged as set in Example 5.1. Considering that the exact solutions of p is inversely proportional 16

10 0

k=0 k=1 k=2 k=3 k=4 k=5

10 0

10 -2

k=0 k=1 k=2 k=3 k=4 k=5

10 -2

10 -4

10 -4 10 -6

10 -6 10 -8 6 6

10 -8

10 -10 1

1

1

0.0313

0.0625

0.125

1

0.25

0.0313

(a) uh − I kh u vs. h

0.0625

0.125

0.25

k p vs. h (b) ph − πh

Figure 4: Errors for the numerical Example 5.2. p −π k p

to the coefficient K, we give the relative errors hph in Table 2. As can be seen from Table 2, the HHO method still has good convergence order and maintains good error precision with only a few DOFs, which has a great attraction in the practical applications. Table 2: The errors with different Kj for Example 5.3.

K

K1

K2

K3

K4

Dof 168 624 2400 9408 168 624 2400 9408 168 624 2400 9408 168 624 2400 9408

uh − I kh u 7.8394e-01 4.1837e-01 2.1493e-01 1.0842e-01 6.2996e-01 3.8684e-01 2.0925e-01 1.0747e-01 6.2650e-01 3.8626e-01 2.0927e-01 1.0760e-01 6.2649e-01 3.8626e-01 2.0927e-01 1.0760e-01

ECR – 0.91 0.96 0.99 – 0.70 0.89 0.96 – 0.70 0.88 0.96 – 0.70 0.88 0.96

k ph −πh p p

1.7907e-01 6.3255e-02 1.8783e-02 5.2092e-03 1.0612e-01 2.9125e-02 7.5027e-03 1.9079e-03 1.0517e-01 2.8591e-02 7.2338e-03 1.7866e-03 1.0517e-01 2.8591e-02 7.2336e-03 1.7864e-03

ECR – 1.50 1.75 1.85 – 1.87 1.96 1.98 – 1.88 1.98 2.02 – 1.88 1.98 2.02

5.4. Example 5.4 We now consider a numerical example to verify the mass conservation for the proposed HHO method. This experiment has been conducted in [64], however, we have made some changes to simulate practical problems with curved interface and nonconforming meshes. Let the computational domain Ω be (0, 2) × (0, 2), and the 3 sin( π2 x) cos(πx) + 1. The ΩS is the free-flow domain with a quadratic inflow interface Γ is defined by y = − 20 profile uS = (0, x(x − 2)) on the top boundary and no-slip boundary conditions uS = 0 on the left and right boundaries. For the Darcy domain ΩD , Neumann boundary condition uD ·n = 0 is imposed on the left and right boundary and Dirichlet boundary condition pD = 0 is imposed on the bottom boundary. What’s more, in order to test different permeability variation coefficients, we divide the Darcy ΩD into five parts. The computational domain and mesh are displayed in Figure 5(a) and Figure 5(b), respectively. As mentioned in Remark 2.1, the original curved interface has been replaced by short straight line segments, which coincide with FhΓ .

17

We choose √ α

√

ν trace(K)

= 1, ν = 1, and K = Kj I (j = 1, · · · , 6), where K1 − K6 take the values as follows:

• K1 = 1, for all ΩDi , i = 1, · · · , 5; • K2 = 10−2 , for all ΩDi , i = 1, · · · , 5; • K3 = 10−4 , for all ΩDi , i = 1, · · · , 5; • K4 = 10−6 , for all ΩDi , i = 1, · · · , 5; • K5 – for ΩDi , i = 1, 3, 5, taking the value of 1; – for ΩD2 , taking the value of 10−4 ; – for ΩD4 , taking the value of 10−6 ; • K6 – for ΩD1 , taking the value of 1; – for ΩDi , i = 2, 3, 4, 5, generated randomly on the interval [10−4 , 10−1 ] such that it takes a different constant value on each mesh element. Setting fS = 0 in ΩS , fD = 0 in ΩD and taking k = 1 for all variables. We show the velocity streamlines for different Kj in Figure 6, where the warmer color indicates higher velocity of the flow. Figure 7 displays the distribution of pressure for different Kj . The boundary condition on left and right sides of the whole domain imply that there is no mass lost there. So we use the difference mass of integrals of normal flux between top and bottom sides, namely inflow and outflow sides, to illustrate the mass conservation of HHO method. It is easy to compute the incoming flux integral equal to 4/3. In Table 3, we show the mass error for different Kj , which well illustrate the mass conservation of the method.

(a) Computational domain

(b) Computational mesh

Figure 5: Computational domain and mesh for Example 5.4.

Table 3: Mass difference between the inflow and outflow sides with different Kj for Example 5.4.

Dof K1 96869 4.7162e-13 149740 4.0810e-13 290440 1.7163e-12

K2 3.8842e-13 5.1859e-13 8.7658e-13

K3 2.1773e-13 3.4923e-12 4.6634e-12

18

K4 5.2830e-12 4.2719e-13 2.8839e-13

K5 4.7348e-13 3.8391e-12 5.3081e-12

K6 3.1141e-12 5.2084e-13 3.8772e-12

Figure 6: The velocity streamlines with different Kj for Example 5.4.

5.5. Example 5.5 In the last example, we consider the flow of a single phase fluid in the free fluid region and fractured porous medium region. The simulation of coupled problem is of great interest for a large number of applications, ranging from energy production to water resources management: oil fields exploitation, geothermal energy, carbon dioxide storage, and groundwater contamination [65, 66, 67]. Those papers modeled the fracture as the interface between the surrounding porous media with some interface conditions, and the dimension of the fracture is reduced to be 1D. While in our example, we treat the fractures as the full 2D medium. Therefore, we do not establish a new model, but set different permeability coefficients to distinguish fractures from surrounding medium. Let the computational domain be Ω = (0, 1) × (1, 2), where ΩS = (0, 1) × (1, 2), ΩD = (0, 1) × (0, 1) and Γ = (0, 1) × {1}. The ΩS is the free-flow domain with a quadratic inflow profile uS = (0, 4x(x − 1)) on the top boundary and no-slip boundary conditions uS = 0 on the left and right boundaries. For the porous medium region ΩD , Neumann boundary condition uD · n = 0 is imposed on the left and right boundary and Dirichlet boundary condition pD = 0 is imposed on the bottom boundary. We consider the following test cases: i) a single fracture, illustrated in Figure 8(a); ii) the intersect fractures, illustrated in Figure 8(d). While, the corresponding computational meshes are given in Figure 8(b) and 8(e), respectively; The magnified meshes for Darcy region are displayed in Figure 8(c) and 8(f), respectively. We assume that the fractures elements have a small width h/2.5, the fractures are given by Γ1 = {x ∈ [0.3, 0.75] : y = cos(3x)3 + 0.4}, 1}, Γ4 = {x ∈ [0.5, 0.8], y = 0.1x + 0.2} Γ2 = {x ∈ [0.3, 0.6] : y = 0.8x + 0.15}, Γ3 = {x ∈ [0.4, 0.55], y = 3x − √ = 1, ν = 1, fS = 0 and fD = 0. and Γ5 = {x ∈ [0.7, 0.85], y = −1.5x + 1.95}. Moreover, we choose √ α ν trace(K)

We assume that the fractures are filled with debris and the flow in the fractures respects Darcy’s law. We distinguish two types of fractures: fractures which have a permeability higher than that in the surrounding medium and those in which the permeability is lower than that in the surrounding medium. We consider the permeability K = KP I and K = KF I in medium and fractures, respectively. We have considered the following cases: • Consider the case of one single fracture (Figure 8(a)):

19

Figure 7: The pressure with different Kj for Example 5.4.

– KP = KF . In this case the pressure is expected to be almost constant across the fractures, because the permeability in the fractures is equal to that of the surrounding medium and the fluid is not affected by the presence of the fractures. Set KP = KF = 1, we show the velocity streamlines in Figure 9(a), where the warmer color indicates higher velocity of the flow. Moreover, we show the pressure in Darcy region in Figure 10(a). Finally, we test the effect of KP and KF on the solution, Figure 9(b) and Figure 10(b) show the simulation results for KP = KF = 10−4 . We can observe that when KP and KF becomes smaller, the flow velocity in porous media is significantly reduced. – KP > KF . In this case the fracture permeability is smaller than that of the surrounding medium, therefore we expect a pressure jump across the fracture. Set KP = 1 and KF = 10−4 , we show the velocity streamlines in Figure 9(c) for a single fracture. The corresponding pressure in Darcy region is shown in Figure 10(c). The fractures are much less permeable than the surrounding medium and the strong pressure jumps across fractures can be observed. So fractures with low permeability which act as barrier to the flow and the fluid in the porous medium are hindered from the fractures with low permeability values. – KP < KF . In this case the fracture pressure is expected to be almost continuous across the fracture and the flow is expected to be directed towards the fracture. This is because that the fracture is very permeable and the fluid is “attracted” by the fracture. Set KP = 10−4 , KF = 1, we show the velocity streamlines in Figure 9(d) for a single fracture. The distribution of Darcy pressure is displayed in Figure 10(d). The fluid has a tendency to flow into the fractures and then along the fractures. Indeed, pressure is practically continuous across the fractures and the maximum value of pressure is slightly lower than the non-fractured case. • Consider the case of the intersect fractures (Figure 8(d)): – Firstly, we take KP = 10−3 and KF = 10−10 , 10−3 , 10−3 , 10−10 , 1 for Γ1 − Γ5 , respectively. From the velocity streamlines in Figure 9(e), we can see that Γ2 and Γ3 , which have the same permeability coefficient as the surrounding porous media, have no effect on the fluid. While the fluid in the porous

20

medium avoid the fractures Γ1 and Γ4 , but seem to be “attracted” by the fractures Γ5 . From the corresponding Darcy pressure in Figure 10(e), we observe that the big pressure jump happens across the fractures Γ1 and Γ4 . – Another case is that we take KP = 10−3 and KF = 1, 1, 1, 10−3 , 10−10 for Γ1 − Γ5 . The velocity streamlines and Darcy pressure are given in Figure 9(f) and Figure 10(f), respectively. We observe the phenomenon similar to the above case.

(a) One fracture

(b) Computational domain mesh

(c) Darcy domain mesh

(d) Intersecting fractures

(e) Computational domain mesh

(f) Darcy domain mesh

Figure 8: Fractured domain, computational domain mesh and Darcy domain mesh for Example 5.5.

6. Conclusions In this paper, the hybrid high-order (HHO) method is derived to solve the coupled Stokes-Darcy problem with Beavers-Joseph-Saffman interface conditions on general meshes. We analyze the existence and uniqueness of discrete solutions to this coupled problem. The energy- and L2 -error estimates of order (k + 1) for velocity and pressure are obtained. In fact, according to the theories of [47, 28], the L2 -error estimates of order (k + 2) for uS and pD should be valid. However, in the current framework of error estimates, the L2 -error estimate of Stokes equation involves the energy-error estimation of Darcy equation, but we cannot obtain a higher order error estimation from this term, which makes us unable to obtain the optimal error estimation. Moreover, from the duality theory, the duality problem of the Darcy equation contains interface term. How to deal with this interface term to achieve the optimal error estimation needs further discussion. Therefore, lots of work remains to be done to obtain the L2 -error estimates, which will be analyzed in our subsequent articles. Additionally, some numerical experiments are provided to illustrate that the proposed method supports various kinds of meshes in two regions, achieves the theoretical convergence order for different polynomial degrees, has the properties of mass conservation and robustness with respect to hydraulic conductivity.

21

Figure 9: The velocity and velocity streamlines with different Kj for Example 5.5.

References [1] P. Chidyagwai, B. Rivi`ere, Numerical modelling of coupled surface and subsurface flow systems, Adv. Water Resour. 33 (1) (2010) 92–105. [2] W. Layton, H. Tran, C. Trenchea, Analysis of long time stability and errors of two partitioned methods for uncoupling evolutionary groundwater-surface water flows, SIAM J. Numer. Anal. 51 (1) (2013) 248–272. [3] S. Stoter, P. M¨ uller, L. Cicalese, M. Tuveri, D. Schillinger, T. Hughes, A diffuse interface method for the Navier-Stokes/Darcy equations: Perfusion profile for a patient-specific human liver based on mri scans, Comput. Methods Appl. Mech. Engrg. 321 (2017) 70–102. [4] N. S. Hanspal, A. N. Waghode, V. Nassehi, R. J. Wakeman, Numerical analysis of coupled Stokes/Darcy flows in industrial filtrations, Transp. Porous Media 64 (1) (2006) 73. [5] N. S. Hanspal, A. N. Waghode, V. Nassehi, R. J. Wakeman, Development of a predictive mathematical model for coupled Stokes/Darcy flows in cross-flow membrane filtration, Chem. Eng. J. 149 (1) (2009) 132–142. [6] W. J. Layton, F. Schieweck, I. Yotov, Coupling fluid flow with porous media flow, SIAM J. Numer. Anal.

22

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10: The pressure of Darcy domain with different KP and KF for Example 5.5.

40 (6) (2002) 2195–2218. [7] Y. Cao, M. Gunzburger, X. Hu, F. Hua, X. Wang, W. Zhao, Finite element approximations for Stokes-Darcy flow with Beavers-Joseph interface conditions, SIAM J. Numer. Anal. 47 (6) (2010) 4239–4256. [8] A. M´arquez, S. Meddahi, F.-J. Sayas, Strong coupling of finite element methods for the Stokes-Darcy problem, IMA J. Numer. Anal. 35 (2015) 969–988. [9] Y. Cao, M. Gunzburger, X. He, X. Wang, Robin-Robin domain decomposition methods for the steady-state Stokes-Darcy system with the Beavers-Joseph interface condition, Numer. Math. 117 (4) (2011) 601–629. [10] M. Discacciati, A. Quarteroni, A. Valli, Robin-Robin domain decomposition methods for the Stokes-Darcy coupling, SIAM J. Numer. Anal. 45 (3) (2007) 1246–1268. [11] I. Babuˇska, G. N. Gatica, A residual-based a posteriori error estimator for the Stokes-Darcy coupled problem, SIAM J. Numer. Anal. 48 (2) (2010) 498–523. [12] G. N. Gatica, R. Oyarz´ ua, F.-J. Sayas, A residual-based a posteriori error estimator for a fully-mixed formulation of the Stokes-Darcy coupled problem, Comput. Methods Appl. Mech. Engrg. 200 (21) (2011) 1877–1891. [13] M. Mu, J. Xu, A two-grid method of a mixed Stokes-Darcy model for coupling fluid flow with porous media flow, SIAM J. Numer. Anal. 45 (5) (2007) 1801–1813. [14] L. Zuo, Y. Hou, A two-grid decoupling method for the mixed Stokes-Darcy model, J. Comput. Appl. Math. 275 (2015) 139–147. [15] B. Rivi`ere, I. Yotov, Locally conservative coupling of Stokes and Darcy flows, SIAM J. Numer. Anal. 42 (5) (2005) 1959–1977. [16] B. Rivi`ere, Analysis of a discontinuous finite element method for the coupled Stokes and Darcy problems, J. Sci. Comput. 22 (1) (2005) 479–500. [17] G. Wang, Y. He, R. Li, Discontinuous finite volume methods for the stationary Stokes-Darcy problem, Int. J. Numer. Meth. Eng. 107 (5) (2016) 395–418. [18] R. Li, Y. Gao, J. Li, Z. Chen, Discontinuous finite volume element method for a coupled non-stationary Stokes-Darcy problem, J. Sci. Comput. 74 (2) (2018) 693–727. 23

[19] J. Galvis, M. Sarkis, Non-matching mortar discretization analysis for the coupling Stokes-Darcy equations, Electron. Trans. Numer. Anal. (2007) 350–384. [20] V. Ervin, E. Jenkins, S. Sun, Coupling nonlinear Stokes and Darcy flow using mortar finite elements, Appl. Numer. Math. 61 (11) (2011) 1198–1222. [21] P. Hessari, Pseudospectral least squares method for Stokes-Darcy equations, SIAM J. Numer. Anal. 53 (3) (2015) 1195–1213. [22] S. M¨ unzenmaier, G. Starke, First-order system least squares for coupled Stokes-Darcy flow, SIAM J. Numer. Anal. 49 (1) (2011) 387–404. [23] Y. Boubendir, S. Tlupova, Stokes-Darcy boundary integral solutions using preconditioners, J. Comput. Phys. 228 (23) (2009) 8627–8641. [24] S. Tlupova, R. Cortez, Boundary integral solutions of coupled Stokes and Darcy flows, J. Comput. Phys. 228 (1) (2009) 158–179. [25] W. Wang, C. Xu, Spectral methods based on new formulations for coupled Stokes and Darcy equations, J. Comput. Phys. 257 (2014) 126–142. [26] D. A. Di Pietro, A. Ern, S. Lemaire, An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators, Comput. Meth. Appl. Math. 14 (4) (2014) 461–472. [27] D. A. Di Pietro, A. Ern, A hybrid high-order locking-free method for linear elasticity on general meshes, Comput. Meth. Appl. Mech. Engrg. 283 (2015) 1–21. [28] D. A. Di Pietro, A. Ern, Arbitrary-order mixed methods for heterogeneous anisotropic diffusion on general meshes, IMA J. Numer. Anal. 37 (1) (2017) 40–63. [29] L. Beir˜ao da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, A. Russo, Basic principles of virtual element methods, Math. Mod. Meth. Appl. Sci. 23 (2013) 199–214. [30] L. Beir˜ ao da Veiga, F. Brezzi, L. D. Marini, A. Russo, Virtual element method for general second-order elliptic problems on polygonal meshes, Math. Mod. Meth. Appl. S. 26 (04) (2016) 729–750. [31] F. Bassi, L. Botti, A. Colombo, D. A. Di Pietro, P. Tesini, On the flexibility of agglomeration based physical space discontinuous Galerkin discretizations, J. Comput. Phys. 231 (1) (2012) 45–65. [32] A. Cangiani, Z. Dong, E. H. Georgoulis, P. Houston, hp-Version discontinuous Galerkin methods on polygonal and polyhedral meshes, 1st Edition, SpringerBriefs in Mathematics, Springer International Publishing, 2017. [33] B. Cockburn, J. Gopalakrishnan, R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2) (2009) 1319–1365. [34] B. Cockburn, Static condensation, hybridization, and the devising of the HDG methods, Springer International Publishing, Cham, 2016, pp. 129–177. [35] J. Wang, X. Ye, A weak Galerkin finite element method for second-order elliptic problems, J. Comput. Appl. Math. 241 (2013) 103–115. [36] L. Mu, J. Wang, X. Ye, Weak Galerkin finite element methods on polytopal meshes, Int. J. Numer. Anal. Model. 12 (1) (2015) 31–53. [37] J. Droniou, R. Eymard, A mixed finite volume scheme for anisotropic diffusion problems on any grid, Numer. Math. 105 (1) (2006) 35–71. [38] R. Eymard, T. Gallou¨et, R. Herbin, Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes sushi: a scheme using stabilization and hybrid interfaces, IMA J. Numer. Anal. 30 (4) (2010) 1009–1043. [39] F. Brezzi, K. Lipnikov, M. Shashkov, Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes, SIAM J. Numer. Anal. 43 (5) (2005) 1872–1896. [40] K. Lipnikov, G. Manzini, M. Shashkov, Mimetic finite difference method, J. Comput. Phys. 257 (2014) 1163–1227. [41] J. Droniou, R. Eymard, T. Gallou¨et, R. Herbin, A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods, Math. Models Methods Appl. Sci. 20 (2) (2010) 265–295. [42] G. Wang, F. Wang, L. Chen, Y. He, A divergence free weak virtual element method for the Stokes-Darcy problem on general meshes, Comput. Methods Appl. Mech. Engrg. 344 (2019) 998–1020. [43] W. Chen, F. Wang, Y. Wang, Weak Galerkin method for the coupled Darcy-Stokes flow, IMA J. Numer.

24

[44] [45] [46] [47] [48] [49] [50]

[51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67]

Anal. 36 (2) (2016) 897–921. R. Li, Y. Gao, J. Li, Z. Chen, A weak Galerkin finite element method for a coupled Stokes-Darcy problem on general meshes, J. Comput. Appl. Math. 334 (2018) 111–127. K. Lipnikov, D. Vassilev, I. Yotov, Discontinuous Galerkin and mimetic finite difference methods for coupled Stokes-Darcy flows on polygonal and polyhedral grids, Numer. Math. 126 (2) (2014) 321–360. D. A. Di Pietro, J. Droniou, A. Ern, A discontinuous-skeletal method for advection-diffusion-reaction on general meshes, SIAM J. Numer. Anal. 53 (5) (2015) 2135–2157. D. A. Di Pietro, A. Ern, A. Linke, F. Schieweck, A discontinuous skeletal method for the viscosity-dependent Stokes problem, Comput. Meth. Appl. Mech. Engrg. 306 (2016) 175–195. L. Botti, D. A. Di Pietro, J. Droniou, A Hybrid high-order discretisation of the Brinkman problem robust in the Darcy and Stokes limits, Comput. Methods Appl. Mech. Engrg. 341 (2018) 278–310. D. A. Di Pietro, J. Droniou, A hybrid high-order method for Leray-Lions elliptic equations on general meshes, Math. Comp. 86 (307) (2017) 2159–2191. D. A. Di Pietro, J. Droniou, W s,p -approximation properties of elliptic projectors on polynomial spaces, with application to the error analysis of a hybrid high-order discretisation of Leray-Lions problems, Math. Models Methods Appl. Sci. 27 (5) (2017) 879–908. F. Chave, D. A. Di Pietro, F. Marche, F. Pigeonneau, A hybrid high-order method for the Cahn-Hilliard problem in mixed form, SIAM J. Numer. Anal. 54 (3) (2016) 1873–1898. B. Cockburn, D. A. Di Pietro, A. Ern, Bridging the hybrid high-order and hybridizable discontinuous Galerkin methods, ESAIM Math. Model Numer. Anal. 50 (3) (2016) 635–650. D. Boffi, D. A. Di Pietro, Unified formulation and analysis of mixed and primal discontinuous skeletal methods on polytopal meshes, ESAIM Math. Model Numer. Anal. 52 (1) (2018) 1–28. D. A. Di Pietro, J. Droniou, G. Manzini, Discontinuous skeletal gradient discretisation methods on polytopal meshes, J. Comput. Phys. 355 (2018) 397–425. D. A. Di Pietro, A. Ern, Hybrid high-order methods for variable-diffusion problems on general meshes, C. R. Acad. Sci Paris Ser. 353 (1) (2015) 31–34. G. S. Beavers, D. D. Joseph, Boundary conditions at a naturally permeable wall, J. Fluid Mech. 30 (1) (1967) 197–207. P. G. Saffman, On the boundary condition at the surface of a porous medium, Stud. Appl. Math. 50 (2) (1971) 93–101. D. A. Di Pietro, A. Ern, Mathematical aspects of discontinuous Galerkin methods, Springer Berlin Heidelberg, 2012. F. Chave, D. A. Di Pietro, L. Formaggia, A hybrid high-order method for Darcy flows in fractured porous media, SIAM J. Sci. Comput. 40 (2) (2018) A1063–A1094. F. Brezzi, M. Fortin, Mixed and hybrid finite element methods, Vol. 15 of Springer Series in Computational Mathematics, Springer-Verlag, New York, 1991. D. A. Di Pietro, J. Droniou, A third Strang lemma and an Aubin-Nitsche trick for schemes in fully discrete formulation, Calcolo 55 (3) (2018) Art. 40, 39. L. Chen, iFEM: an innovative finite element methods package in MATLAB, University of California at Irvine 2009. C. Talischi, G. H. Paulino, A. Pereira, I. F. M. Menezes, Polymesher: a general-purpose mesh generator for polygonal elements written in matlab, Struct. Multidiscip. Optim. 45 (3) (2012) 309–328. G. Kanschat, B. Rivi`ere, A strongly conservative finite element method for the coupling of Stokes and Darcy flow, J. Comput. Phys. 229 (17) (2010) 5933–5943. V. Martin, J. Jaffr´e, J. E. Roberts, Modeling fractures and barriers as interfaces for flow in porous media, SIAM J. Sci. Comput. 26 (5) (2005) 1667–1691. P. Angot, F. Boyer, F. Hubert, Asymptotic and numerical modelling of flows in fractured porous media, ESAIM Math. Model Numer. Anal. 43 (2) (2009) 239–275. L. Formaggia, P. F. Antonietti, A. Scotti, M. Verani, N. Verzotti, Mimetic finite difference approximation of flows in fractured porous media, ESAIM Math. Model Numer. Anal. 50 (3) (2016) 809–832.

25

Conflict of Interest Statement The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.