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 its way into the water supply. This coupling is also important in technological applications involving filtration. We refer to the nice overview  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 , modified by Saffman in , and later mathematically justified in    .
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    - , 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,    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    and the references therein).
In , 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 condition, obtaining as a result optimal accuracy with respect to solution regularity. Numerical experiments are also presented, which confirm the excellent stability and optimal performance of our method. The difference between our paper and the reference  is that our discretization is nonconforming in both the Stokes domain and Darcy domain (in , or 3). As a result, additional terms are included in the priori error analysis that measures the non-conformity of the method. One essential difficulty in choosing the unified discretization is that, the Stokes side velocity is in while the Darcy side velocity is only in . Thus, we introduce a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element space (larger than the space used in  ). The choice of [see (34)] is more natural than the one introduced in  since the space approximates only and not , 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 , or 3, notations and the weak formulation. Section 3 is devoted to the finite element discretization and the error estimation.
In Section 4, we present the results of numerical experiments to verify the predicted rates of convergence. Finally, we offer our conclusion and the further works in Section 5.
2. Preliminaries and Notation
2.1. Model Problem
We consider the model of a flow in a bounded domain ( or 3), 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 and Figure 2 give a schematic representation of the geometry.
For any function v defined in , since its restriction to or to could play a different mathematical roles (for instance their traces on ), we will set and .
Figure 1. A sketch of the geometry of the problem (case: ).
Figure 2. A sketch of the geometry of the problem (case ).
In , we denote by the fluid velocity and by p the pressure. The motion of the fluid in is described by the Stokes equations
while in the porous medium , by Darcy’s law
Here, is the fluid viscosity, the deformation rate tensor defined by
and a symmetric and uniformly positive definite tensor representing the rock permeability and satisfying, for some constants ,
is a term related to body forces and a source or sink term satisfying the compatibility condition
Finally we consider the following interface conditions on :
Here, Equation (3) represents mass conservation, Equation (4) the balance of normal forces, and Equation (5) the Beavers-Joseph-Saffman conditions. Moreover, denotes an orthonormal system of tangent vectors on , , and is a parameter determined by experimental evidence.
Equations (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 and m is a non negative integer, the Sobolev space is defined in the usual way with the usual norm and semi-norm . In particular, and we write for . Similarly we denote by the or inner product. For shortness if W 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 semi-norm 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:
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:
equipped with the norm
Multiplying the first equation of (1) by a test function and the second one by , integrating by parts over the terms involving and , 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 function and the second on by , integrating by parts over the terms involving , 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 have (weakly) continuous normal components on (see , Theorem 2.5), the mixed variational formulation of the coupled problem (1)-(5) can be stated as follows : Find that satisfies
where the bilinear forms and are defined on and , respectively, as:
By last, the linear forms L and G are defined as:
Theorem 2.1. If and , there exists a unique solution to the problem (17).
Proof. We will establish the existence and uniqueness of a weak solution by using the classical theory of mixed methods so-called Brezzi conditions (see, e.g., , Theorem and Corollary 4.1 in Chapter I).
· It is easy to prove that and are continuous. It is also clear that F and G are continuous and bounded. In summary, the triangular inequality and Cauchy-Schwarz inequality lead to:
· Now we define the null space of i.e. . From the classical Korn’s inequality ( , Section 5), it follows that there is a positive constant C such that:
Let . Because , then there exists unique such that . Thus, we obtain:
Hence, the coercivity of the bilinear form ,
· The second Brezzi condition that we need to verify is the inf-sup condition. Let , then there exists such that,
with the estimation
Thus, from (19) we have
Using the estimate (20), we obtain
Hence, the inf-sup condition,
Remark 2.1. Note that if g is of mean zero, (17) directly implies that (1), (2) 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  ), with the particularity that, in our case, for any , we have that .
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 . Indeed, we observe that taking the second equation of Darcy’s problem (2) we can write, for any ,
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 , it was (weakly) continuous normal components on , the variational form of the modified Stokes-Darcy problem can be written as follows: Find satisfying
where the bilinear forms and are defined on , , respectively, as:
By last, the linear forms and G are defined as:
Then, applying the classical theory of mixed methods it follows the well-posedness of the continuous formulation (27).
Theorem 2.2. There exists a unique solution to modified formulation (27). In addition, there exists a positive constant , depending on the continuous inf-sup condition constant for , the coercivity constant for and the boundedness constants for and , such that:
We end this section with some notation. In 2D, the curl of a scalar function w is given as usual by while in 3D, the curl of a vector function is given as usual by . Finally, let be the space of polynomials of total degree not larger than k. In order to avoid excessive use of constants, the abbreviations and stand for and , respectively, with positive constants independent of or .
3. 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 be a family of triangulations of with nondegenerate elements (i.e. triangles for and tetrahedrons for ). For any , we denote by the diameter of T and 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 such that , for all . We also assume that the triangulation is conform with respect to the partition of into and , namely each is either in or in (see Figures 3-5).
Let and be the corresponding induced triangulations of and . For any , we denote by (resp. ) the set of its edges ( ) or faces ( ) (resp. vertices) and set , . For we define
Notice that can be split up in the form
where . Note that is included in .
Figure 3. Isotropic element T in 2d.
Figure 4. Example of conforming mesh in 2d.
Figure 5. Example of nonconforming mesh in 2d.
With every edges , we associate a unit vector such that is orthogonal to E and equals to the unit exterior normal vector to if . For any and any piecewise continuous function , we denote by its jump across E in the direction of :
For , we set (see Figure 6 in 2D for illustration):
The triplet with is finite element [ , Page 83]. The local basis functions are defined by:
where for each , is barycentric coordonates of .
In classical reference element , the basis fonctions are given in by:
Figures 7-9 represent the surface of these functions in space.
The measure of an element or edge/face is denoted by and , respectively.
Based on the above notation, we introduce a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element space (larger than the space used in  )
and piecewise constant function space
where is the space of the restrictions to T of all polynomials of degree less than or equal to m. The space is equipped with the norm while
Figure 6. -nonconforming finite element T in 2d.
Figure 7. .
Figure 8. .
Figure 9. .
the norm on will be specified later on. The choice of is more natural than the one introduced in  since the space approximates only and not , while our a priori error analysis is only valid in this larger space.
Let us introduce the discrete divergence operator by
Then, we can introduce two bilinear forms
Then the finite element discretization of (27) is to find such that
This is the natural discretization of the modified weak formulation (27) except that the penalizing term is added. This bilinear form is defined by following the decomposition (37) of :
Here, is the length ( ) or diameter ( ) of E. Note that each element of only contributes with one jump term in .
Remark 3.1 The Equation (36) have the matrix representation
where (resp. ) denote the coefficients of (resp. ) expanded with respect to a basis for (resp. ).
We are now able to define the norm on (see  ):
In the sequel, we will denote by and 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 and :
Lemma 3.1. (Continuity of forms) There holds:
Lemma 3.2. (Coercivity of Ah) There is an such that:
Proof. Let . We have
We introduce the local space
and for , we define
with the semi-norm:
Using Young’s inequality and Green formula, we have:
· Estimate ( or ). We have by Cauchy-Schwarz inequality:
Also, we have:
Hence we deduce
· Now we estime the term . By Cauchy-Schwarz, we obtain:
Thus we deduce the estimation:
We apply Korn’s discrete inequality  and we get:
The estimates (49), (50), (51) and (52), lead to (42). The proof is complete. o
In order to verify the discrete inf-sup condition, we define the space:
We define also the Crouzeix-Raviart interpolation operator by:
Lemma 3.3. The operator is bounded: there is a constant depending on , and N such that
Proof. The proof is similar to ( , Lemma 4.5, Page 2695). o
Then, we have the following result
Lemma 3.4. (Inf-Sup condition) There exists a positive constant depending on , and N such that
Proof. We use Fortin argument , i.e. for each , we find such that:
Let . Then from ( , Corollary 2.4, Page 24), there exist vectoriel function satisfying
, hence . We take and we have:
Thus, we obtain
Using the system (58), we have:
From (59) et (60), we deduce:
The Inf-Sup condition holds and the proof is complete. o
From Lemma 3.2 and Lemma 3.4 we have the following result:
Theorem 3.1. There exists a unique solution to the problem (27).
3.2. A convergence Analysis
We now present an a priori analysis of the approximation error: The use of nonconforming finite element leads to , so the approximation error contains some extra consistency error terms. In fact, the abstract error estimates from  give the following result:
Lemma 3.5. Let be the solution of problem (27) and be the solution of the discrete problem (36). Then we have
where and are the consistency error terms define by:
Note that , thus .
For estiming the approximation error, we assume that the solution of problem (27) satisfies the smoothness assumptions:
1) , , ;
2) , , .
We begin with an estimate for the terms and .
Lemma 3.6. (Ref.  ) There hold:
Finally, let us consider the term . The smoothness assumption of implies , thus . Clearly,
Thus, we have
In order to evaluate the four face integrals, let us introduce two projections operators in the following.
For any and , denote by the constant space of the restrictions to E and the projection operator from on to such that
The operator has the property :
For any , we let be the function in such that
Using inequality (69), we obtain
Then we have the following lemma:
Lemma 3.7. (Estimation the four face integrals) There holds:
1) Estimate (71): We begin with an estimate for the first term . For any face , there exists at least one element such that . Then, from condition (68), Höder’s inequality and inequality (70), it follows that
2) Estimate (72):
We have , hence .
Furthermore, summing on faces, we obtain the estimate:
3) For the terms and , we use the same techniques as in the proof of the bounds for , , and we obtain:
The proof is complete. o
From Lemma 3.5, Lemma 3.6 and Lemma 3.7, now we derive the following convergence theorem:
Theorem 3.2 Let the solution of problem (27) satifies the smoothness assumption (Assumption 3.1). Let be the solution of the discrete problem (36). Then there exists a positive constant C depending on and such that:
4. 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++   in isotropic coupled mesh of Figure 10.
The solutions have been represented by Mathematica software. For simplicity we choose each domain , as the unit square, , and the permeability tensor is taken to be the identity. The interface , is the line , i.e. like the show Figure 11.
We consider the application on the open domain . In , we define and we obtain:
We choose quadratic pressure by
Figure 10. Isotropic mesh on coupled domain .
Figure 11. The domain Ω in 2d.
The exact solution satisfies the following condition:
and the Beavers-Joseph-Saffman interface conditions on [ ]:
Furthermore, we obtain the right-hand terms define by
Thus, in leads to
and in , is given by:
We study now the convergence of the speed and the pressure in each subdomain on a triangulation for the non-conforming finite elements of Crouzeix-Raviart. The results obtained are consistent, as shown in Figures 14-17. The order of convergence in norm approximately equal to 2 for the speed in each domain , and of order 1 for the pressure. Then, we represent on the same triangulation, the curves of isovalours of each component of the speed in , (see Figures 18-21).
Finally, we study the structure of the stiffness matrix in on a uniform triangulation. We find that, when the discretization step h becomes more and more infinitely small, the matrix becomes more and more sparse (cf. Table 1), which partly shows the effectiveness of the method digital used.
The exact solutions and p are represented by Figures 22-24, while the second members , , and are represented respectively by Figures 25-28.
Figure 12. Example of isotropic mesh in 2d.
Figure 13. Example of anisotropic mesh in 2d.
Figure 14. Error for the velocity in (log/log plot).
Figure 15. Error for the pressure in (log/log plot).
Figure 16. Error for the velocity in (log/log plot).
Figure 17. Error for the pressure in (log/log plot).
Figure 18. The isovalue of the first velocity component in .
Figure 19. The isovalue of the second velocity component in .
Figure 20. The isovalue of the first velocity component in .
Figure 21. The isovalue of the second velocity component in .
Table 1. Structure of rigidity Matrix on 10 iterations.
Figure 22. Component in .
· In this contribution, we have investigated a new mixed finite element method to solve the Stokes-Darcy fluid flow model without introducing any Lagrange multiplier. We have proposed a modification of the Darcy problem which allows us to apply a slight variant nonconforming Crouzeix-Raviart element to
Figure 23. Component in .
Figure 24. Pressure p in .
Figure 25. Right-hand term in .
Figure 26. Right-hand term in .
Figure 27. Right-hand term in .
Figure 28. Right-hand term in .
the whole coupled Stokes-Darcy problem. The proposed method is probably the cheapest method for Discontinuous Galerkin (DG) approximation of the coupled system, has optimal accuracy 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.
· Further works: Further it is well known that an internal layer appears at the interface as the permeability tensor degenerates, in that case anisotropic meshes have to be used in this layer (see for instance   ). Hence we intend to extend our results to such anisotropic meshes.
· bounded domain
· : the porous medium domain
· (resp. ) the unit outward normal vector along (resp. )
· : the fluid velocity
· p: the fluid pressure
· In 2D, the curl of a scalar function w is given as usual by
· In 3D, the curl of a vector function is given as usual by namely,
· : the space of polynomials of total degree not larger than k
· : triangulation of
· : the corresponding induced triangulation of ,
· For any , is the diameter of T and is the diameter of the largest ball inscribed into T
· : the set of all the edges or faces of the triangulation
· : the set of all the edges ( ) or faces ( ) of a element T
· : the set of all the vertices of a element T
· For ,
· For , we associate a unit vector such that is orthogonal to E and equals to the unit exterior normal vector to
· For , is the jump across E in the direction of
· In order to avoid excessive use of constants, the abbreviations and stand for and , respectively, with positive constants independent of or .
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. We are thankful to the editor and the anonymous reviewers for many valuable suggestions to improve this paper. This article has been posted since April 05, 2019 according to the link https://arxiv.org/abs/1908.01892.
There are no data underlying the findings in this paper to be shared.
This work has received no funding.
 Discacciati, M. and Quarteroni, A. (2009) Navier-Stokes/Darcy Coupling: Modeling, Analysis, and Numerical Approximation. Revista Matemática Complutense, 22, 315-426.
 Jäger, W. and Mikelić, A. (1996) On the Boundary Conditions of the Contact Interface between a Porous Medium and a Free Fluid. Annali della Scuola Normale Superiore. Classe di Scienze, 23, 403-465.
 Jäger, W., Mikelić, A. and Neuss, N. (2001) Asymptotic Analysis of the Laminar Visous Flow over a Porous Bed. SIAM Journal on Scientific Computing, 22, 2006-2028.
 Payne, L. and Straughan, B. (1998) Analysis of the Boundary Condition at the Interface between a Viscous Fluid and a Porous Medium and Related Modeling Questions. Journal de Mathématiques Pures et Appliquées, 77, 317-354.
 Houédanou, K.W. and Ahounou, B. (2016) A Posteriori Error Estimation for the Stokes-Darcy Coupled Problem on Anisotropic Discretization. Mathematical Methods in the Applied Sciences, 40, 3741-3774.
 Nicaise, S., Ahounou, B. and Houédanou, W. (2015) A Residual-Based Posteriori Error Estimates for a Nonconforming Finite Element Discretization of the Stokes-Darcy Coupled Problem: Isotropic Discretization. Afrika Matematika, 27, 701-729.
 Mu, M. and Xu, J. (2007) A Two-Grid Method of a Mixed Stokes-Darcy Model for Coupling Fluid Flow with Porous Media Flow. SIAM Journal on Numerical Analysis, 45, 1801-1813.
 Rui, H. and Zhang, R. (2009) A Unified Stabilized Mixed Finite Element Method for Coupling Stokes and Darcy Flows. Computer Methods in Applied Mechanics and Engineering, 198, 2692-2699.
 Arbogast, T. and Brunson, D. (2007) A Computational Method for Approximating A Darcy-Stokes System Governing a Vuggy Porous Medium. Computational Geosciences, 11, 207-218.
 Gatica, G.-N., Oyarzùa, R. and Sayas, F.-J. (2011) Convergence of a Family of Galerkin Discretizations for the Stokes-Darcy Coupled Problem. Numerical Methods for Partial Differential Equations, 27, 721-748.
 Chen, W. and Wang, Y. (2014) A Posteriori Error Estimate for H(div) Conforming Mixed Finite Element for the Coupled Darcy-Stokes System. Journal of Computational and Applied Mathematics, 255, 502-516.
 Babuška, I. and Gatica, G. (2010) A Residual-Based a Posteriori Error Estimator for the Stokes-Darcy Coupled Problem. SIAM Journal on Numerical Analysis, 48, 498-523.
 Gatica, G., Oyarzùa, R. and Sayas, F.-J. (2011) A Residual-Based a Posteriori Error Estimator for a Fully-Mixed Formulation of the Stokes-Darcy Coupled Problem. Computer Methods in Applied Mechanics and Engineering, 200, 1877-1891.
 Armentano, M.G. and Stockdale, M.L. (2018) A Unified Mixed Finite Element Approximations of the Stokes-Darcy Coupled Problem. Computers and Mathematics with Applications, 77, 2568-2584.
 Yu, J., Mahbub, M.A.A., Shi, F. and Zheng, H. (2018) Stabilized Finite Element Method for the Stationary Mixed Stokes-Darcy Problem. Advances in Difference Equations, 2018, Article No. 346.
 Li, R., Li, J., Chen, Z. and Gao, Y. (2015) A Stabilized Finite Element Method Based on Two Local Gauss Integrations for a Coupled Stokes-Darcy Problem. Journal of Computational and Applied Mathematics, 292, 92-104.
 Houédanou, K.W., Adetola, J. and Ahounou, B. (2017) Residual-Based a Posteriori Error Estimates for a Conforming Finite Element Discretization of the Navier-Stokes/Darcy Coupled Problem. Journal of Pure and Applied Mathematics: Advances and Applications, 18, 37-73.
 Gatica, G.N., Meddahi, S. and Oyarzùa, R. (2009) A Conforming Mixed Finite Element Method for the Coupling of Fluid Flow with Porous Media Flow. IMA Journal of Numerical Analysis, 29, 86-108.
 Girault, V. and Raviart, P.-A. (1986) Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms. Springer Series in Computational Mathematics, Volume 5, Springer, Berlin.
 Crouzeix, M. and Raviart, P. (1973) Conforming and Nonconforming Finite Element Methods for Solving the Stationary Stokes Equations I. ESAIM: Mathematical Modelling and Numerical Analysis, 7, 33-75.
 Houédanou, K.W. (2015) Analyse d’erreur a-posteriori pour quelques méthodes d’éléments finis mixtes pour le problème de transmission Stokes-Darcy: Discrétisations isotrope et anisotrope. Thèse de Doctorat, Université d’Abomey-Calavi, Godomey, 210 p.