A new unified stabilized mixed finite element method of the Stokes-Darcy coupled problem: Isotropic discretization

In this paper we develop an a priori error analysis of a new unified mixed finite element method for the coupling of fluid flow with porous media flow in $\mathbb{R}^N$, $N\in\{2,3\}$ on isotropic meshes. Flows are governed by the Stokes and Darcy equations, respectively, and the corresponding transmission conditions are given by mass conservation, balance of normal forces, and the Beavers-Joseph-Saffman law. The approach utilizes a modification of the Darcy problem which allows us to apply a variant nonconforming Crouzeix-Raviart finite element to the whole coupled Stokes-Darcy problem. The well-posedness of the finite element scheme and its convergence analysis are derived. Finally, the numerical experiments are presented, which confirm the excellent stability and accuracy of our method.

into the water supply. This coupling is also important in technological applications involving filtration. We refer to the nice overview [9] 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 [4], modified by Saffman in [31], and later mathematically justified in [20,21,23,28].
It is well known that the discretization of the velocity and the pressure, for both Stokes and Darcy problems and the coupled of them, has to be made in a compatible way in order to avoid instabilities. Since, usually, stable elements for the free fluid flow cannot been successfully applied to the porous medium flow, most of the finite element formulations developed for the Stokes-Darcy coupled problem are based on appropriate combinations of stable elements for the Stokes equations with stable elements for the Darcy equations. In [1-3, 6, 12, 13, 15, 18-20, 22-27, 29, 30, 32-34], and in the references therein, we can find a large list of contributions devoted to numerically approximate the solution of this interaction problem, including conforming and nonconforming methods.
There are a lot of papers considering different finite element spaces in each flow region (see, for example, [8,13,14] and the references therein). In contrast to this, other articles use the same finite element spaces in both regions by, in general, introducing some penalizing terms (ref. for examples [2,27,30] and the references therein).
In [2], a conforming unified finite element has been proposed for the modified coupled Stokes-Darcy problem in a plane domain, which has simple and straightforward implementations. The authors apply the classical Mini-element to the whole coupled Stokes-Darcy coupled problem. An a priori error analysis is performed with some numerical tests confirming the convergence rates.
In this article, we propose a modification of the Darcy problem which allows us to apply a variant nonconforming finite element to the whole coupled Stokes-Darcy problem. We use a variant nonconforming Crouzeix-Raviart finite element method that has so many advantages for the velocities and piecewise constant for the pressures in both the Stokes and Darcy regions, and apply a stabilization term penalizing the jumps over the element edges of the piecewise continuous velocities. We prove that the formulation satisfies the discrete inf-sup conditions, obtaining as a result optimal accuracy with respect to solution regularity. Numerical experimants are also presented, which confirm the excellent stability and optimal performance of our method. The difference between our paper and the reference [2] is that our discretization is nonconforming in both the Stokes domain and Darcy domain (in Ω ⊂ R N , N = 2 or 3). As a result, additional terms are included in the priori error analysis that measure the non-conformity of the method. One essential difficulty in choosing the unified discretization is that, the Stokes side velocity is in H 1 while the Darcy side velocity is only in H(div). Thus, we introduce a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element space (larger than the space H h used in [30]). The choice of H h [see (28)] is more natural than the one introduced in [30] since the space H h approximates only H(div, Ω d ) and not [H 1 (Ω d )] N , while our a priori error analysis is only valid in this larger space.
The rest of the paper is organized as follows. In Section 2 we present the modified coupled Stokes-Darcy problem in Ω ⊂ R N , N = 2 or 3, notations and the weak formulation. Section 3 is devoted to the finite element discretization and the error estimation. Finally, in Section 4, we present the results of numerical experiments to verify the predicted rates of convergence.
Each interface and boundary is assumed to be polygonal (N = 2) or polyhedral (N = 3). We denote by n s (resp. n d ) the unit outward normal vector along ∂Ω s (resp. ∂Ω d ). Note that on the interface Γ I , we have n s = −n d . The Figures 1 and 2 give a schematic representation of the geometry.  For any function v defined in Ω, since its restriction to Ω s or to Ω d could play a different mathematical roles (for instance their traces on Γ I ), we will set v s = v |Ωs In Ω, we denote by u the fluid velocity and by p the pressure. The motion of the fluid in Ω s is described by the Stokes equations while in the porous medium Ω d , by Darcy's law Here, µ > 0 is the fluid viscosity, D the deformation rate tensor defined by and K a symmetric and uniformly positive definite tensor representing the rock permeability and satisfying, for some constants 0 < K * K * < +∞, f ∈ [L 2 (Ω)] N is a term related to body forces and g ∈ L 2 (Ω) a source or sink term satisfying the compatibility condition Ω g(x)dx = 0.
Finally we consider the following interface conditions on Γ I : Here, Eq. (3) represents mass conservation, Eq. (4) the balance of normal forces, and Eq. (5) the Beavers-Joseph-Saffman conditions. Moreover, {τ j } j=1,...,N −1 denotes an orthonormal system of tangent vectors on Γ I , κ j = τ j · K · τ j , and α 1 is a parameter determined by experimental evidence. Eqs. (1) to (5) consist of the model of the coupled Stokes and Darcy flows problem that we will study below.
2.2. New weak formulation. We begin this subsection by introducing some useful notations. If W is a bounded domain of R N 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 )] N ×N 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 (., .) l = (·, ·) Ω l , for l = s, d.
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: We also define the special vector-valued functions space To give the variational formulation of our coupled problem we define the following two spaces for the velocity and the pressure: and Q = L 2 0 (Ω) := q ∈ L 2 (Ω) : Multiplying the first equation of (1) by a test fonction v ∈ H and the second one by q ∈ Q, integrating by parts over Ω s the terms involving div D(u) and ∇p, yield the variational form of Stokes equations: Using interface conditions (4) and (5) in (11), we obtain: We apply a similar treatment to the Darcy equations by testing the first equation of (2) with a smooth fonction v ∈ H and the second on by q ∈ Q, integrating by parts over Ω d the terms involving ∇p d , yield the variational form of Darcy equations: Now, incorporating the first boundary interface condition (3) and taking into account that the vector valued functions in H have (weakly) continuous normal components on Γ I (see [16,Theorem 2.5]), the mixed variational formulation of the coupled problem (1)-(5) can be stated as follows [27]: Find (u, p) ∈ H × Q that satisfies where the bilinear forms a(·, ·) and b(·, ·) are defined on H × H and H × Q, respectively, as: By last, the linear forms L and G are defined as: It is easy to prove that a et b are continuous, b satisfies the continuous inf-sup condtion and a is coercive on the null space of b. It is also clear that F and G are continuous and bounded. Then, using the classical theory of mixed methods (see, e.g., [16, Theorem and Corollary 4.1 in Chapter I]) it follows the well-posedness of the continuous formulation (17) and so the following theorem holds [27]: (17).
and (3) hold ( the differential equations being understood in the distributional sense), while the interface conditions (4) and (5) are imposed in a weak sense. Also, we observe that the mixed variational formulation of the coupled problem (1)-(5) is equivalent to weak formulation (2.4) (and also (2.5) of [33]), with the particularity that, in our case, for any v ∈ H, we have that v s − v d , n s p s Γ I = 0. Now we introduce a modification to the Darcy equation, with the purpose in mind of the development of a unified discretization for the coupled problem, that is, the Stokes and Darcy parts be discretized using the same finite element spaces. The modification that we apply to the Darcy equation follows the idea (same argument) given in [2]. Indeed, we observe that taking the second equation of Darcy' problem (2) we can write, for any v ∈ H, Then, by adding this equation to the first equation of the variational form in (15), we get: From now on, we work with this modified variational form of Darcy equations.
In the same way that before, incorporating the boundary conditions (3) and remambering that, since v ∈ H, it was (weakly) continuous normal components on Γ I , the variational form of the modified Stokes-Darcy problem can be written as where the bilinear formsã(·, ·) and b(·, ·) are defined on H × H, H × Q, respectively, as: By last, the linear formsL and G are defined as: Then, applying the classical theory of mixed methods it follows the well-posedness of the continuous formulation (21).
There exists a unique (u, p) ∈ H × Q solution to modified formulation (21). In addition, there exists a positive constantC, depending on the continuous inf-sup condition constant for b, the coercivity constant forã and the boundedness constants forã and b, such that: We end this section with some notation. In 2D, the curl of a scalar function w is given as usual by curl w := ( ∂w ∂x 2 , − ∂w ∂x 1 ) while in 3D, the curl of a vector function w is given as usual by curl w := ∇ × w. Finally, let P k be the space of polynomials of total degree not larger than k. In order to avoid excessive use of constants, the abbreviations x y and x ∼ y stand for x cy and c 1 x y c 2 x, respectively, with positive constants independent of x, y or T h .

