A ug 2 01 9 A posteriori error analysis for a new fully-mixed isotropic discretization to the stationary Stokes-Darcy coupled problem

In this paper we develop an a posteriori error analysis for the stationary Stokes-Darcy coupled problem approximated by conforming finite element method on isotropic meshes in R, d ∈ {2, 3}. The approach utilizes a new robust stabilized fully mixed discretization developed in [58]. The a posteriori error estimate is based on a suitable evaluation on the residual of the finite element solution plus the stabilization terms. It is proven that the a posteriori error estimate provided in this paper is both reliable and efficient. Mathematics Subject Classification [MSC]: 74S05,74S10,74S15, 74S20,74S25,74S30.


Introduction
There are many serious problems currently facing the world in which the coupling between groundwater and surface water is important. These include questions such as predicting how pollution discharges into streams, lakes, and rivers making its way into the water supply. This coupling is also important in technological applications involving filtration. We refer to the nice overview [24] and the references therein for its physical background, modeling, and standard numerical methods. One important issue in the modeling of the coupled Darcy-Stokes flow is the treatement of the interface condition, where the Stokes fluid meets the porous medium. In this paper, we only consider the so-called Beavers-Joseph-Saffman condition, which was experimentally derived by Beavers and Joseph in [8], modified by Saffman in [51], and later mathematically justified in [36][37][38]47].
There are three popular formulations of the coupled Darcy-Stokes flow, namely the primal formulation, the mixed formulation in the Darcy region or the fully mixed formulation, see for examples [2,26,29,30,41,48,49,57,58] for some mathematical analysis. The authors in [57] studied two different mixed formulations: the first one imposes the weak continuity of the normal component of the velocity field on the interface, by using a Lagrange multiplier; while the second one imposes the strong continuity in the functional space. Later on we call these two mixed formulations, the weakly coupled formulation and the strongly coupled formulation respectively. The weakly coupled formulation gives more freedom in the choice of the discretization in the Stokes side and the Darcy side separately. The works in [26,29,30,35,44,50,57] are based on the weakly coupled formulation. Researches on the strongly coupled formulation have been focused on the development of an unified discretization, that is, the Stokes side and the Darcy side are discretized using the same finite element. This approach simplifies the numerical implementation, only if the unified discretization is not significantly more complicated than the commonly used discretizations for the Darcy and the Stokes problems. In [2,3], a conforming, unified finite element has been proposed for the strongly coupled mixed formulation. Superconvergence analysis of the finite element methods for the Stokes-Darcy system was studied in [15]. Other less restrictive discretizations as the non-conforming unified approach [41,50] or the discontinuous Galerkin (DG) approach have been proposed in [39,48,49]. Due to its discontinuous nature, some (DG) discretizations for the coupled Darcy-Stokes problem may break the strong coupling in the discrete level [48,49], as they impose the normal continuity across the interface via interior penalties.
A posteriori error estimators are computable quantities, expressed in terms of the discrete solution and of the data that measure the actual discrete errors without the knowledge of the exact solution. They are essential to design adaptive mesh refinement algorithms which equi-distribute the computational effort and optimize the approximation efficiency. Since the pioneering work of Babuska and Rheinboldt [6], adaptive finite element methods based on a posteriori error estimates have been extensively investigated.
An a priori error analysis is performed with some numerical tests confirming the convergence rates.
A posteriori error estimations have been well-established for both the mixed formulation of the Darcy flow [10,12,42], and the Stokes flow [1,7,14,22,25,34,46,52,55,56]. However, only few works exist for the coupled Darcy-Stokes problem, see for instance [5,16,21,28,44]. The paper [16,44] concern the strongly coupled mixed formulation where a H(div) conforming and nonconforming finite element methods have been used and [5,21] concern the weakly coupled mixed formulation while [21] uses the primal formulation on the Darcy side. The authors in [28] employ a fully-mixed formulation where Raviart-Thomas elements have been used to approximate the velocity in both the Stokes domain and Darcy domain, and constant piecewise for approximate the pressure.
In [58], a stabilized finite element method for the stationary mixed Stokes-Darcy problem has been proposed for the fully-mixed formulation. The authors have used the well-know MINI elements (P 1b − P 1) to approximate the velocity and pressure in the conduit for Stokes equation. To capture the fully mixed technique in the porous medium region linear Lagrangian elements, P 1 have been used for hydraulic (piezometric) head and Brezzi-Douglas-Marini (BDM 1) piecewise constant finite elements have been used for Darcy velocity. To our best knowledge, there is no a posteriori error estimation for the fully-mixed discretization proposed in [58]. Here we develop such a posteriori error analysis. The a posteriori error estimate is based on a suitable evaluation on the residual of the finite element solution. We further prove that our a posteriori error estimator is both reliable and efficient. The difference between our paper and the reference [28] is that our discretization uses MINI elements (P 1b − P 1) to approximate the velocity and pressure in the conduit for Stokes equations, P 1-Lagrange elements to approximate hydraulic (piezometric) head and Brezzi-Douglas-Marini (BDM 1) piecewise constant finite elements have been used for Darcy velocity. As a result, additional term is included in the error estimator that measure the stability of the method. In order to treat appropriately this stability term, we further need a special Helmholtz decomposition [44,  The paper is organized as follows. Some preliminaries and notation are given in section 2. The efficiency result is derived using the technique of bubble function introduced by R. Verfürth [54] and used in similar context by C. Carstensen [12,13]. In section 3, the a posteriori error estimates are derived.

Preliminaries and Notations
2.1. Model problem. We consider the model of a flow in a bounded domain Ω ⊂ R d (d = 2 or 3), consisting of a porous medium domain Ω p , where the flow is a Darcy flow, and an open region Ω f = Ω Ω p , where the flow is governed by the Stokes equations. The two regions are separated by an interface Γ = ∂Ω p ∩ ∂Ω f . Let Γ l = ∂Ω l Γ, l = f, p. Each interface and boundary is assumed to be polygonal (d = 2) or polyhedral (d = 3). We denote by n f (resp. n p ) the unit outward normal vector along ∂Ω f (resp. ∂Ω p ). Note that on the interface Γ, we have n f = −n s . The Figure 1 shows a sketch of the problem domain, its boundaries and some other notations. The fluid velocity and pressure u f (x) and p(x) are governed by the Stikes equations in Ω f :

Ωp: Porous Medium
where T = −pI + 2νD(u f ) denotes the stress tensor, and D(u f ) = 1 2 ∇u f + (∇u h ) T represents the deformation tensor. The porous media flow is governed by the following Darcy equations on Ω p through the fluid velocity u p and the piezometric head φ(x): We impose impermeable boundary conditions, u p · n p = 0 on Γ p , on the exterior boundary of the porous media region, and no slip conditions, u f = 0 on Γ f , in the Stokes region. Both selections of boundary conditions can be modified. On Γ the interface coupling conditions are conservation of mass, balance of forces and a tangential condition on the fluid region's velocity on the interface. The correct tangential condition is not competely understood (possibly due to matching a pointwise velocity in the fluid region with an averaged or homogenized velocity in the porous region). In this paper, we take the Beavers-Joseph-Saffman (-Jones), see [8,[36][37][38]47,51], interfacial coupling: This is a simplification of the original and more physically relistic Beavers-Joseph conditions (in which u f · τ j in (2.8) is replaced by (u f − u p ) · τ j ; see [8] ).
Here we denote f f , f p -body forces in the fluid region and source in the porous region, K-symmetric positive define (SPD) hydraulic conductivity tensor, α-constant parameter. We shall also assume that all material and fluid parameters defined above are uniformly positive and bounded, i.e., 0 ≤ k min ≤ λ(K) ≤ k max < ∞.

2.2.
Notations and the weak formulation. In this part, we first introduce some Sobolev spaces [43] and norms. If W is a bounded domain of R d and m is a non negative integer, the Sobolev space H m (W ) = W m,2 (W ) is defined in the usual way with the usual norm · m,W and semi-norm |.| m,W . In particular, H 0 (W ) = L 2 (W ) and we write · W for · 0,W . Similarly we denote by (·, ·) W the L 2 (W ) [L 2 (W )] N or [L 2 W )] d×d inner product. For shortness if W is equal to Ω, we will drop the index Ω, while for any m ≥ 0, · m,l = · m,Ω l , |.| m,l = |.| m,Ω l and (., For a connected open subset of the boundary Γ ⊂ ∂Ω s ∪ ∂Ω d , we write ., . Γ for the L 2 (Γ) inner product (or duality pairing), that is, for scalar valued functions λ, σ one defines: By setting the space we introduce the following spaces: For the spaces X f and X p , we define the following norms: Ωp , ∀v p ∈ X p . The variational formulation of the steady-state Stokes-Darcy problem (1)-(5) reads as: Find (u f , p; u p , φ) ∈ (X f , Q f ; X p , Q p ) satisfying: where the bilinear forms are defined as: the weak formulation (8)- (11) can be equivalentl rewritten as follows: Find U ∈ H satisfying It is easy to verify that this variational fromulation is well-posedness.
To end this section, we recall the following Poincaré, Korn's and the trace inequalities, which will be used in the later analysis; There exist constant C p , Besides, there exists a constantC v that only depends on Ω p such that for all ψ ∈ Q p , 2.3. Fully-mixed isotropic discretization. First, we consider the family of triangulations T h of Ω, consisting of T f h and T p h , which are regular triangulations of Ω f and Ω p , respectively, where h > 0 is a positive parameter. We also assume that on the interface Γ the two meshes of T f h and T p h , which form the regular triangulation T h := T f h ∪ T p h , coincide. The domain of the uniformly regular triangulation Ω f ∪ Ω p is such that Ω = {∪K : K ∈ T h } and h = max K∈T h h K . There exist positive constants c 1 and c 2 satisfying c 1 h ≤ h K ≤ c 2 ρ K . To approximate the diameter h K of the trangle (or tetrahedral) K, ρ K is the diameter of the greatest ball included in K. Based on the subdivisions T f h and T p h , we can define finite element spaces We consider the well-known MINI elements (P 1b − P 1) to approximate the velocity and the pressure in the conduit for Stokes equations [4]. To capture the fully-mixed technique in the porous medium region linear Lagrangian elements, P 1 are used for hydraulic (piezometric) head and Brezzi-Douglas-Marini (BDM 1) piecewise constant finite elements are used for Darcy velocity [11]. In the fluid region, we select for the Stokes problem the finite element spaces (X f h , Q f h ) that satisfy the velocity-pressure inf-sup condition: There exists a constant β f > 0, independent of h, such that, In the porous region, we use the finite element spaces (X ph , Q ph ) that also satisfy a standard inf-sup condition: There exist a constant β p > 0 such that for all Then the finite element discretization of (14) is to find This is the natural discretization of the weak formulation (14) except that the stabilized term J Γ (U h , V h ) is added. This bilinear form J Γ (., .) is defined by We are now able to define the norm on H h : We have the following results (see [58, (18) and if the solution U ∈ H of the continuous problem (14) is smooth enough, then we have: Below, in order to avoid excessive use of constants, the abbreviation x y stand for x cy, with c a positive constant independent of x, y and T h .  (14) to (18) to obtain the Galerkin orthogonality relation: thus, we have the relation: where here and below, the errors in the velocity and in the pressure of Stokes equations, and erors in the hydraulic and Darcy velocity equations are respectively defined by:

A posteriori error analysis
3.1. Some technical results. Our a posteriori analysis requires some analytical results that are recalled. We define the space Ωp . The first one concerns a sort of Helmholtz decomposition of elements of H. Recall first that if d = 3, H 0 (curl, Ω p ) = {ψ ∈ L 2 (Ω p ) 3 : curl ψ ∈ L 2 (Ω p ) 3 and ψ × n = 0 on ∂Ω p }.
The second result that we need is a regularity result for the solution U = (u f , p, u p , φ) ∈ H of (14) is the following theorem:  (14). If K ∈ [C 0,1 (Ω p )] d×d , then there exists ǫ > 0 such that Let us finish this section by an estimation of the stability error (see [44,Theorem 3.3]): 3.2. Error estimator. In order to solve the Stokes-Darcy coupled problem by efficient adaptive finite element methods, reliable and efficient a posteriori error analysis is important to provide appropriated indicators. In this section, we first define the local and global indicators and then the lower and upper error bounds are derived in Section 3.3.
3.2.1. Error equations. The general philosophy of residual error estimators is to estimate an appropriate norm of the correct residual by terms that can be evaluated easier, and that involve the data at hand. Thus we define the error equations: Let U = (u f , p, u p , φ) ∈ H be the exact solution and U h = (u f h , p h , u ph , φ h ) ∈ H h be the finite element solution. Then for any (1) Error indicators in Stokes domain : For K ∈ T s h , we set:

Residual Error Estimators.
(2) Error indicators in Darcy domain: For K ∈ T p h , we set: The residual error estimator is locally defined by: The global residual error estimator is given by: Furthermore denote the local and global approximation terms by ∀K ∈ T p h , and Remark 3.1. The residual character of each term on the right-hand sides of (28) and (29) is quite clear since if (u f h , p h , u ph , φ h ) ∈ H h would be the exact solution of (14), then they would vanish.

Analytical tools.
(1) Inverse inequalities: In order to derive the lower error bounds, we proceed similarly as in [12] and [13] (see also [27]), by applying inverse inequalities, and the localization technique based on simplex-bubble and face-bubble functions. To this end, we recall some notation and introduce further preliminary results. Given K ∈ T h , and E ∈ E(K), we let b K and b E be the usual simplexe-bubble and face-bubble functions respectively (see (1.5) and (1.6) in [54]). In particular, b We also recall from [53] that, given k ∈ N, there exists an extension operator L : C(E) −→ C(K) that satisfies L(p) ∈ P k (K) and L(p) |E = p, ∀p ∈ P k (E). A corresponding vectorial version of L, that is, the componentwise application of L, is denoted by L. Additional properties of b K , b E and L are collected in the following lemma (see [53]) Lemma 3.1. Given k ∈ N * , there exist positive constants depending only on k and shape-regularity of the triangulations (minimum angle condition), such that for each simplexe K and E ∈ E(K) there hold q K qb 1/2 K K q K , ∀q ∈ P k (K) (33) |qb K | 1,K h −1 K q K , ∀q ∈ P k (K) (34) (2) Continuous trace inequality where ∆(K) := ∪ {K ′ ∈ T h : K ′ ∩ K = ∅} and ∆(E) := ∪ {K ′ ∈ T h : K ′ ∩ E = ∅}.

Efficiency and Reliability.
Theorem 3.4. Let U = (u f , p, u p , φ) ∈ H be the exact solution of (14) and U h = (u f h , p h , u ph , φ h ) ∈ H h be the finite element solution of (18). Then there holds: