A posteriori error analysis of a fully-mixed formulation for the Brinkman–Darcy problem

We develop the a posteriori error analysis for a mixed finite element method applied to the coupling of Brinkman and Darcy equations in 3D, modelling the interaction of viscous and non-viscous flow effects across a given interface. The system is formulated in terms of velocity and pressure within the Darcy subdomain, together with vorticity, velocity and pressure of the fluid in the Brinkman region, and a Lagrange multiplier enforcing pressure continuity across the interface. The solvability of a fully-mixed formulation along with a priori error bounds for a finite element method have been recently established in Álvarez et al. ( Comput Methods Appl Mech Eng 307:68–95, 2016). Here we derive a residual-based a posteriori error estimator for such a scheme, and prove its reliability exploiting a global inf-sup condition in combination with suitable Helmholtz decompositions, and interpolation properties of Clément and Raviart–Thomas operators. The estimator is also shown to be efficient, following a localisation strategy and appropriate inverse inequalities. We present numerical tests to confirm the features of the estimator and to illustrate the performance of the method in academic and application-oriented problems.


Introduction
We have recently introduced a mixed finite element method to numerically approximate the flow patterns of a viscous fluid within a highly permeable medium described by Brinkman equations, and its interaction with pure porous media flow under Darcy's law [1]. There, the system is formulated in terms of velocity and pressure of the nonviscous flow, together with vorticity, velocity and pressure of the Brinkman region. The tangential vorticity vanishes on the boundary of the Brinkman domain, whereas slip velocity conditions are assumed on the overall boundary. The corresponding mixed variational formulation leads to a Lagrange multiplier enforcing pressure continuity across the interface, while mass balance results from essential boundary conditions on each domain. As a consequence, a classical saddle-point operator equation is obtained, whose invertibility hinges on the well-known Babuška-Brezzi theory. A similar treatment is used to establish the solvability of the discrete problem associated to the Galerkin method. The needed continuous and discrete inf-sup conditions can be guaranteed thanks to the so-called T -coercivity argument (cf. [18,26] and the references therein), where one defines adequate injective operators delivering lower bounds of the corresponding suprema. As the stability of the Galerkin scheme requires that the curl of the discrete vorticity space is contained in the discrete Brinkman velocity space, we specify Raviart-Thomas and Nédélec spaces for the approximation of the global velocity and the Brinkman vorticity, respectively.
On the other hand, the derivation of adaptive schemes for transmission free flow-porous media problems has been extensively studied in recent years. In particular, we refer to [6,[9][10][11]14,16,27], which are focused on Stokes-Darcy and Navier-Stokes/Darcy couplings, and where the interface conditions are treated in different ways, from both mathematical and numerical perspective. For instance, in [6,14,16,27], Beavers-Joseph-Saffman-type conditions are considered on the interface, whereas in [10,11], similar transmission conditions to those employed in [1] are assumed. Also, an interesting feature of the proof of reliability in [6], which differs from the approaches in the other works, is the utilisation of intermediate inf-sup inequalities that are obtained along the proof of the global inf-sup condition. Differently from the above, and similarly as in [14,16,27], the efficiency estimates in [6] follow from usual arguments based on inverse inequalities and the localisation technique employing triangle-bubble and edge-bubble functions. In turn, the assumption of a smallness condition on the data is the distinctive feature of the approach in [16], where a reliable and efficient residual-based a posteriori error estimator for the three dimensional version of the augmented-mixed method introduced in [17], is derived. Furthermore, an a posteriori error estimator for a conforming and nonconforming vorticity-based finite element method of a Stokes-Darcy coupled problem was derived in [10,11], respectively, but the resulting estimate in [11] is not optimal. In addition, even though in [10,11] the model problem is addressed for both two and three spatial dimensions, the corresponding a posteriori error analysis is explicitly derived only for the 2D case.
According to the previous discussion, and as a natural continuation of the a priori error analysis developed in [1], our goal in the present paper is to provide a reliable and efficient residual-based a posteriori error estimator for the finite element method introduced and analyzed in that reference. In this way, we aim to improve the accuracy of the discrete scheme from [1] in different scenarios, including presence of singularities or high gradients of the solution. Indeed, in contrast with the methodology developed in [10,11], and following the approaches in [16,27], we highlight that the derivation of our error estimator is based on a global inf-sup condition in combination with suitable Helmholtz decompositions adapted from [16,31], and local approximation properties of Clément, Raviart-Thomas, and Clément-type Nédélec interpolators. Then, similarly as in [14,16,27], the associated efficiency estimates are consequence of suitable inverse inequalities and local bounds for tetrahedron-bubble and facet-bubble functions.
The remainder of the paper is structured in the following manner. General preliminary notation is presented in the last part of this section. The model problem and the mixed variational formulation are outlined in Sect. 2, where we also recall its unique solvability and the mixed finite element discretisation. The core of the present analysis is contained in Sect. 3, where we define the error estimator and provide a detailed derivation of its reliability and efficiency. Finally, Sect. 4 gives two numerical tests aimed to illustrate the features of the method and the proposed estimator.
Some recurrent notation to be employed throughout the paper includes the following. If S ⊆ R 3 is a domain or a Lipschitz surface, and r ∈ R, we set vectorial Sobolev spaces as H r (S) := [H r (S)] 3 , adopt the convention H 0 (S) ≡ L 2 (S), and denote the corresponding norms by · r,S (for both H r (S) and H r (S)). In general, given a generic Hilbert space H, we will employ H to denote its vectorial counterpart H 3 . We also recall the definition of the Hilbert spaces where, for any vector field v : In addition we will use the space endowed with the the usual norm of L 2 (S). In turn, for each integer k ≥ 0 we denote by P k (S) the space of polynomials in S of total degree ≤ k, and we set P k (S) = [P k (S)] 3 . Finally, the symbol 0 will stand for a generic null vector (including the null functional and operator), and C (indistinguishably c, with or without subscripts, bars, tildes or hats) will denote generic constants independent of the discretisation parameters.
2 Governing equations and a mixed variational formulation