A priori error analysis
3.1. Finite element discretization. In this subsection, we will use a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element approximation for the velocity and piecewise constant approximation for the pressure.
Let {T h } h>0 be a family of triangulations of Ω with nondegenerate elements (i.e. triangles for N = 2 and tetrahedrons for N = 3). For any T ∈ T h , we denote by h T the diameter of T and ρ T the diameter of the largest ball inscribed into T and set We assume that the family of triangulations is regular, in the sense that there exists σ 0 > 0 such that σ h σ 0 , for all h > 0. We also assume that the triangulation is conform with respect to the partition of Ω into Ω s and Ω d , namely each T ∈ T h is either in Ω s or in Ω d (see Fig. 3 where With every edges E ∈ E h , we associate a unit vector n E such that n E is orthogonal to E and equals to the unit exterior normal vector to ∂Ω if E ⊂ ∂Ω. For any E ∈ E h and any piecewise continuous function ϕ, we denote by [ϕ] E its jump across E in the direction of n E : For i ∈ {0, · · · , N }, we set: The where for each i ∈ {0, · · · , N }, λ i (T ) is barycentric coordonates of T ∈ T h . In classical reference element T , the basis fonctions are given by:   Based on the above notation, we introduce a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element space (larger than the space H h used in [30]) and piecewise constant function space is the space of the restrictions to T of all polynomials of degree less than or equal to m. The space Q h is equipped with the norm · while the norm on H h will be specified later on. The choice of H h is more natural than the one introduced in [30] since the space H h approximates only H(div, Ω d ) and not [H 1 (Ω d )] N , while our a priori error analysis is only valid in this larger space.
Let us introduce the discrete divergence operator div h ∈ L(H h ; Q h ) ∩ L(H; Q) by Then, we can introduce two bilinear forms Then the finite element discretization of (21) is to find This is the natural discretization of the modified weak formulation (21) except that the penalizing term J(u h , v h ) is added. This bilinear form J(., .) is defined by following the decomposition (32) of E h : where Here, h E is the length (N = 2) or diameter (N = 3) of E. Note that each element of E h only contributes with one jump term in J(u, v).
Remark 3.1. The Eq. (31) have the matrix representation where U (resp. P) denote the coefficients of u h (resp. p h ) expanded with respect to a basis for H h (rep. Q h ).
We are now able to define the norm on H h (see [30]): In the sequel, we will denote by α, β and C i various constants independent of h. For the sake of convenience, we will define the bilinear form: From Hölder's inequality, we derive the boundedness of A h (·, ·) and b h (·, ·): Lemma 3.1. (Continuity of forms) There holds: Theorem 3.1. (Coercivity of A h ) There is an α > 0 such that: and for ψ ∈ [H 1 (T )] N , we define with the semi-norm Using Young's inequality and Green formula, we have: We have by Cauchy-Schwarz inequality: Also, we have: Then, Hence we deduce • Now we estime the term By Cauchy-Schwarz, we obtain: Thus we deduce the estimation: Then, . We apply Korn's discrete inequality [5] and we get: Thus Hence, We have, The estimates (44), (45), (46) and (47), lead to (37). The proof is complete.
In order to verify the discrete inf-sup condition, we define the space: We define also the Crouzeix-Raviart interpolation operator r h : W → H h by: Lemma 3.2. The operator r h is bounded: there is a constant C 5 > 0 depending on σ, µ and N such that Proof. The proof is similar to [30].
Then, we have the following result Proof. We use Fortin argument i.e. for each [H 1 0 (Ω)] N ⊂ W , hence v ∈ W . We take v h = r h v ∈ H h and we have: = 0 ( from the identities (49) and (50)).
Thus, we obtain Using the system (53), we have: From (54) et (55), we deduce: The Inf-Sup condition holds and the proof is complete.
From Theorem 3.1 and Theorem 3.2 we have the following result: There exists a unique solution (u h , p h ) ∈ H h × Q h to the problem (21).

A convergence analysis.
We now present an a priori analysis of the approximation error: The use of nonconforming finite element leads to H h H, so the approximation error contains some extra consistency error terms. In fact, the abstract error estimates give the following result: (21) and (u h , p h ) ∈ H h × Q h be the solution of the discrete problem (31). Then we have where E 1h and E 2h are the consistency error terms define by: For estiming the approximation error, we assume that the solution (u, p) of problem (21) satisfies the smoothness assumptions: We begin with the estimates for the terms: inf Finally, let us consider the term Thus, we havẽ where In order to evaluate the four face integrals, let us introduce two projections operators in the following. For any T ∈ T h and E ∈ E(T ), denote by P 0 (E) the constant space of the restrictions to E and π E the projection operator from L ( E) on to P 0 (E) such that The operator π E has the property [7]: For any v ∈ [L 2 (E)] N , we let Π E v be the function in [P 0 (E)] N such that Using inequality (64), we obtain Then we have the following lemma: Lemma 3.5. (Estimation the four face integrals) There holds: Proof.
(1) Estimate (66): We begin with an estimate for the first term R 1 (v h ). For any face E ∈ E h (Ω + s ), there exists at least one element T ∈ T s h such that E ∈ E(T ). Then, from condition (63), Höder's inequality and inequality (65), it follows that (2) Estimate (67): Thus, Furthermore, summing on E ∈ E h (Ω + s ) faces, we obtain the estimate: (3) For the terms R 3 (v h ) and R 4 (v h ), we use the same techniques as in the proof of the bounds for R i (v h ), i ∈ {1, 2}, and we obtain: The proof is complete.
From Lemma 3.3, Lemma 3.4 and Lemma 3.5, now we derive the following convergence theorem: Theorem 3.4. Let the solution (u, p) of problem (21) satifies the smoothness assumption (Assumption 3.1). Let (u h , p h ) be the solution of the discrete problem (31). Then there exists a positive constant C depending on N , µ, K * , K * , α 1 and σ such that: (72)

Numerical experiments
In this section we present one test case to verify the predicted rates of convergence. The numerical simulations have been performed on the finite element code FreeFem++ [11,17] in isotropic coupled mesh of Fig. 13. The solutions have been represented by Mathematica software. For simplicity we choose each domain Ω l , l ∈ {s, d} as the unit square, α 1 = µ = 1, and the permeability tensor K is taken to be the identity. The interface Γ I , is the line x = 1, i.e. Ω = [0, 1[∪{1}∪]1, 2] like the show the Figure 10. We consider the application φ : In Ω, we define u = (u 1 , u 2 ) = curl φ − ∂φ ∂y , ∂φ ∂x and we obtain: We choose quadratic pressure p ∈ L 2 (Ω) by Thus, Ω p(x, y)dxdy = 0 and ∇p = (2x − 2y, −2x + y) .

Conclusion
In this contribution, we investigated a new mixed finite element method to solve the Stokes-Darcy fluid flow model without introducing any Lagrange multiplier.    We proposed a modification of the Darcy problem which allows us to apply a slight variant nonconforming Crouzeix-Raviart element to the whole coupled Stokes-Darcy problem. The proposed method is probably one the cheapest method for Discontinuous Galerkin (DG) approximation of the coupled system, has optimal accuracy        Figure 28. Right-hand term k 2 in Ω d .
with respect to solution regularity, and has simple and straightforward implementations. Numerical experiments have been also presented, which confirm the excellent stability and accuracy of our method.

Acknowledgment
The author thanks Professor Emmanuel Creusé (University of Lille 1, France) for having sent us useful documents and for fruitful discussions concerning the numerical tests.