Abstract

In this paper, we develop an a posteriori error analysis for the stationary Stokes-Darcy coupled problem approximated by conforming the finite element method on isotropic meshes in , . The approach utilizes a new robust stabilized fully mixed discretization. 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.

1. 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 [1] 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 treatment 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 [2], modified by Saffman in [3], and later mathematically justified in [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 example, the works [816] for some mathematical analysis. Two different mixed formulations have been studied by Layton et al. in [15]. The first one enforces the weak continuity of the normal component of the velocity field on the interface, by employing a Lagrange multiplier, while the second one imposes the strong continuity in the functional space. We can 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 singly. The authors in [911, 15, 1719] employ the weakly coupled formulation. The works on the strongly coupled formulation have been based on the development of a unified discretization, namely, 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 Stokes problems. The authors in [8, 20] have proposed a conforming, unified finite element for the strongly coupled mixed formulation. Superconvergence analysis of the finite element methods for the Stokes-Darcy system was studied by Chen et al. in [21]. Other less restrictive discretizations as the nonconforming unified approach [12, 19] or the discontinuous Galerkin (DG) approach have been analyzed in [13, 14, 22]. Due to its discontinuous nature, some (DG) discretizations for the coupled Darcy-Stokes problem may break the strong coupling in the discrete level [13, 14], as they impose the normal continuity across the interface via interior penalties.

The adaptive techniques have become indispensable tools and unavoidable in the field of study behavior of the error committed during solving partial differential equations (PDE). 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 equidistribute the computational effort and optimize the approximation efficiency. Since the pioneering work of Babuška and Rheinboldt [23], adaptive finite element methods based on a posteriori error estimates have been extensively investigated.

A posteriori error estimations have been widely studied for both the mixed formulations of the Darcy flow [2426] and the Stokes flow [2736]. However, only a few works exist for the coupled Darcy-Stokes problem, see for instance [18, 3740]. The works in [18, 38] concern the strongly coupled mixed formulation where a conforming and nonconforming finite element methods have been employed. The papers [37, 39] concern the weakly coupled mixed formulation while [39] uses the primal formulation on the Darcy side. The authors in [40] employ a fully mixed formulation where Raviart-Thomas elements have been used to approximate the velocity in both the Stokes domain and the Darcy domain, and constant piecewise for approximating the pressure.

In [16], 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-known MINI elements () to approximate the velocity and pressure in the conduit for the Stokes equation. To capture the fully mixed technique in the porous medium region linear Lagrangian elements, has been used for hydraulic (piezometric) head and Brezzi-Douglas-Marini (BDM1) piecewise constant finite elements have been used for Darcy velocity. An a priori error analysis is performed with some numerical tests confirming the convergence rates. However, to our best knowledge, they did not talk about the adaptative method for the fully mixed discretization proposed in [16]. In this case, our main objective is to perform an a posteriori error analysis by constructing reliable and efficient indicator errors. The a posteriori error estimate is based on a suitable evaluation on the residual of the finite element solution. We prove that our indicator errors are efficiency and reliability and then are optimal. The difference between our paper and the reference [40] is that our discretization uses MINI elements () to approximate the velocity and pressure in the conduit for Stokes equations; -Lagrange elements to approximate hydraulic (piezometric) head and Brezzi-Douglas-Marini () piecewise constant finite elements have been used for Darcy velocity. As a result, an additional term is included in the error estimator that measures the stability of the method. In order to treat appropriately this stability term, we further need a special Helmholtz decomposition ([18], Theorem 3), a regularity result ([18], Theorem 4), and an estimate of the stability error ([18], Theorem 5).

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 Verfürth [41] and used in a similar context by Carstensen [25, 42]. In Section 3, the a posteriori error estimates are derived. We offer our conclusion and the further works in Section 4.

2. Preliminaries and Notations

2.1. Model Problem

We consider the model of a flow in a bounded domain , consisting of a porous medium domain , where the flow is a Darcy flow, and an open region , where the flow is governed by the Stokes equations. The two regions are separated by an interface Let , . Each interface and boundary is assumed to be polygonal or polyhedral . We denote by (resp., ) the unit outward normal vector along (resp., ). Note that on the interface , we have . Figure 1 shows a sketch of the problem domain, its boundaries, and some other notations.

The fluid velocity and pressure and are governed by the Stokes equations in : where denotes the stress tensor, and represents the deformation tensor. The porous media flow is governed by the following Darcy equations on through the fluid velocity and the piezometric head :

We impose impermeable boundary conditions, on , on the exterior boundary of the porous media region, and no-slip conditions, on , 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 completely 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 [27], interfacial coupling:

This is a simplification of the original and more physically realistic Beavers-Joseph conditions (in which in (2.8) is replaced by ; see [2]). Here, we denote and are the body forces in the fluid region and source in the porous region, is the symmetric positive define (SPD) hydraulic conductivity tensor, and is the constant parameter.

We shall also assume that all material and fluid parameters defined above are uniformly positive and bounded, i.e.,

2.2. Notations and the Weak Formulation

In this part, we first introduce some Sobolev spaces [43] and norms. If is a bounded domain of and is a nonnegative integer, the Sobolev space is defined in the usual way with the usual norm and seminorm . In particular, and we write for . Similarly, we denote by the or inner product. For shortness, if is equal to , we will drop the index , while for any , , and , for . The space denotes the closure of in . Let be the space of vector-valued functions with components in . The norm and the seminorm on are given by

For a connected open subset of the boundary , we write for the 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 and , we define the following norms:

The variational formulation of the steady-state Stokes-Darcy problem (1)–(5) reads as find satisfying where the bilinear forms are defined as

After introducing, for and , the weak formulation (12)–(15) can be equivalently rewritten as follows: find satisfying

It is easy to verify that this variational formulation is well-posedness [16].

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 , , and , only depending on such that for all ,

Besides, there exists a constant that only depends on such that for all ,

2.3. Fully Mixed Isotropic Discretization

First, we consider the family of triangulations of , consisting of and , which are regular triangulations of and , respectively, where is a positive parameter. We also assume that on the interface the two meshes of and , which form the regular triangulation , coincide.

The domain of the uniformly regular triangulation is such that and . There exist positive constants and satisfying . To approximate the diameter of the triangle (or tetrahedral) , is the diameter of the greatest ball included in . Based on the subdivisions and , we can define finite element spaces , , , and . We consider the well-known MINI elements to approximate the velocity and the pressure in the conduit for Stokes equations [44]. To capture the fully mixed technique in the porous medium region linear Lagrangian elements, are used for hydraulic (piezometric) head and Brezzi-Douglas-Marini (BDM1) piecewise constant finite elements are used for Darcy velocity [45].

In the fluid region, we select for the Stokes problem the finite element spaces that satisfy the velocity-pressure inf-sup condition: there exists a constant , independent of , such that,

In the porous region, we use the finite element spaces that also satisfy a standard inf-sup condition: there exist a constant such that for all ,

Then, the finite element discretization of (18) is to find such that

This is the natural discretization of the weak formulation (18) except that the stabilized term is added. This bilinear form is defined by

We are now able to define the norm on :

We have the following results (see [16], Theorem 2 and Theorem 3):

Theorem 1. There exists a unique solution to problem (23) and if the solution of the continuous problem (18) is smooth enough, then we have Below, in order to avoid excessive use of constants, the abbreviation stand for , with a positive constant independent of , , and .

Remark 2. (Galerkin orthogonality relation). Let be the exact solution and be the finite element solution. Then, for any , and using technical regularity result Theorem 4 below, we can subtract (18) to (23) 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 errors in the hydraulic and Darcy velocity equations are, respectively, defined by

3. A Posteriori Error Analysis

3.1. Some Technical Results

Our a posteriori analysis requires some analytical results that are recalled. We define the space with the norm

The first one concerns a sort of Helmholtz decomposition of elements of . Recall first that if ,

Theorem 3. (see Ref. [18], page 708). Any admits the Helmholtz-type decomposition where but satisfying , where if , while if , with the estimate

The second result that we need is a regularity result for the solution of (18) is the following theorem:

Theorem 4. (see [18], page 710). Let be the unique solution of (18). If , then there exists such that

Let us finish this section by an estimation of the stability error (see [18], Theorem 5):

Theorem 5. For any , we have

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 be the exact solution and be the finite element solution. Then, for any and , using the Helmholtz decomposition (Theorem 3), we have where

3.2.2. Residual Error Estimators

Definition 6 (a posteriori error indicators). The residual error estimator is locally defined by where The global residual error estimator is given by Furthermore, denote the local and global approximation terms by where the global function is defined by while in , we take for all , as the unique element of such that

Remark 7. The residual character of each term on the right-hand sides of (43) and (44) is quite clear since if would be the exact solution of (18), then they would vanish.

3.2.3. Analytical Tools

(1) Inverse Inequalities. In order to derive the lower error bounds, we proceed similarly as in [25, 42] (see also [46]), 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 , and , we let and be the usual simplex-bubble and face-bubble functions, respectively (see (1.5) and (1.6) in [41]). In particular, satisfies , , , and . Similarly, , , , and . We also recall from [47] that, given , there exists an extension operator that satisfies and . A corresponding vectorial version of , that is, the componentwise application of , is denoted by . Additional properties of , , and are collected in the following lemma (see [47]).

Lemma 8. Given , there exist positive constants depending only on and shape-regularity of the triangulations (minimum angle condition), such that for each simplex and , there hold

(2) Continuous Trace Inequality.

Lemma 9 (continuous trace inequality). There exists a positive constant depending only on such that

(3) Clément Interpolation Operator. In order to derive the upper error bounds, we introduce the Clément interpolation operator that approximates optimally nonsmooth functions by continuous piecewise linear functions:

In addition, we will make use of a vector-valued version of , that is, , which is defined componentwise by The following lemma establishes the local approximation properties of (and hence of ), for a proof see [48], Section 3.

Lemma 10. There exist constants , independent of , such that for all , there hold where and .

3.3. Optimality of
3.3.1. Reliability Result

Theorem 11 (reliability of ). Let be the exact solution of (18) and be the finite element solution of (23). There exist a constant such that the following estimate holds:

Proof. We take in error equation (38)–(41) and we obtain where The inf-sup condition of leads to Now, using the error equation (57)–(59), Cauchy-Schwarz inequality, and the Clément operator of Lemma 10, we deduce the estimate (56). The proof is complete.

3.3.2. Efficiency Result

To prove local efficiency for , let us denote by where

The main result of this subsection can be stated as follows.

Theorem 12 (efficiency of ). Under the assumptions of Theorem 4, the following lower error bound holds: where is a finite union of neighboring elements of .

Proof. We begin by bounding each term of the residuals separately. (i)Element residual in : to estimate , we choose in error equation (57)–(59) for each , and with on ,for obtained, Noted that and vanish on . Because in this case, then we have The first inverse inequality (49) and the Cauchy-Schwarz inequality lead to Using inverse inequality (50), we deduce the estimate: (ii)Divergence element residual in (estimation of ): for each , we have,(iii)Element residual in : we have for each ,(iv)Curl element residual in : for , we set and . Hence, we notice that belongs to and is divergence free; therefore, by Equations (57)–(59), we obtain with and , . The first inverse inequality (49) and the Cauchy-Schwarz inequality lead toAgain, the inverse inequality (49) allows to get (v)Interface elements on : we fix an edge included in and for a constant fixed later on a unit vector , we considerthat clearly belongs to . Hence, by residual equation (57)–(59), we obtain with , where (resp., ) is the unique triangle/tetrahedron included in (resp., ) having as edge/face, and Taken in , in and for each , with . We have Hence, Inverse inequalities (50) and (51) and Cauchy-Schwarz inequality lead to with .
Taken in , in and for each , with . As before the identities (72)–(74) and the inverse inequalities (50) and (51) lead to (vi)Piezometric head jump in : for each edge/face , we consider . As , we setUsing the residual equation (57)–(59), we obtain with and where : Therefore, Cauchy-Schwarz inequality and inverse inequalities (51) and (52) lead to In additionally, since by regularity Theorem 4, the jump of is zero through all the edges of ; hence, we clearly have . Thus, (vii)Finally, for , we haveThus, The estimates (66), (67), (70), (76), (77), (81), and (83) provide the desired local lower error bound.

4. Summary

In this paper, we have discussed a posteriori error estimates for a finite element approximation of the Stokes-Darcy system. A residual type a posteriori error estimator is provided that is both reliable and efficient. Many issues remain to be addressed in this area, let us mention other types of a posteriori error estimators or implementation and convergence analysis of adaptive finite element methods.

5. Nomenclature

(i) bounded domain(ii) the porous medium domain(iii) the fluid region(iv)(v)(vi) (vii) (resp., ) the unit outward normal vector along (resp., )(viii): the fluid velocity(ix): the fluid pressure(x)In 2D, the of a scalar function is given as usual by(xi)In 3D, the of a vector function is given as usual by , namely,(xii): the space of polynomials of total degree not larger than (xiii): triangulation of (xiv): the corresponding induced triangulation of , (xv)For any , is the diameter of and is the diameter of the largest ball inscribed into (xvi) and (xvii): the set of all the edges or faces of the triangulation(xviii): the set of all the edges () or faces () of a element (xix)(xx): the set of all the vertices of a element (xxi)(xxii)For , (xxiii)For , we associate a unit vector such that is orthogonal to and equals to the unit exterior normal vector to (xxiv)For , is the jump across in the direction of (xxv)In order to avoid excessive use of constants, the abbreviations and stand for and , respectively, with positive constants independent of , , or

Data Availability

There are no data underlying the findings in this paper to be shared.

Conflicts of Interest

The authors declare that they have no conflicts of interest.