TWO-GRID FINITE ELEMENT METHOD FOR THE STABILIZATION OF MIXED STOKES-DARCY MODEL

. A two-grid discretization for the stabilized ﬁnite element method for mixed Stokes-Darcy problem is proposed and analyzed. The lowest equal- order velocity-pressure pairs are used due to their simplicity and attractive computational properties, such as much simpler data structures and less com- puter memory for meshes and algebraic system, easier interpolations, and convenient usages of many existing preconditioners and fast solvers in simulations, which make these pairs a much popular choice in engineering practice; see e.g., [4, 27]. The decoupling methods are adopted for solving coupled systems based on the signiﬁcant features that decoupling methods can allow us to solve the submodel problems independently by using most appropriate numerical tech- niques and preconditioners, and also to reduce substantial coding tasks. The main idea in this paper is that, on the coarse grid, we solve a stabilized ﬁ- nite element scheme for coupled Stokes-Darcy problem; then on the ﬁne grid, we apply the coarse grid approximation to the interface conditions, and solve two independent subproblems: one is the stabilized ﬁnite element method for Stokes subproblem, and another one is the Darcy subproblem. Optimal error estimates are derived, and several numerical experiments are carried out to demonstrate the accuracy and eﬃciency of the two-grid stabilized ﬁnite ele- ment algorithm.

1. Introduction. In this paper, we consider the coupling problem of incompressible fluid flow with porous media flow. The incompressible fluid flow is described by the Stokes (or Navier-Stokes) equations, while the porous media flow is governed by the Darcy equations, and two flows are coupled through certain interface conditions. This coupled model has received wide applications over the past years, for example, the interaction between surface and subsurface flows [14,15,19,33], the subsurface flow in karst aquifers [8,9,25], industrial filtrations [21], and so on.
The two-grid method was firstly introduced by Mu and Xu in [39] to solve the mixed Stokes-Darcy problem. Later on, Cai and his coworkers proposed a decoupled and linearized two-grid algorithm to solve a coupled Navier-Stokes/Darcy model in [7]. However, the optimal error estimates for ∇u and p in L 2 norm are not theoretically proven, although their numerical tests suggested the optimal convergence results of order 2 by using stable finite element pair (P 2 , P 1 , P 2 ). To accord with the numerical results, Hou derived the optimal error estimates for fluid flow in [24] recently. In [49,47], some modified two-grid methods were developed for the mixed Stokes-Darcy model with optimal error estimates.
Most of the existing methods, based on the coupled/decoupled schemes for the Stokes-Darcy systems, require different suitable finite elements for fluid and porous media flow domains to meet the so-called inf-sup condition for stability. To overcome this problem, the lowest equal-order mixed finite element method accompanied with a well-developed stabilized method is chosen, since these finite element triples are computationally convenient in parallel processing and multigrid settings, for instance, the linear element pairs can ensure us to use much simpler data structures and less computer memory for meshes and algebraic system, easier interpolations, and also many existing preconditioners and fast solvers in simulations, see [2,41,35,48]. To solve the coupled steady Stokes-Darcy problem, Li et al presented and analyzed a stabilized finite element method in [34], by using the lowest equalorder finite element triples. In addition, the decoupling methods have attracted more and more attentions for solving coupled systems. As indicated by Mu and Xu in [39], decoupling methods allow us to solve the submodel problems independently by using most appropriate numerical techniques and preconditioners.
In this paper, by combining the two-grid technique and stabilized method, we develop a two-grid discretization based on the stabilized finite element method for the mixed Stokes-Darcy problem. On the coarse mesh, (P 1 , P 1 , P 1 ) triples are used to solve a stabilized coupling finite element approximation. On the fine mesh, two decoupled subproblems are solved independently: the Stokes problem is solved by a stabilized finite element solver (UMPACK solver) using (P 1 , P 1 ) pairs and the Darcy problem is solved by the standard finite element solver (Conjugate Gradient method). This natural combination of two-grid technique and stabilized methods retains the best features (such as, significant saving CPU time nearly without losing approximate accuracy, much simpler data structure for meshes and algebraic system) of both methods and overcomes many of their defects (larger algebraic system and much harder solving than decoupled Stokes and Darcy systems).
The rest of the paper is organized as follows. A mixed Stokes-Darcy model is described in the next section. In section 3, a two-grid discretization of the stabilized finite element method for the mixed Stokes-Darcy problem is given, while the optimal error estimates are derived. Then, several numerical experiments are presented in section 4. Finally, we draw some conclusions in section 5.

2.
Mixed Stokes-Darcy model. Consider the Stokes-Darcy model in a bounded domain Ω ⊂ R d (d = 2 or 3), which consists of a fluid region Ω f and a porous media The fluid flow in Ω f is governed by the Stokes equations: and the porous media flow in Ω p is described by the Darcy equations: Here, u and p represent the fluid velocity and the pressure, ν is the kinetic viscosity, f 1 is the external force on the fluid domain, f 2 is the source term, u p and ϕ denote the Darcy velocity and the piezometric head, n is the volumetric porosity, and K = {K ij } d×d stands for the hydraulic conductivity tensor which characterizes permeability of the rock. Throughout this paper, we always assume that the matrix K is symmetric positive definite. Combining Darcy's law (3) with the continuity equation (4), we obtain the following equation: On the interface Γ, the following three interface conditions are imposed: u · n f + u p · n p = 0 on Γ, p − νn f ∂u ∂n f = ρgϕ on Γ, where n f and n p denote the unit outward normal vectors on ∂Ω f and ∂Ω p , respetively. In particular, n f = −n p on the interface Γ. And also g denotes the gravitational acceleration, α is a positive parameter depending on the properties of the porous medium, and τ i (i = 1, ..., d − 1) denote the unit tangential vectors on the interface Γ. The first two interface conditions (6) and (7) state that the mass and normal force across the interface conserve respectively. The condition (8) is referred as the Beavers-Joseph-Saffman interface condition [44], which is a simplified Beavers-Joseph interface condition [1], and also approximately ensures that the slip velocity is proportional to the shear stress along the interface, see for instance [39].
For the sake of simplicity, we assume the homogeneous Dirichlet boundary conditions on Γ f and Γ p : These boundary conditions are not essential for the present algorithm and its numerical analysis. For any domain E ⊂ R d , we denote the L 2 (E) inner product and its norm by (·, ·) E and · 0,E , respectively. The norm of the Hilbert spaces H k (E) are denoted by · k,E , for scalar and vector-valued functions without confusion.
Define the following spaces: The space W is equipped with the following norm: ∀ u = (u, ϕ) ∈ W , We shall assume that n, ρ, g, ν and α are constants, and f 1 and f 2 are smooth enough. Then the weak formulation of the above Stokes-Darcy model reads as follows: find u = (u, ϕ) ∈ W and p ∈ Q, such that where The well-posedness of the mixed Stokes-Darcy model (12)-(13) can be found in [14] and [16]. Throughout this paper, we shall use the letter c (with or without subscripts) to denote a generic positive constant which may depend the physical parameters and is independent of the mesh parameters. It may stand for different values at its different occurrences.
3.1. Two-grid stabilized decoupled scheme. Let τ h (Ω) = {T } be a regular triangulation of Ω with mesh-size function h(x) whose value is the diameter h T of the element T containing x. We assume the triangulations τ h (Ω f ) and τ h (Ω p ) induced on the sub-domains Ω f and Ω p are compatible on the interface Γ. Let W h = H f,h × H p,h ⊂ W and Q h ⊂ Q be the finite element subspaces. Consider the lowest equal-order mixed finite element pairs for velocity and the pressure spaces, i.e., the finite element spaces H f,h , Q h and H p,h are selected as: Since the lowest equal-order finite element pairs do not satisfy the discrete inf-sup condition, some stabilization techniques are usually required. Thus, we introduce a stabilization term G(p h , q h ) (see [35,34,26] for details) and define it as following: where I is the identity operator, and Π : Then the finite element scheme for (12)-(13) based on (W h , Q h ) reads: The stability and error estimate for (14)- (15) can be found in [34]. For this system, one needs to solve a coupled problem, which usually results in difficulties in solving the resulting linear systems, such as larger scale of the algebraic system, and much hard solving than the decoupled Stokes and Darcy subsystems.
In [39], Mu and Xu proposed a two-grid scheme to decouple the mixed Stokes-Darcy problem, and obtained a convergence rate of O(h + H 3 2 ) for ∇u and p in L 2 norms for the finite element triples (P 1b −P 1 −P 1 ). Recently, in [24], Hou derived the optimal error estimate of order O(h + H 2 ). Here we present our modified two-grid stabilized method as follows: Algorithm 2 (Two-grid stabilized scheme): Step 1. Solve the coupled problem on a coarse grid: Step 2. On the fine grid, solve a modified problem: It's easy to verify that system (18)-(19) is equivalent to two decoupled subproblems: one is a Stokes problem on Ω f and another one is a Darcy problem on Ω p . In fact, on the fine grid, we solve the discrete Stokes problem in Ω f : find This can be viewed as a stabilization of mixed Stokes problem; see [2] for reference. And in the porous flow domain Ω p , we only need to solve a discrete Darcy problem: 3.2. Optimal error estimates. In this section, we present some error analysis for the proposed two-grid stabilized method for the Stokes-Darcy equations, and derive the optimal error estimates. Since the analysis of Algorithm 1 was well derived in [34], we carry out the present analysis through a comparison of Algorithm 2 and Algorithm 1.
Due to [34], by assuming the regularity u ∈ (H 2 (Ω f )) 2 , ϕ ∈ H 2 (Ω p ) and p ∈ H 2 (Ω f ), the following error estimates hold for the solution (u h , p h , ϕ h ) of Algorithm 1: and Consequently, on the coarse mesh, we have Theorem 3.1. Let (u h , p h , ϕ h ) and (u h , p h , ϕ h ) be the solutions of Algorithm 1 and Algorithm 2, respectively. Then we have the following error estimates Proof. Following the idea in [39], we compare Algorithm 1 and Algorithm 2 on the same fine grid h, and derive Let v h = (0, ψ h ) in (28), we have Then by choosing ψ h = ϕ h − ϕ h , we obtain Similar to [39], we introduce an auxiliary problem. Let θ ∈ H 1 (Ω f ) satisfy the following equations: Following [36], let H 1/2 00 (Γ) = [L 2 (Γ), H 1 0 (Γ)] 1/2 be the interpolation space. We have following inequalities: For the first term in the right hand side for the last term, thank to the definition G(p, q) and the properties of Π, we have

Combining all the above inequalities, we obtain
which yields the error estimate (26).
To get the optimal error estimate for the velocity and pressure in the fluid flow region, we introduce the following auxiliary problem in the porous media region: for any given It is easy to derive (28) and (29) yield For the right hand side of (36), we have Thanks to (26), we obtain the estimate for the first term in the right hand side of the above equality: It is also obvious that By (38) and (40), we have For (36) and (37), thank to [2], we obtain the following coercive property for the stabilization of Stokes problem: Hence the error estimate (27) is also proved.
By the triangle inequality, estimes (23)- (24) and Theorem 3.1, it is straightforward to get the following corollary.
In particular, when choosing H = O(h 1/2 ), we have 4. Numerical experiments. In this section, we will give fourth numerical examples to validate the stability and efficiency of the proposed two-grid stabilized method for the coupled Stokes-Darcy system. We use the lowest equal-order finite element triples (P 1 , P 1 , P 1 ) for the couple problem on coarse mesh, and lowest equal-order finite element pairs (P 1 , P 1 ) for the Stokes problem and P 1 -elements for the Darcy flow on fine mesh. For comparison, the well-known stable Mini elements (P 1b , P 1 ) are used for the Stokes problem and the linear Lagrangian elements P 1 are used for the Darcy flow. The relation of the coarse mesh H and the fine mesh h is chosen as h = H 2 . Moreover, the software package FreeFEM++ [23] is employed to implement the algorithms.
The force term and the boundary condition are determined by this exact solution.
In Tables 1, 2 and 3, by choosing the mesh scales H = h 1/2 , we list the errors between the exact solution and the finite element approximate solutions by SFEM (only computing on the fine mesh), TGM, and StbTGM, respectively. Note that, the error for p in L 2 -norm for SFEM is a little worse than O(h), but the errors for ∇ϕ, ∇u and p in the L 2 -norm by SFEM, TGM and StbTGM are all around the order of O(h). The computational time for the four schemes shows that the twogrid scheme can save significant amount of CPU time than SFEM and TGM, along with getting very similar approximation accuracy. Besides, by fixing the coarse mesh with H = 1/8, we execute Algorithm 2 (StbTGM) for different fine meshes with h = 1/16, 1/64 and 1/256. The approximate results are reported in Table  4. It can be observed that, along with the decreasing of fine mesh scale h, the approximate accuracy of the numerical solution may not always reduce. Therefore, suitable choices of the pairs of coarse and fine meshes are important in two-grid method.
In order to validate the effectiveness of our proposed stablization process, we solve the problem by two-grid algorithm in Mu and Xu [39] for finite element pair (P 1 , P 1 , P 1 ) (referred as TGM-(P 1 , P 1 , P 1 )), without using our stabilization method, and the results are displayed in Table 5. It is obvious that both the errors of ∇ϕ and ∇u are as good as the other methods, but the errors of p indicate the failure of this method which is a strong evidence that the stabilization is of significant importance for the mixed couple problem. Without stabilization, the two-grid scheme doesn't even work for fine mesh with h = 1/256. This is another reason why the stabilization technique is necessary. Example 2. In the second example, we focus on examining the effects of the physical parameters in our two-grid stabilized method. In this case, we take K = kI, and the same physical parameters ρ, g, n, ν, α are set to be 1, and choose the exact solution The force term and boundary condition are determined by the exact solution.
In Tables 6, 7 and 8, we list the approximation errors between the exact solution and the finite element approximate solutions by our StbTGM corresponding to k = 0.1, 0.01, and 0.001, respectively. The errors for ϕ and u increase along with the decreasing of parameter k, and the errors for p remain good accuracy. For these different values of k, our StbTGM can obtain optimal error orders of convergence, which validates the theoretical analysis of Corollary 1. and Ω p = (0, 1) × (0, 1) with interface Γ = (0, 1) × {1}. A modified driven cavity flow with the Dirichlet boundary conditions for the flow region is used as follows; See also [34], u = (sin(πx), 0), on (0, 1) × {1.25} u = (0, 0), on {0, 1} × (1, 1.25). The matrix K and other physical parameters ρ, g, n, ν, α are set to be 1, and f 1 = 0, f 2 = 0, and ϕ = 0 on Γ p are used.
The contour plots of the pressure obtained by three methods TGM, StbTGM and TGM-(P 1 , P 1 , P 1 ) without stabilization, with uniform meshes of H = 1/8 and  Figure 1. All these three methods capture smooth numerical pressure in the porous domain Ω p . Moreover, it can be easily observed that the numerical pressure in fluid domain are stable by TGM and StbTGM. However, without stabilization, oscillation appears for the pressure computed by TGM-(P 1 , P 1 , P 1 ) (Right one in the Figure 1). Furthermore, the velocity streamlines derived by the three methods are shown in Figure 2, from which on can conclude that these three methods can obtain nearly the same velocity. From all the three examples above, it can be concluded that our StbTGM is stable and efficient for the mixed Stokes-Darcy problem.   y), 0) at the inflow. For simplicity, the matrix K and other physical parameters ρ, g, n, ν, α are set to be 1, and f 1 = 0, f 2 = 0, and ϕ = 0 on Γ p are used. Since the interface between these two domains have two parts, we consider the following two situations as our testing cases: The velocity streamlines by StbTGM for two cases are depicted in Figure 3. Here the fluid and porous media regions are marked in yellow and green colors, respecitively. From this figure, the mass conservation is seen only across the interface (0, 1) × {0.5} in case 1 due to Dirichlet condition on {1} × (0, 0.5); while for case 2, the mass conservation can be found in both (0, 1) × {0.5} ∪ {1} × (0, 0.5). Our StbTGM captures these details very well. This test shows that our StbTGM is stable and valid for the mixed Stokes-Darcy problem.

Conclusion.
In this paper, we presented a two-grid discretization for the stabilized finite element method for the mixed Stokes-Darcy problem base on the lowest equal order finite element triples, and studied the optimal error estimates. Numerical tests showed that the present method is stable and efficient. In future work, this technique will be applied for studying the Navier-Stokes/Darcy coupled problems, and the related time-dependent coupled problems and so on. Besides, two grid methods [30,31,32] or decoupled preconditioning techniques [5] can be used to deal with large scale problems, particularly three dimensional problems with complex interfaces.