The continuous model
We first let B and D be bounded and simply connected polyhedral Lipschitz domains in R 3 such that ∂ B ∩ ∂ D =: = ∅ and B ∩ D = ∅, and set := B ∪ ∪ D with boundary = ∂ split into B ⊆ ∂ B and D ⊆ ∂ D . Then, given source terms f D ∈ L 2 ( D ) and f B ∈ L 2 ( B ), we are interested in the Brinkman-Darcy coupled problem 1) which is formulated in terms of the Brinkman velocity u B , the Brinkman pressure p B , the Brinkman vorticity ω B , the Darcy velocity u D , and the Darcy pressure p D . Here n stands for the outward normal at B and D , ν > 0 is the kinematic viscosity of the fluid, and κ D , κ B > 0 are the absolute permeabilities of the Darcy and Brinkman subdomains, respectively.
The boundary conditions on the Brinkman and Darcy subdomains suggest the following spaces In addition, the pressure continuity across the interface allows us to define its trace via the auxiliary unknown λ := p D | = p B | ∈ H 1/2 ( ), where ·, · denotes the duality pairing of H −1/2 ( ) and H 1/2 ( ) with respect to the L 2 ( )-inner product. In turn, the continuity of normal velocities across is imposed in a weak manner as Then, a fully-mixed formulation for (2.1) reads as follows: The well-posedness of (2.2) has been established in [1] using the classical Babuška-Brezzi theory:

Discretisation using a finite element method
Let T B h and T D h be respective partitions of B and D by shape-regular tetrahedra T of diameter h T . We assume that these tetrahedrisations match on the interface so that is a regular family of triangulations of = B ∪ ∪ D , with meshsize h := max{h T : T ∈ T h }. We denote by h the triangulation on induced by T h , which is formed by triangles F of diameter h F , and set h := max{h F : F ∈ h }. Next we introduce the finite-dimensional spaces where ∈ {B, D}, and for any T ∈ T h we denote by RT 0 (T ) := P 0 (T ) ⊕ P 0 (T ) x the local Raviart-Thomas space of lowest order. In addition, we set where for any T ∈ T B h , ND 1 (T ) := P 0 (T ) ⊕ P 0 (T ) × x is the local edge space of Nédélec type The approximation of the interface unknown will occur on an independent triangulation h of , by elements F of maximum diameter h := max h F : F ∈ h , where we define the space (2.8) In this way the Galerkin scheme associated to (2.2) reads: Find

A residual-based a posteriori error estimator
In this section we derive a reliable and efficient a posteriori error estimator for the Galerkin scheme (2.9). Most of the present proofs make extensive use of estimates available in [1,3,5,6,8,15,21,23,24,27].

Preliminaries
Given a tetrahedron T ∈ T h , we let E(T ) and F(T ) be the sets of its edges and faces, respectively. In addition, we denote by E h and F h be the sets of all edges and faces of T h , respectively, so that F h is subdivided as follows: Also, for each face F ∈ F h ( ) we fix a unit normal n F to F, so that given v ∈ L 2 ( ) such that v| T ∈ C(T ) on each T ∈ T h , and given F ∈ F h ( ), we let v × n F be the corresponding jump of the tangential traces across F, that is v × n F := (v| T − v| T )| F × n F , where T and T are the elements of T h having F as a common face. In addition, for each edge E of a tetrahedron T ∈ T h , we fix a unit tangential vector t E along E. When no confusion arises, we simple write n instead of n F , and t instead t E .
We now recall from [13] the tangential curl operator curl s : H 1/2 ( ) → L(H −1/2 ( )), with L(H −1/2 ( )) denoting the tangential vector fields of order −1/2, which is defined by curl s (χ ) := ∇χ × n, for any sufficiently smooth function χ . This is a linear and continuous map (see [13,Propositions 3.4 and 3.6]) which will be required in the sequel. We will also make use of the Raviart-Thomas interpolator of lowest order (see [22]) h : H 1 ( ) → H h , ∈ {B, D}, which according to its characterisation given by the identity where P h is the L 2 ( )-orthogonal projector onto P 0 ( ). In addition, we recall the Clément operator onto the space of the continuous piecewise linear functions and let I h : H 1 ( ) → X h be its vectorial counterpart defined component-wise. These maps satisfy the following properties (see [12,20,22], respectively) where T F is a tetrahedron of T h containing a face F on its boundary.
Furthermore, following [21] we define the Clément-type Nédélec interpolator N h : and λ E is the standard basis function for the lowest order Nédelec element, which satisfies where δ E,E is the Kronecker delta. The approximation properties of N h are summarised in the following Lemma (see [ We will also require stable Helmholtz decompositions for H (div; ) with * ∈ {B, D}. A technical assumption is that lies on the "convex part" of , signifying that there exists a convex domain containing , whose boundary contains . More precisely, introducing the space

Lemma 4 Assume that there exists a convex domain
such that ⊆ and ⊆ ∂ . Then, given v ∈ H (div; ) there exist w ∈ H 2 ( ) and β ∈ H 1 ( ) such that v = ∇w + curl β in and w 2, where C is a positive constant independent of all the foregoing variables.
In turn, a decomposition for H 0 (curl; B ) is given as follows.
Proof See [31, Lemma 2.2 and §5] We end this section with an estimate (in terms of local quantities) for the H −1/2 ( ) norm of functions in a particular subspace of H −1/2 ( ) ∩ L 2 ( ). According to the definition of Q h (cf. 2.8), we introduce the following orthogonal-type space

Lemma 6
Assume that for each F ∈ h there exists F ∈ h such that F ⊆ F and h ≤ C 1 h , with a constant C 1 > 0 independent of h and h . Then, there exists C > 0 independent of the aforementioned meshsizes, such that

Defining the proposed estimator
the unique solution of (2.9), we define for each T ∈ T B h , the local a posteriori error indicator B,T as follows: and for each T ∈ T D h , we define the local a posteriori error indicator D,T as It is not difficult to see that each term defining 2 B,T and 2 D,T is residual. Hence a global residual error estimator for (2.9) can be defined as (3.7) The remainder of this section advocates to establish the existence of positive constants C eff and C rel , independent of the meshsizes and the continuous and discrete solutions, such that where h.o.t stands, eventually, for one or several terms of higher order. The upper and lower bounds in (3.8), are derived below in Sects. 3.3 and 3.4, respectively.

Preliminary estimates
We begin by recalling that the first inequality in the continuous dependence result for all (w, r) ∈ H × Q 0 . This allows to establish a first estimate for the total error as follows.

10)
and satisfies Proof Applying (3.9) to the error (w, r) := (u, p) − (u h , p h ) and using (3.10) we arrive at Then, noting that obviously and applying the supremum in (3.12), we find that Next, employing the second equation of (2.2) and the definition of b, we deduce that Finally, from (3.10) and the first equation of (2.9), we obtain (3.11), and the proof concludes.
The next step consists in deriving suitable upper bounds for the residual term u B h · n − u D h · n −1/2, and for E H . We begin with the following result. Lemma 7 There exists C 4 > 0, independent of the meshsizes, such that Proof Taking ξ h ∈ Q h and then p h = (0, 0, ξ h ) ∈ Q h,0 in the second equation of (2.9), we find that In this way, (3.13) follows from a direct component-wise application of (3.4) (cf. Lemma 6).
We now aim to estimate E H . To this end, we first rewrite the functional as follows

Upper bound for
Given v B ∈ H B (div; B ), we consider its Helmholtz decomposition established in Lemma 4. More precisely, we let w ∈ H 2 ( B ) and β ∈ H 1 B ( B ) be such that v B = ∇w + curl β in B , and . (3.17) Consequently, in what follows we derive suitable upper bounds for the module of the two expressions on the right hand side of (3.17), which are provided by the following two lemmas.

Lemma 8
There exists C > 0, independent of meshsizes, such that for each w ∈ H 2 ( B ) there holds Proof Using the definition of the functional E 1 (cf. 3.14), the identity (3.2), the fact that p B h | F ∈ P 0 (F) for each F ∈ F h ( ), and the characterisation of B h given in (3.1), we find that In turn, the fact that ∇w ∈ H 1 ( B ) guarantees that (∇w − B h (∇w)) · n ∈ L 2 ( ), and hence which, together with (3.19), gives In this way, employing the Cauchy-Schwarz inequality, and the approximation properties of B h given in Lemma 1, we deduce from the above expression that which yields (3.18) and completes the proof.

Lemma 9
There exists C > 0, independent of meshsizes, such that for each β ∈ H 1 ( B ) there holds Proof Given β ∈ H 1 ( B ), we deduce from (3.14) and the identity div{curl(β − I B h β)} = 0, that In turn, thanks to the identities given in [30, Chapter I, Eq. (2.17) and Theorem 2.11], we find that which gives Now, integrating by parts in the first term on the right hand side of the last equation, we obtain (3.21) Applying Cauchy-Schwarz inequality, Lemma 2, and the uniform boundedness of the number of tetrahedra of the macro-elements B (T ) and B (F), we deduce from (3.21) that which implies (3.20) and ends the proof.
The following Lemma concludes the upper bound for E 1 H B (div; B ) .

Lemma 10 Assume that there exists a convex domain B such that B ⊆ B and
Then, there exists C 1 > 0, independent of meshsizes, such that where 2 B,T := 2 1,T + 2 2,T , that is Proof It follows from (3.18), (3.20), and the stability of the Helmholtz decomposition (3.16).

Upper bounds for E 2 H 0 (curl; B ) and E 3 H D (div; D )
We first establish the upper bound for E 3 H D (div; D ) , which is basically a "mirror reflection" through of Lemma 10.

Lemma 11 Assume that there exists a convex domain D such that D ⊆ D and
Then, there exists C 3 > 0, independent of the meshsizes, such that Proof It proceeds exactly as in the proofs of Lemmas 8, 9, and 10, by replacing B , B , and H B (div; B ) by D , D , and H D (div; D ), respectively. We omit further details.
The upper bound for E 2 H 0 (curl; B ) is provided next. Indeed, the derivation of this bound hinges on the Helmholtz decomposition given in Lemma 5, integration by parts, and the approximation properties of the Clément operators I B h and N h established in Lemmas 2 and 3, respectively.

Lemma 12
There exists C 2 > 0, independent of the meshsizes, such that Proof Given z B ∈ H 0 (curl; B ), we know from Lemma 5 that there exist ϕ ∈ H 1 0 ( B ) and χ ∈ H 1 0 ( B ), such that and Next, employing the operators I B h and N h defined in Sect. 3.1, we introduce the following discrete Helmholtz decomposition which clearly belongs to H B h,0 . In this way, and recalling from (3.15) that E 2 (z B h ) = 0, it follows that from which, according to the definition of E 2 (cf. 3.14), we find that Then, integrating by parts on each T , and noting that div ω B h is zero on T (cf. 2.6, 2.7), we have In this way, applying Cauchy-Schwarz inequality, the approximation properties of I B h and N h given in Lemmas 2 and 3, respectively, the fact that the number of tetrahedra of the macro-elements B (T ) and B (F) is uniformly bounded, and the stability estimate (3.23), we get (3.22) and finish the proof.
We end this section by concluding that the reliability of , that is the upper bound in (3.8), is a straightforward consequence of Theorem 2 and Lemmas 7, 10-12.

Efficiency
We now devote our attention to the derivation of upper bounds depending on the actual errors associated to the local indicators on each subdomain. For clarity of the analysis we will restrict ourselves to piecewise polynomial forcing terms f B and f D , but we remark that if they are otherwise sufficiently smooth, the error committed from suitable polynomial approximation would produce additional higher order terms in (3.8), explaining the eventual appearance of h.o.t in that inequality.
First, and thanks to the incompressibility condition in B (respectively D ), one has that The remaining terms in 2 B,T and 2 D,T can be treated very much in the same way as done in [24,27,28], where the analysis is based on inverse inequalities found in [19], together with the localisation technique based on tetrahedron-bubble and facetbubble functions [34]. Such a theory requires further notation and preliminary results collected in what follows.
Given T ∈ T h and F ∈ F(T ), let ψ T and ψ F denote tetrahedron-bubble and face-bubble functions, respectively (see [33,Eqs. (1.4) and (1.6)]), which satisfy: In addition, there exists an extension operator L : C(F) → C(T ) that satisfies L( p) ∈ P k (T ) and L( p)| F = p ∀ p ∈ P k (F), for a given k ≥ 0 (see [32]). The vectorial counterpart of L will be denoted L. Moreover, the following properties hold (where a proof can be found in [32,Lemma 4.1]).

Lemma 13
Given k ∈ N ∪ {0}, there exist c 1 , c 2 , c 3 > 0, depending only on k and the shape regularity of the triangulations, such that for each T ∈ T h and F ∈ F(T ), there hold The following inverse estimate is also required (see a proof in [ Finally we give two technical lemmas before tackling the derivation of the required upper bounds.

Lemma 15
Let ζ h ∈ L 2 ( ) be an element-wise polynomial of degree k ≥ 0, and let ζ ∈ L 2 ( ) be such that curl(ζ ) = 0 in . Then, there exist c 5 , c 6 > 0, independent of the meshsize, such that where the set ω F is given by Proof See [24, Lemmas 4.9 and 4.10, respectively].

Lemma 16
Let ζ h ∈ L 2 ( ) be an element-wise polynomial of degree k ≥ 0, and let ζ ∈ L 2 ( ) be such that div(ζ ) = 0 in . Then, there exist c 7 , c 8 > 0, independent of the meshsize, such that Proof Indeed, applying the first inequality given in (3.25), using that div(ζ ) = 0 in , integrating by parts, and then employing the Cauchy-Schwarz inequality, we get Now, using the inverse inequality (3.26), and the fact that 0 ≤ ψ T ≤ 1 in T , we find that which together with (3.31) gives (3.29). The proof of (3.30) corresponds to a slight adaptation of the proof of [7,Lemma 4.6], which makes use of (3.29).
After these preliminary results, we are ready to give local efficiency estimates for several terms associated to the interface.

Lemma 18
There exist constants c i > 0, i ∈ {12, 13}, independent of the meshsizes, such that where, given F ∈ F h ( ), T F is the tetrahedron of T B h (respectively T D h ) having F as a face.
Proof The proofs of (a) and (b) follow after a straightforward adaptation of that of [25,Lemma 20], and recalling from [13,Lemma 3.6] that the operator curl s is bounded.
We remark that estimates (a) and (b) provided by the previous lemma are the only nonlocal bounds of the efficiency analysis. However, under an additional regularity assumption on λ we are able to prove the following local bounds.

Lemma 19
Assume that λ| F ∈ H 1 (F), for each F ∈ F h ( ). Then there exist c 14 , c 15 > 0, independent of the meshsizes, such that for each F ∈ F h ( ) there hold Proof The derivation of these estimates follows as in the proof of [25, Lemma 21].
The following three lemmas provide the corresponding upper bounds for the remaining terms defining 2 B,T (cf. 3.5) and 2 D,T ( cf. 3.6).

Lemma 20
There exist positive constants c i , i ∈ {16, 17, 18, 19}, independent of the meshsizes, such that h having F as a face, and to derive (a) and (c) it suffices to apply the estimates (3.27) and (3.28) (cf. Lemma 15), respectively, to On the other hand, reasoning similarly as in the proof of [29, Lemma 5.14], we get from which, it is easy to deduce the estimate (b). For the proof of (d), we set ζ and ζ h as in (3.32). Given F ∈ F h ( ) we denote χ F := ζ h × n on F. Then, applying the second inequality given in (3.25), and the extension operator L : C(F) → C(T ), we find that

Now, integrating by parts, it follows that
Next, applying the Cauchy-Schwarz inequality, the inverse estimate (3.26), and the preliminary bound for curl(ζ h ) 0,T F (cf. 3.27), we deduce that In turn, recalling that 0 ≤ ψ F ≤ 1 in F, and employing the third inequality in (3.25), we can write Finally, using (3.35), the definitions of ζ and ζ h (cf. 3.32), the preliminary estimate (3.33), and the fact that h F ≤ h T F , we deduce from (3.34) that which gives (d), and ends the proof.

Lemma 21
There exist c i > 0, i ∈ {20, 21, 22, 23}, independent of the meshsizes, such that Proof Thanks to the fact that curl( f D − κ −1 D u D ) = curl(∇ p D ) = 0 in D , (a) and (c) can be obtained applying (3.27) and (3.28), respectively, to ζ : The remaining estimates follow analogously to the proofs of (b) and (d) in Lemma 20.
We next turn to the derivation of local efficiency estimates for the residual expressions defining 2 B,T .

Lemma 22
There exist c i > 0, i ∈ {24, 25, 26}, independent of the meshsizes, such that Proof Regarding (a), let us denote χ T := ω B h −curl u B h in a generic T ∈ T h . Applying the first estimate of (3.25) to χ T , and then using that curlu B = ω B in B , we find that Next, integrating by parts in the first term on the right hand side of the last identity, and recalling that ψ T vanishes on ∂ T , we obtain Then, applying Cauchy-Schwarz inequality and (3.26), we deduce from (3.36) that In this way, using that 0 ≤ ψ T ≤ 1 in T , we get which gives (a). Estimate (b) can be derived by adapting the arguments in the proof of [2,Lemma 4.15]. Finally, since div(ω B ) = div(curl u B ) = 0 in B , for the derivation of (c), it suffices to apply (3.30) to ζ := ω B and ζ h := ω B h .
We end this section by observing that the term h F λ − λ h 2 0,F appearing in Lemma 17 (items (a) and (b)), is bounded as follows: Therefore the efficiency of is a direct consequence of (3.24) and Lemmas 17-21.

Numerical results
In this section we provide two computational tests aimed at illustrating the properties of the estimator introduced in Sect. 3.2. All linear systems are solved with the distributed multifrontal direct solver MUMPS.
The domains are discretised into a series of nested uniform triangulations, where errors, experimental convergence rates, and effectivity indexes will be defined as with e andê denoting errors associated to two consecutive meshes of sizes h andĥ, and being associated to methods having N andN degrees of freedom, respectively. The first two parts of Table 1 show optimal convergence for all fields under either adaptive or uniform mesh refinement. Secondly, we regard the same domains but manufacture an exact pressure that is singular near one wall of D (the singularity being located at (x a , x b , x c ) = (0, 0, −0.55)): We expect that the convergence is hindered by the lower regularity of the exact solution. This is indeed evidenced in the third block of Table 1, where we see an oscillating effectivity index and a very low convergence, especially so for the Darcy pressure and the Lagrange multiplier. The optimal character of the error decay is however restored when we use an adaptive mesh refinement strategy (see the last section of the table). We also confirm that the error indicator performs well even if the fluid viscosity ν has a considerable variation (see Table 2). Intermediate adapted meshes and some components of the approximate solution are displayed in Fig. 1.
Example 2 Next we turn to the simulation of the flow behaviour within a composite domain = (0, 2) × (0, 0.2) × (0.75). A smooth interface exists between the Darcy and Brinkman subdomains, where the Brinkman part is on top (see related test cases in [1,4,10,14]). For this problem we assume a uniform current flow on the x 1 −direction and the presence of gravity, so f B = f D = (0.25, 0, −0.1) T . In addition, we take a dimensional parameters specified as κ −1 B = 1, ν = 0.01, and κ −1 where η is the sum of characteristic functions on 20 balls of radius 1E-4, located randomly in D and representing obstacles of much lower permeability.
The boundary conditions are set as follows: on the face x 1 = 0 we impose a unitary normal inflow velocity u B · n = 1, on the bottom and top faces we set slip velocity conditions u B · n = 0 and u D · n = 0 (for the Brinkman and Darcy boundaries, respectively), and on the remaining parts of the boundary we do not force velocity nor pressure. On the interface we impose zero tangential vorticity and the transmission conditions analyzed in the paper. We now use the method based on the RT 1 − ND 2 − RT 1 − P 1 − P 1 − P 2 family, and a penalisation approach is used to impose zero-mean value of the Brinkman pressure. In Fig. 2 we present the sketch of the domains and Table 1 Test 1: error history and convergence rates achieved by the lowest-order method