Convergence of the MAC scheme for the Stokes/Darcy coupling problem

全文

(1)J Sci Comput https://doi.org/10.1007/s10915-018-0660-7. Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem Ming-Cheng Shiue1 Ming-Chih Lai1. · Kian Chuan Ong1 ·. Received: 6 October 2017 / Revised: 22 January 2018 / Accepted: 30 January 2018 © Springer Science+Business Media, LLC, part of Springer Nature 2018. Abstract In this paper, we extend the MAC scheme for Stokes problem to the Stokes/Darcy coupling problem. The interface conditions between two separate regions are discretized and well-incorporated into the MAC grid setting. We first perform the stability analysis of the scheme for the velocity in both Stokes and Darcy regions and establish the stability for the pressure in both regions by considering an analogue of discrete divergence problem. Following the similar analysis on stability, we perform the error estimates for the velocity and the pressure in both regions. The theoretical results show the first-order convergence of the scheme in discrete L 2 norms for both velocity and the pressure in both regions. Moreover, in fluid region, the first-order convergence for the x-derivative of velocity component u and the y-derivative of velocity component v is also obtained in discrete L 2 norms. However, numerical tests show one order better for the velocity in Stokes region and the pressure in Darcy region. Keywords Stokes–Darcy flow · MAC scheme · Stability · Convergence · Finite difference method · Staggered grids. 1 Introduction The coupling of incompressible fluid flow with porous media flow has been an active research topic in recent years due to various applications of the filtration in biological and environmental engineering. The mathematical modeling of such physical processes consists of Stokes or. B. Ming-Cheng Shiue [email protected] Kian Chuan Ong [email protected] Ming-Chih Lai [email protected]. 1. Department of Applied Mathematics, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu 300, Taiwan. 123.

(2) J Sci Comput. Navier–Stokes description in the incompressible fluid region and Darcy’s law in the porous media region. These two flow regions are coupled at the fluid–porous interface through some physical interface conditions which we shall describe later. A detailed modeling, analysis and numerical approximation for the problem can be found in a recent review [13]. For past years, numerical methods for Stokes and Darcy coupling problems have been investigated mainly in the framework of finite element method such as in [1,5–10,12–15,19,20,23], just to name a few. Among those finite element discretizations, the fully coupled system can be either solved as a whole [5,9,10] or be decoupled into two separate subproblems with iteratively updating solution information across the interface [12,20]. There are other numerical approaches for the Stokes/Darcy coupling problem, such as spectral method (or pseudospectral method) [18,30], and boundary integral method [2] etc. Recently, an augmented immersed interface method based on finite difference scheme for Stokes/Darcy coupling problem with complex interface has been proposed in [21]. However, unlike most of the finite element methods, there is lack of convergence analysis of the method. This may be due to the absence of the variational formulations. In this paper, we propose a MAC (marker and cell) scheme for the Stokes/Darcy coupling problem and give the convergence proof for the scheme. To the best of our knowledge, this result is new. The MAC scheme proposed by Harlow and Welch [16] has been a popular finite difference scheme for Stokes and Navier–Stokes equations. The scheme adopts a nice grid layout in finite difference setting in which the velocity and pressure are located at different locations of a grid cell. More precisely, in 2D, the x component velocity u and the y component velocity v are defined at the middle points of vertical and horizontal edges, respectively; while the pressure p is defined at the cell center as depicted in Fig. 1. Although the MAC scheme has been developed in the 1960s, the first analysis and convergence for Stokes equations were carried out until 1992 by Nicolaides [24] simply because only limited mathematical tools are available for finite difference method. The author showed that the vorticity and the pressure are both first-order accurate. Later Han and Wu [17] proved that the MAC scheme can be obtained from a new mixed finite element method and showed that the first-order convergence for both the velocity (in the H 1 norm) and the pressure (in L 2 norm). Several similar convergent results using different finite element discretization can be found in the [25]. Until very recently, Li and Sun [22] proved the superconvergence (both velocity and pressure are second-order convergent in L 2 norm) for the MAC scheme on non-uniform grids using finite difference approach. Under the assumption of second-order convergence for the pressure, the authors were able to prove the second-order convergence of the velocity. However, the second-order convergence for the pressure is not exactly proved in the MAC framework. Rui and Li [25] established the inf-sup condition and the stability for both velocity and pressure; thus, the superconvergence for the MAC scheme can be obtained on nonuniform grids. For unstructured grids, we refer the interested reader to [3,4] for more details. For block-centered finite difference methods that is another type of MAC scheme, the secondorder convergence of the block-center scheme for incompressible and compressible DarcyForchheimer problems has proven in [27,28] and a two-grid block-centered finite difference scheme is also studied in [26]. In this paper, we first develop a finite difference discretization based on the MAC scheme for the Stokes/Darcy coupling equations. The interface conditions are discretized and can be well-incorporated into the MAC grid setting. Following the similar spirit used in [22,25] for Stokes equations, we conduct a stability and convergence analysis for the scheme on uniform grids. We would like to emphasize that the extension from Stokes problems to Stokes/Darcy coupling ones is not standard, especially in the proof of stability and convergence of the scheme. The major difficulty comes from the estimates of those relevant terms near the. 123.

(3) J Sci Comput. fluid–porous interface where three different interface conditions are imposed. Our analysis shows the first-order convergence in discrete L 2 norms for the velocity and the pressure in both incompressible fluid and porous regions. Moreover, in fluid region, the first-order convergence for the x-derivative of velocity component u and the y-derivative of velocity component v is also obtained in discrete L 2 norms. The rest of paper is organized as follows. In Sect. 2, we present the problem with the interface conditions. In Sect. 3, we present the MAC scheme for the Stokes/Darcy coupling equations and the discretization of the interface conditions. The major stability and error analysis are given in Sects. 4 and 5, respectively. Two numerical tests are given in Sect. 6 showing better convergence results than the theory. Concluding remarks are made in Sect. 7.. 2 The Stokes/Darcy Coupling Problem In this paper, the model under consideration consists of Stokes flow in the fluid region f and Darcy’s law in the porous media domain p , where these bounded domains f and p ⊂ R2 are assumed to be rectangular and separated by an interface as illustrated in Fig. 1. Let the boundary f ( p ) be ∂ f \ (∂ p \) respectively and n f (n p ) be the unit outward normal vector of the domain f ( p ) respectively. Let us denote u = (u, v) and p by the fluid velocity and pressure in f and u p and φ by the fluid velocity and pressure in p . In the region f , the Stokes flow (u, p) satisfies the following equations − νu + ∇ p = f 1 , in f ,. (2.1). ∇ · u = 0, in f ,. (2.2). u = 0, in f ,. (2.3). where ν is the viscosity and f 1 = ( f u , f v ) is the external force; in the region p , the Darcy’s flow (u p , φ) satisfies the following equations u p = −K ∇φ, in p ,. (2.4). ∇ · u p = f 2 , in p ,. (2.5). u p · n p = 0, in p ,. (2.6). where K is a symmetric and positive definite tensor representing the rock permeability divided by the fluid viscosity and f 2 is the external source. For simplicity, we choose K = κ I2 where κ is a positive constant and I2 is 2 × 2 identity matrix. By combining equations (2.4) and (2.5), we obtain − ∇ · (κ∇φ) = f 2 , in p . Here, the source f 2 is assumed to satisfy the solvability condition f 2 dx = 0, p. (2.7). (2.8). which is due to the no-slip (2.3) and no-flow boundary condition (2.6) on the boundaries f and p , respectively; and the mass conservation u p · n p + u · n f = 0 across the interface . In the present setting, this mass conservation across the interface results in (2.10). More detailedly, we have. 123.

(4) J Sci Comput. . f. ∇ · u dx +. p. ∇ · u p dx =. p. f 2 dx = 0.. The pressures ( p, φ) are assumed to satisfy the condition for the uniqueness of the solutions as follows: p dx + φ dx = 0. (2.9) f. p. To complete the problem (2.1)–(2.6), three conditions across the interface should be satisfied; namely, the mass conservation, the balance of normal forces, and the Beavers– Joseph–Saffman (BJS) condition, see the detailed physical meanings of those conditions in [13,20]. Readers who are interested in the well-posedness of the problem with above three interface conditions can refer to the review article [13]. Since the considered interface is a straight line in this paper, those three interface conditions can be simplified into ∂φ , ∂y ∂v p − φ = 2ν , ∂y ∂v k˜ ∂u u= + , α1 ∂ y ∂x v = −κ. where k˜ = νκ and α1 are positive constants. Here, the form of meaning of friction coefficient.. (2.10) (2.11) (2.12) ˜ 1 has the physical k/α. 3 Finite Difference Discretization Based on the MAC Scheme In this section, the finite difference discretization based on the MAC scheme for solving the problem (2.1)–(2.12) is presented. We start with the mesh description. Let = f ∪ p . For simplicity in presentation, the domain is assumed to be [0, L x ] × [−L y , L y ] where f = [0, L x ] × [0, L y ] , p = [0, L x ] × [−L y , 0] and L x , L y are positive constants. Given positive integers M and N , the mesh widths x and y are equal to L x /M and L y /N respectively. Let the nodal points (xi , y j ), 0 ≤ i ≤ M + 1, −N ≤ j ≤ N + 1 be defined as follows: 1 1 xi = i − x, y j = j − y. 2 2 For possible integers i, j, 0 ≤ i ≤ M, −N ≤ j ≤ N , we define xi+1/2 = (xi + xi+1 )/2 and y j+1/2 = (y j + y j+1 )/2. Here, the staggered grids are applied. Namely, the pressure is defined on one set of grid points while the velocities are defined on another set of grid points. We let u i+1/2, j , vi, j+1/2 and pi, j denote discrete approximations of the flow velocity u(xi+1/2 , y j ), v(xi , y j+1/2 ) and the pressure p(xi , y j ) respectively; let φi, j denote the discrete approximation of the pressure φ(xi , y j ) (see Fig. 1). To discretize (2.1), we employ central differences and derive u i+1/2, j+1 − 2u i+1/2, j + u i+1/2, j−1 u i+3/2, j − 2u i+1/2, j + u i−1/2, j + −ν ( x)2 ( y)2 pi+1, j − pi, j u + (3.1) = f i+1/2, j , 1 ≤ i ≤ M − 1, 1 ≤ j ≤ N , x. 123.

(5) J Sci Comput. y(j). (a). (b). j =N. vi,j+ 12 j+. Ωf. 1 2. ui− 12 ,j. j j−. 1 2. pi,j. ui+ 12 ,j. vi,j− 12. j=1. Γ. j=. 1 2. u. v. p. φ. j=0. Ωp. j = −N +1. i=1. i− 12 i i+ 21. i=M. x(i). Fig. 1 a Schematic representation of the finite difference discretization within the staggered grid framework. p and φ are defined at the cell centres, while u and v are defined at the centre of the cell faces. The interface is defined along j = 1/2. b Staggered arrangement of the variables. and . vi+1, j+1/2 − 2vi, j+1/2 + vi−1, j+1/2 vi, j+3/2 − 2vi, j+1/2 + vi, j−1/2, −ν + ( x)2 ( y)2 pi, j+1 − pi, j + = f i,v j+1/2 , 1 ≤ i ≤ M, 1 ≤ j ≤ N − 1, y. . (3.2). where u f i+1/2, j. 1 = x y. f i,v j+1/2 =. 1 x y. [xi ,xi+1 ]×[y j−1/2 ,y j+1/2 ]. . [xi−1/2 ,xi+1/2 ]×[y j ,y j+1 ]. f u (x, y) d xd y, f v (x, y) d xd y.. To discretize (2.2), we derive its discrete approximation at the mesh points (xi , y j ). That is, we obtain u i+1/2, j − u i−1/2, j vi, j+1/2 − vi, j−1/2 + = 0, 1 ≤ i ≤ M, 1 ≤ j ≤ N . x y. (3.3). 123.

(6) J Sci Comput. As for the finite difference discretization of the equation (2.7), similarly, we have φi, j+1 − 2φi, j − φi, j−1 φi+1, j − 2φi, j + φi−1, j −κ = ( f 2 )i, j , ( x)2 ( y)2 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0,. −κ. where ( f 2 )i, j. 1 = x y. (3.4). [xi−1/2 ,xi+1/2 ]×[y j−1/2 ,y j+1/2 ]. f 2 (x, y) d xd y.. The discrete approximations (3.1)–(3.4) are employed to determine the pressure and the velocity at the interior points in except the ones on the interface boundary . For certain of these equations, the information about the boundary conditions and the ghost points should be provided. For boundary conditions, we have u 1/2, j = u M+1/2, j = 0, j = 1, . . . , N , vi,N +1/2 = 0, i = 1, . . . , M, φ0, j = φ1, j , φ M+1, j = φ M, j j = −N + 1, . . . , 0, φi,−N = φi,−N +1 , i = 1, . . . , M.. (3.5) (3.6) (3.7) (3.8). As for the extrapolation of the ghost points on the boundary f ∪ p , we use linear extrapolation and the boundary condition which lead to u i+1/2,N +1 = −u i+1/2,N , i = 0, . . . , M,. (3.9). v0, j+1/2 = −v1, j+1/2 , j = 0, . . . , N ,. (3.10). v M+1, j+1/2 = −v M, j+1/2 , j = 0, . . . , N .. (3.11). Regarding the interface conditions on , we introduce the values of φi,1 , 1 ≤ i ≤ M (the pressure on the Darcy region) and u i+1/2,0 , 1 ≤ i ≤ M − 1 (the velocity on Stokes region) defined on the ghost points (xi , y1 ) and (xi+1/2 , y0 ) respectively. Furthermore, the interface conditions (2.10)–(2.12) can be approximated by choosing the values of φi,1 and u i+1/2,0 such that the following discretizations hold: φi,1 − φi,0 , i = 1, . . . , M, y vi,3/2 − vi,1/2 , i = 1, . . . , M, pi,1 − φi,0 = 2ν y u i+1/2,1 + u i+1/2,0 vi+1,1/2 − vi,1/2 k˜ u i+1/2,1 − u i+1/2,0 = + , 2 α1 y x i = 1, . . . , M − 1. vi,1/2 = −κ. (3.12) (3.13). (3.14). The idea of the discrete approximations (3.12) and (3.14) is to approximate the differential operators on the interface as follows: u i+1/2,1 − u i+1/2,0 ∂v vi+1,1/2 − vi,1/2 ∂φ φi,1 − φi,0 ∂u ≈ , ≈ , ≈ , ∂y y ∂y y ∂x x and apply the linear approximation to obtain the value of u on the interface as follows: u≈. 123. u i+1/2,1 + u i+1/2,0 . 2.

(7) J Sci Comput. We remark that the discrete approximation (3.13) of (2.11) is applied by using one-sided first order finite difference methods. The convergence for the unknowns (u, p) and φ in discrete L 2 norm is shown to be of first order in Sect. 4. However, the first order approximation on the interface condition does not contaminate the convergence of the second order for the velocity field which is illustrated by numerical experiments in Sect. 5. To discretize (2.9), we apply direct integration as follows: N M . M . pi, j x y +. 0 . φi, j x y = 0.. (3.15). i=1 j=−N +1. i=1 j=1. Remark that the following compatibility condition should be satisfied directly due to the definition of ( f 2 )i, j : M . 0 . ( f 2 )i, j x y = 0.. (3.16). i=1 j=−N +1. Now, we introduce the following standard forward and backward difference operators − Dx+ , Dx− , D + y and D y as follows: u i+3/2, j − u i+1/2, j u i+1/2, j − u i−1/2, j , Dx− u i+1/2, j = , x x u i+1/2, j+1 − u i+1/2, j u i+1/2, j − u i+1/2, j−1 , D− . = y u i+1/2, j = y y. Dx+ u i+1/2, j = D+ y u i+1/2, j. − Similarly, we can define the notations for Dx+ vi, j+1/2 , Dx− vi, j+1/2 , D + y vi, j+1/2 , D y vi, j+1/2 , + − + − Dx pi, j , Dx pi, j , D y pi, j and D y pi, j . These notations are also applied for Darcy’s flow. In summary, we rewrite the finite difference scheme for (3.1) to (3.4) as follows: − − u −ν Dx+ Dx− u i+1/2, j − ν D + y D y u i+1/2, j + Dx pi+1, j = f i+1/2, j ,. 1 ≤ i ≤ M − 1, 1 ≤ j ≤ N , −ν Dx+ Dx− vi, j+1/2. (3.17). − − ν D+ y D y vi, j+1/2. +. D− y pi, j+1. =. f i,v j+1/2 ,. 1 ≤ i ≤ M, 1 ≤ j ≤ N − 1,. (3.18). Dx− u i+1/2, j + D − y vi, j+1/2 = 0, 1 ≤ i ≤ M, 1 ≤ j ≤ N ,. (3.19). −κ Dx+ Dx− φi, j. (3.20). − − κ D+ y D y φi, j. = ( f 2 )i, j , 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0.. 4 Stability Analysis In this section, the stability analysis for the scheme (3.17)–(3.20) will be presented. Assume that the discrete solutions u i+1/2, j , vi, j+1/2 and φi, j satisfy the boundary conditions and interface conditions described in (3.5)–(3.14). We begin with introducing the following discrete norms that will be used later: u2 =. N M−1 . |u i+1/2, j |2 x y,. (4.1). i=1 j=1. v2 =. −1 M N i=1 j=1. |vi, j+1/2 |2 x y +. M i=1. |vi,1/2 |2. x y , 2. (4.2). 123.

(8) J Sci Comput. p2 =. N M . | pi, j |2 x y,. (4.3). i=1 j=1. φ2 =. M . 0 . |φi, j |2 x y, f 2 2 =. i=1 j=−N +1. f u 2 =. N M−1 . M . 0 . |( f 2 )i, j |2 x y,. (4.4). i=1 j=−N +1. u 2 | f i+1/2, j | x y,. (4.5). i=1 j=1. f v 2 =. −1 M N . | f i,v j+1/2 |2 x y +. i=1 j=1. Dx− u2 =. N M . M . v | f i,1/2 |2. i=1. x y , 2. (4.6). |Dx− u i+1/2, j |2 x y,. (4.7). i=1 j=1 2 D − y u =. −1 M−1 N. 2 |D − y u i+1/2, j+1 | x y. i=1 j=1. +. M−1 . 2 − 2 |D − y u i+1/2,N +1 | + |D y u i+1/2,1 |. x y 2. i=1. Dx− v2 =. −1 M N . ,. (4.8). |Dx− vi, j+1/2 |2 x y. i=2 j=1. +. N −1 . |Dx− v M+1, j+1/2 |2 + |Dx− v1, j+1/2 |2. x y. j=1. +. M . |Dx− vi,1/2 |2. i=2. 2. x y x y − + |Dx v1,1/2 |2 + |Dx− v M+1,1/2 |2 , 2 4 (4.9). 2 D − y v =. N M . 2 |D − y vi, j+1/2 | x y,. (4.10). i=1 j=1. Dx− φ2 =. M−1 . 0 . |Dx− φi+1, j |2 x y,. (4.11). i=1 j=−N +1 2 D − y φ =. M . 0 . i=1 j=−N +2. 2 |D − y φi, j | x y +. M i=1. 2 |D − y φi,1 |. x y . 2. (4.12). u v Notice that by the definition of f i+1/2, j ( f i, j+1/2 and ( f 2 )i, j ) respectively, we have the u u v v inequality f ≤ f L 2 ( f ≤ f L 2 and f 2 ≤ f 2 L 2 ) respectively where · L 2 denotes its corresponding L 2 norm. This implies that the discrete L 2 norms of the forcing terms are independent of the mesh widths. Before proceeding the stability analysis, we need the following discrete Poincare inequalities:. 123.

(9) J Sci Comput. Lemma 4.1 (4.13) u2 ≤ L 2x Dx− u2 , y 2 v2 ≤ L y + L y D − (4.14) y v , 2.

(10) 1/2 0 M . − 2 − 2 2 2. ( f 2 )i, j φi, j x y ≤ 2 L x + L y f 2 Dx φ + D y φ .. i=1 j=−N +1. (4.15) Proof To show (4.13), we have u2 =. N M−1 . |. i . Dx− u m+1/2, j x|2 x y. i=1 j=1 m=1. ≤. N M−1 M−1 i=1 j=1. ≤. |Dx− u m+1/2, j |2 M( x)3 y. m=1. L 2x Dx− u2 .. (4.16). As for (4.14), by applying similar techniques, we have v2 =. −1 M N N i=1 j=1. ≤. i=1 j=1. ≤. Ly +. 2 x y +. M N . m= j+1. −1 M N N . . D− y vi,m+1/2 y. i=1. . 2 |D − N ( y)3 x + y vi,m+1/2 |. m=2. y 2. . D− y vi,m+1/2 y. 2. m=1. x y 2. L y y 2 D − y v 2. 2 L y D − y v .. (4.17). As for (4.15), for any i ≥ i and j ≥ j we have ( f 2 )i, j φi, j − ( f 2 )i, j φi , j. j i . ≤ |( f 2 )i, j | |Dx− φm, j | x + |D − φ | y . y i ,m m=i +1. (4.18). m=0. Multiplying (4.18) by x y and summing all i and j, applying the condition (3.16) and using Cauchy–Schwarz inequality, (4.15) is obtained. Then, the proof of Lemma 4.1 is complete.

(11). Now, in order to make our presentation of the stability analysis clear, we need the following lemmas. Lemma 4.2 − ν Dx+ Dx− u i+1/2, j − ν Dx+ D − y vi, j+1/2 = 0, 1 ≤ i ≤ M − 1, 1 ≤ j ≤ N , (4.19) − + − −ν D + y Dx u i+1/2, j − ν D y D y vi, j+1/2 = 0, 1 ≤ i ≤ M, 1 ≤ j ≤ N − 1.. (4.20). The proof of Lemma 4.2 is established by directly applying (3.19).. 123.

(12) J Sci Comput. Lemma 4.3 N M−1 . u i+1/2, j Dx− pi+1, j x y +. i=1 j=1. =−. vi, j+1/2 D − y pi, j+1 x y. i=1 j=1. M . vi,1/2 pi,1 x,. i=1. −κ. −1 M N . M . 0 . . (4.21). − D φ Dx+ Dx− φi, j + D + y y i, j φi, j x y. i=1 j=−N +1. = κDx− φ2 + κ. M . 0 . 2 (D − y φi, j ) x y − κ. i=1 j=−N +2. M . D− y φi,1 φi,0 x. (4.22). i=1. Proof To show (4.21), we have N M−1 . N M−1 . u i+1/2, j Dx− pi+1, j x y =. i=1 j=1. u i+1/2, j ( pi+1, j − pi, j ) y. i=1 j=1. =. N M−1 j=1. =. =−. M−1 . i=1. N M j=1. u i+1/2, j pi+1, j −. i=1. u i−1/2, j pi, j −. i=2. M N . u i+1/2, j pi, j y. M−1 . u i+1/2, j pi, j y. i=1. Dx− u i+1/2, j pi, j x y,. (4.23). j=1 i=1. where the boundary conditions (3.5) are applied. By applying the same technique, we have −1 M N . vi, j+1/2 D − y pi, j+1 x y = −. i=1 j=1. N M . D− y vi, j+1/2 pi, j x y −. i=1 j=1. M . vi,1/2 pi,1 x,. i=1. (4.24) where the boundary conditions (3.6) are applied. (4.21) can be derived by summing (4.23) and (4.24) and applying (3.19). Similar procedures can be applied to obtain (4.22).

(13). Lemma 4.4 −ν. N M−1 . u i+1/2, j Dx+ Dx− u i+1/2, j x y = νDx− u2 ,. (4.25). i=1 j=1. −ν. N M−1 . − − 2 u i+1/2, j D + y D y u i+1/2, j x y = νD y u. i=1 j=1. +ν. M−1 i=1. 123. u i+1/2,1 + u i+1/2,0 − D y u i+1/2,1 x, 2. (4.26).

(14) J Sci Comput −1 M N . −ν. vi, j+1/2 Dx+ Dx− vi, j+1/2 x y = ν. −1 M−1 N. i=1 j=1. +. −ν. i=1 j=1. N −1 . ν 2. |Dx− vi, j+1/2 |2 x y. . |Dx− v M+1, j+1/2 |2 + |Dx− v1, j+1/2 |2 x y,. (4.27). j=1 −1 M N . − − 2 vi, j+1/2 D + y D y vi, j+1/2 x y = νD y v. i=1 j=1. +ν. M . vi,1/2 D − y vi,3/2 x,. (4.28). i=1. −ν. N M−1 . u i+1/2, j Dx+ D − y vi, j+1/2 x y. i=1 j=1. =ν. −1 M−1 N. − D− y u i+1/2, j+1 Dx vi+1, j+1/2 x y. i=1 j=1. +ν. M−1 . u i+1/2,1 Dx− vi+1,1/2 x,. (4.29). i=1. −ν. −1 M N . − vi, j+1/2 D + y Dx u i+1/2, j x y. i=1 j=1. =ν. −1 M−1 N. − D− y u i+1/2, j+1 Dx vi+1, j+1/2 x y.. (4.30). i=1 j=1. Proof To show (4.25), we observe that −ν. N M−1 . u i+1/2, j Dx+ Dx− u i+1/2, j x y. i=1 j=1. = −ν. N M−1 . u i+1/2, j Dx+ u i+1/2, j. − u i+1/2, j Dx− u i+1/2, j. y. i=1 j=1. = −ν. N M−1 j=1. = −ν. = −ν. u i−1/2, j Dx− u i+1/2, j. u i+1/2, j Dx− u i+1/2, j. y. −. M−1 . u i+1/2, j Dx− u i+1/2, j. y. i=1. −|Dx− u i+1/2, j |2 x y + u M+1/2, j Dx− u M+1/2, j y. i=1. −u 1/2, j Dx− u 3/2, j y = νDx− u2 ,. M−1 i=1. i=2. N M j=1. −. i=1. N M j=1. u i+1/2, j Dx+ u i+1/2, j. . (4.31). 123.

(15) J Sci Comput. where boundary conditions (3.9) are applied. Then (4.25) is derived. As for (4.26), we have −ν. M−1 N . − u i+1/2, j D + y D y u i+1/2, j x y. i=1 j=1. = −ν. M−1 N i=1. = −ν. j=1. u i+1/2, j−1 D − y u i+1/2, j −. 2 − −|D − y u i+1/2, j | x y + u i+1/2,N D y u i+1/2,N +1 x. −u i+1/2,1 D − y u i+1/2,1 x. . 2 |D − y u i+1/2, j | x y +. i=1 j=2. +ν. M−1 i=1. u i+1/2, j D − y u i+1/2, j x. j=1. j=2. N M−1 . u i+1/2, j D − y u i+1/2, j x. N . j=2. M−1 N i=1. =ν. N j=1. M−1 +1 N i=1. = −ν. u i+1/2, j D + y u i+1/2, j −. M−1 ν 2 − 2 |D − x y u | + |D u | i+1/2,N +1 i+1/2,1 y y 2 i=1. u i+1/2,1 + u i+1/2,0 − D y u i+1/2,1 x 2. 2 = νD − y u + ν. M−1 i=1. u i+1/2,1 + u i+1/2,0 − D y u i+1/2,1 x, 2. (4.32). where the linear extrapolation conditions (3.9) are used. Then (4.26) is obtained. (4.27) and (4.28) can be derived by using similar techniques. The proofs are omitted. To prove (4.29), we observe − + Dx+ D − y vi, j+1/2 = D y Dx vi, j−1/2 .. (4.33). Notice that −ν. N M−1 . + u i+1/2, j D − y Dx vi, j−1/2 x y. i=1 j=1. = −ν. M−1 N i=1. = −ν. j=1. i=1. u i+1/2, j−1 Dx+ vi, j−1/2 −. j=2. M−1 N . u i+1/2, j Dx+ vi, j−1/2 x. N . u i+1/2, j Dx+ vi, j−1/2 x. j=1 + + −D − y u i+1/2, j Dx vi, j−1/2 x y + u i+1/2,N Dx vi,N +1/2 x. j=2. −u i+1/2,1 Dx+ vi,1/2 x. 123. N j=1. M−1 +1 N i=1. = −ν. u i+1/2, j Dx+ vi, j+1/2 −. .

(16) J Sci Comput. =ν. N M−1 . + D− y u i+1/2, j Dx vi, j−1/2 x y + ν. i=1 j=2. M−1 . u i+1/2,1 Dx+ vi,1/2 x. (4.34). i=1. Then, (4.29) is proven. Again, (4.30) can be done by using the same technique which proof is omitted. Finally, the proof of Lemma 4.4 is complete..

(17). Now, we are in a position to state and show the following theorem for the stability analysis of the scheme (3.1)–(3.14) as follows: Theorem 4.1 Given the mesh widths x and y satisfying. νκ 2α22 y ≤ min , , 2L y L y where α2 =. √ 2 κ˜ α1 ,. (4.35). we have. ν ν ν ν u2 + v2 v2 + Dx− u2 + D − 2 2 2L x 8L y + 2νκ 2 8 y −1 M−1 N. +ν. − D− y u i+1/2, j+1 + Dx vi+1, j+1/2. 2. x y. i=1 j=1 M−1 . +ν. 2 |D − y u i+1/2,N +1 |. i=1. +. x y 2. N −1. ν − |Dx v1, j+1/2 |2 + |Dx− v M+1, j+1/2 |2 x y 2 j=1. κ κ 2 + Dx− φ2 + D − y φ 2 2 L 2y + L2 ≤ C f := x f u 2 + 2ν 2ν. νκ 4. f v 2 +. L 2x + L 2y κ. f 2 2 .. (4.36). Proof By multiplying (3.17), (3.18), (3.20), (4.19) and (4.20) by u i+1/2, j x y, vi, j+1/2 x y, φi, j x y, u i+1/2, j x y and vi, j+1/2 x y respectively, summing the resulting equations for all i and j, and applying Lemmas (4.3) and (4.4), we obtain 2νDx− u2 + ν. −1 M−1 N. − D− y u i+1/2, j+1 + Dx vi+1, j+1/2. 2. x y. i=1 j=1. +ν. M−1 . 2 − 2 |D − y u i+1/2,N +1 | + |D y u i+1/2,1 |. i=1. +. x y 2. N −1. ν − |Dx v1, j+1/2 |2 + |Dx− v M+1, j+1/2 |2 x y 2 j=1. 2 − 2 − 2 +2νD − y v + κDx φ + κD y φ M − = φ φ − 2νv D v vi,1/2 pi,1 + κ D − i,1/2 y i,3/2 x y i,1 i,0 i=1. 123.

(18) J Sci Comput. −ν. M−1 i=1. M−1 u i+1/2,1 + u i+1/2,0 − D y u i+1/2,1 x − ν u i+1/2,1 Dx− vi+1,1/2 x 2 i=1. M . +κ. 2 |D − y φi,1 |. i=1. +. N M−1 . x y 2. u u i+1/2, j f i+1/2, j x y +. i=1 j=1. +. M . 0 . −1 M N . vi, j+1/2 f i,v j+1/2 x y. i=1 j=1. ( f 2 )i, j φi, j x y. i=1 j=−N +1. := I1 + I2 + I3 + I4 + I5 .. (4.37). Now, we estimate the terms Ii , i = 1, . . . , 5. For the term I1 , by applying the interface conditions (3.12) and (3.13), we obtain I1 =. M . vi,1/2 ( pi,1 − φi,0 ) − 2νvi,1/2 D − v y i,3/2 x = 0.. (4.38). i=1. To estimate the term I2 , we have I2 = −ν. M−1 i=1. −ν. u i+1/2,1 + u i+1/2,0 − D y u i+1/2,1 + Dx− vi+1,1/2 x 2. M−1 i=1. u i+1/2,1 − u i+1/2,0 − Dx vi+1,1/2 x 2. 2 M−1 M−1 u i+1/2,1 − u i+1/2,0 να1 u i+1/2,1 + u i+1/2,0 x − ν = −√ Dx− vi+1,1/2 x, 2 2 κ˜ i=1 i=1 (4.39) where the interface conditions (3.14) are applied. To estimate the second term on the righthand side of (4.39), we have −( y − α2 ) α2 y u i+1/2,1 + D − vi+1,1/2 , y + α2 y + α2 x √ ˜ where the constant α2 is equal to (2 κ)/(α 1 ). This leads to u i+1/2,0 =. u i+1/2,1 − u i+1/2,0 y α2 y = D − vi+1,1/2 . u i+1/2,1 − 2 y + α2 2( y + α2 ) x Plugging (4.41) into the second term on the right-hand side of (4.39), we obtain −ν. M−1 i=1. = −ν. u i+1/2,1 − u i+1/2,0 − Dx vi+1,1/2 x 2 M−1 i=1. 123. y α2 y u i+1/2,1 − Dx− vi+1,1/2 Dx− vi+1,1/2 x y + α2 2( y + α2 ). (4.40). (4.41).

(19) J Sci Comput. = −ν. M−1 . M−1 y α2 y |D − vi+1,1/2 |2 x u i+1/2,1 Dx− vi+1,1/2 x + ν y + α2 2( y + α2 ) x. i=1. i=1. = J1 + J2 .. (4.42). To estimate the term J1 , we have J1 = −ν. M−1 y u i+1/2,1 vi+1,1/2 − vi,1/2 y + α2 i=1. y = −ν y + α2 =ν. M. u i−1/2,1 vi,1/2 −. M−1 . i=2. u i+1/2,1 vi,1/2. i=1. M y − Dx u i+1/2,1 vi,1/2 x y + α2 i=1. =ν. M N y − Dx u i+1/2,1 − D− v y x i, j+1/2 y y + α2 i=1. j=1. = (due to Young’s inequality) N 2 M y L y α2 − ≤ν x y |Dx− u i+1/2,1 |2 + D y vi, j+1/2 y + α2 4α2 Ly i=1. ≤ν. L y y 4α2 ( y + α2 ). j=1. M . |Dx− u i+1/2,1 |2 x y + ν. i=1. α2 2 D − y v y + α2. L y y 2 ≤ν Dx− u2 + νD − y v . 4α22. (4.43). To estimate the term J2 , we have J2 = ν. M−1 i=1. =ν. α2 y |D − vi+1,1/2 |2 x 2( y + α2 ) x. 2 M−1 2 u i+1/2,1 + u i+1/2,0 α2 y u x − D− i+1/2,1 y 2( y + α2 ) α2 2 i=1. M−1 . 4 u i+1/2,1 + u i+1/2,0 2 2 + |D − y u i+1/2,1 | 2 2 α 2 i=1 4 u i+1/2,1 + u i+1/2,0 − − D y u i+1/2,1 x α2 2 M−1 2 y 2 u i+1/2,1 + u i+1/2,0 2 x ≤ν 1+. 2 ( y + α2 ) α2. =ν. α2 y 2( y + α2 ). i=1. +ν. . α2 1+ y + α2 2. M−1 i=1. 2 |D − y u i+1/2,1 |. x y , 2. (4.44). where is a positive parameter which will be determined later.. 123.

(20) J Sci Comput. Combining (4.39), (4.43) and (4.44), we obtain M−1 2ν 2 y u i+1/2,1 + u i+1/2,0 2 α2 − x I2 ≤ − α2 ( y + α2 ). 2 i=1. M−1 . − α2 x y 1+ |D y u i+1/2,1 |2 +ν y + α2 2 2 i=1. L y y 2 Dx− u2 + νD − +ν y v . 4α22 Now, set =. 2 y α2 ,. I2 ≤ ν. (4.45). we can infer from (4.45) that M−1 . 2 |D − y u i+1/2,1 |. i=1. L y y x y 2 Dx− u2 + νD − +ν y v . 2 4α22. (4.46). To estimate the term I3 , we have I3 =. 2 M M N x y x y 1 − 1 ≤ |vi,1/2 |2 |D y vi, j+1/2 | y κ 2 κ 2 i=1. i=1. j=1. L y y 2 D − (4.47) y v . 2κ To estimate the term I4 , by applying discrete Poincare inequalities and Young’s inequality, we have ≤. I4 ≤ u f u + v f v . y v D − ≤ Ly Ly + y v f 2 2 L y L y + y 2 ν ν L ≤ Dx− u2 + D − v2 + x f u 2 + (4.48) f v 2 . 2 2 y 2ν 2ν To estimate the term I5 , by applying Lemma 4.1, we have

(21) 2 I5 ≤ 2 L 2x + L 2y f 2 Dx− φ2 + D − y φ 2 + L2 L x y κ 2 Dx− φ2 + D − + ≤ (4.49) f 2 2 . y φ 2 κ Combining the estimates of Ii , i = 1, . . . , 5, we can infer from (4.37) that −1 M−1 2 N L y y 3 − 2 − D− u + ν x y − νD x y u i+1/2, j+1 + Dx vi+1, j+1/2 2 2 4α2 i=1 j=1 L x Dx− u f u +. +ν. M−1 i=1. +. x y 2. N −1 ν − |Dx v1, j+1/2 |2 + |Dx− v M+1, j+1/2 |2 x y 2 j=1. 123. 2 |D − y u i+1/2,N +1 |. .

(22) J Sci Comput. L y y ν κ κ 2 − 2 − 2 + − D − y v + Dx φ + D y φ 2 2κ 2 2 y L 2x + L 2y L 2x u 2 L y L y + 2 ≤ f + f v 2 + f 2 2 . 2ν 2ν κ By taking y to satisfy the following conditions. (4.50). L y y L y y 3 ν ν ≥ 1, and − − ≥ , 2 2 2κ 4 4α22 which implies. y ≤ min. νκ 2α22 , 2L y L y. (4.51). ,. (4.52)

(23). the proof of Theorem 4.1 is complete by using Lemma 4.1. To show the boundedness of the pressure p and φ, we need the following lemmas: Lemma 4.5 Let the mesh widths satisfy (4.35) and y = x, we have M−1 u i+1/2,1 + u i+1/2,0 2 Cf α2 8α2 L 2x x ≤ K 1 = +2 Cf + . 2 ν Ly ν i=1. Proof To begin with, we rewrite the summation for Dx− v, apply the triangle inequality and have 2 M−1 M−1 N y − 2 − − − D y vi+1, j+1/2 + D y vi, j+1/2 (Dx vi+1,1/2 ) x y = x y x i=1. i=1. j=1. 4L y y 2 D − ≤ y v . (N y = L y ) ( x)2. (4.53). Then, we use (4.41) and have y. M−1 . 2 (D − y u i+1/2,1 ) x y. i=1. ≤. M−1 M−1 2α22 y − 8 y 2 u x y + (Dx vi+1,1/2 )2 x y. i+1/2,1 ( y + α2 )2 ( y + α2 )2 i=1. i=1. 8 y 2 ≤ 2 u2 + 8L y D − y v , α2. (4.54). where the condition y = x and the inequality (4.53) are applied. To show the desired inequality, we rewrite (4.37) in Theorem 4.1 as follows: 2νDx− u2 + ν. −1 M−1 N. − D− y u i+1/2, j+1 + Dx vi+1, j+1/2. 2. x y. i=1 j=1. +ν. M−1 i=1. 2 |D − y u i+1/2,N +1 |. N −1. ν − x y + |Dx v1, j+1/2 |2 + |Dx− v M+1, j+1/2 |2 x y 2 2 j=1. 123.

(24) J Sci Comput 2 − 2 − 2 +2νD − y v + κDx φ + κD y φ M − vi,1/2 pi,1 + κ D − = φ φ − 2νv D v i,1/2 y i,3/2 x y i,1 i,0 i=1 M−1 . −ν. M−1 u i+1/2,1 + u i+1/2,0 − D y u i+1/2,1 x − ν u i+1/2,1 Dx− vi+1,1/2 x 2. i=1. −. ν 2. +κ. i=1. M−1 . 2 |D − y u i+1/2,1 | x y. i=1 M . 2 |D − y φi,1 |. i=1. +. N M−1 . x y 2. u u i+1/2, j f i+1/2, j x y +. i=1 j=1. + :=. M . −1 M N . vi, j+1/2 f i,v j+1/2 x y. i=1 j=1. 0 . ( f 2 )i, j φi, j x y i=1 j=−N +1 I1 + I2 + I3 + I4 + I5 .. (4.55). The estimates for the terms Ii , i = 1, 3, 4, 5 remain the same as in (4.38), (4.47), (4.48) and (4.49). To estimate the term I2 , we have I2. M−1 M−1 2ν u i+1/2,1 + u i+1/2,0 2 ν − =− x − D y u i+1/2,1 Dx− vi+1,1/2 x y α2 2 2 i=1. − =−. ν 2. M−1 . i=1. 2 |D − y u i+1/2,1 | x y. i=1. M−1 2ν u i+1/2,1 + u i+1/2,0 2 x α2 2 i=1. ν − α2 ≤−. ν α2. M−1 i=1. M−1 i=1. u i+1/2,1 + u i+1/2,0 2 u i+1/2,1 + u i+1/2,0 2. . D− y u i+1/2,1 x y. 2 x +. M−1 2 ν y − D y u i+1/2,1 x y, 4α2 i=1. (4.56) where Cauchy–Schwarz and Young’s inequalities are applied from the second line to the third line. Combining the estimates (4.38), (4.47), (4.48), (4.49) and (4.56) and applying (4.35), we obtain M−1 M−1 2 ν u i+1/2,1 + u i+1/2,0 2 ν y − D y u i+1/2,1 x y. x ≤ C f + (4.57) α2 2 4α2 i=1. i=1. The desired inequality is obtained due to (4.57), (4.54), (4.35) and Theorem 4.1.. 123.

(25).

(26) J Sci Comput. In addition to Lemma 4.5, we need the following lemma extended from the discrete divergence problem for finite difference scheme for Stokes equations to the present Stokes/Darcy coupling equations on a staggered grid. Lemma 4.6 For any given pi, j , 1 ≤ i ≤ M, 1 ≤ j ≤ N and φi, j , 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0 satisfying (3.15), there exist two vectors u˜ i+1/2, j , 0 ≤ i ≤ M, 0 ≤ j ≤ N + 1 and v˜i, j+1/2 , 0 ≤ i ≤ M + 1, 0 ≤ j ≤ N satisfying the following properties: u˜ 1/2, j = u˜ M+1/2, j = 0, 0 ≤ j ≤ N + 1,. (4.58). u˜ i+1/2,N +1 = −u˜ i+1/2,N , 0 ≤ i ≤ M,. Dx− u˜ i+1/2, j. +. (4.59). v˜i,N +1/2 = 0, 0 ≤ i ≤ M + 1,. (4.60). v˜0, j+1/2 = −v˜1, j+1/2 , v˜ M+1, j+1/2 = −v˜ M, j+1/2 , 0 ≤ j ≤ N ,. (4.61). D− y v˜i, j+1/2. = pi, j , 1 ≤ i ≤ M, 1 ≤ j ≤ N ,. (4.62). and there exists a vector φ˜ i, j , 0 ≤ i ≤ M + 1, −N ≤ j ≤ 1 satisfying Dx− φ˜ 1, j = Dx− φ˜ M+1, j = 0, −N + 1 ≤ j ≤ 0, −˜ ˜ D− y φi,−N +1 = 0, −κ D y φi,1 = v˜i,1/2 , 1 ≤ i ≤ M, −κ Dx+ Dx− φ˜ i, j. −˜ − κ D+ y D y φi, j. (4.63) (4.64). = φi, j , 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0.. (4.65). Moreover, we have Dx− u ˜ 2 + D − ˜ 2 + Dx− v ˜ 2 + D − ˜ 2 y u y v ˜ 2 + D − ˜ 2 +u ˜ 2 + v ˜ 2 + Dx− φ y φ +. ≤ C˜ d p2 + φ2 ,. M−1 i=1. u˜ i+1/2,1 + u˜ i+1/2,0 2. 2 x (4.66). where C˜ d = (1 + 2y + (1 + κ12 )(L 2x + L 2y ))Cd is a constant independent of the mesh widths x and y and Cd is a constant defined in Lemma 4.7. L. To prove Lemma 4.6, we need the following lemma related to the finite difference scheme for the Stokes problem on the whole region with homogenous Dirichlet boundary conditions: Lemma 4.7 For any given pi, j , 1 ≤ i ≤ M, 1 ≤ j ≤ N and φi, j , 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0 satisfying (3.15), there exist a positive constant Cd independent of the mesh widths x and y and two vectors Ui+1/2, j , 0 ≤ i ≤ M, −N ≤ j ≤ N + 1 and Vi, j+1/2 , 0 ≤ i ≤ M + 1, −N + 1 ≤ j ≤ N satisfying U1/2, j = U M+1/2, j = 0, −N ≤ j ≤ N + 1,. (4.67). Ui+1/2,N +1 = −Ui+1/2,N , Ui+1/2,−N = −Ui+1/2,−N +1 , 0 ≤ i ≤ M,. (4.68). Vi,N +1/2 = Vi,−N +1/2 = 0, 0 ≤ i ≤ M + 1,. (4.69). V0, j+1/2 = −V1, j+1/2 , VM+1, j+1/2 = −VM, j+1/2 , −N + 1 ≤ j ≤ N ,. (4.70). Dx− Ui+1/2, j + D − y Vi, j+1/2 = pi, j , 1 ≤ i ≤ M, 1 ≤ j ≤ N ,. (4.71). Dx− Ui+1/2, j + D − y Vi, j+1/2 = φi, j , 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0, − 2 − 2 2 2 Dx U + D y U 2 + Dx− V 2 + D − y V ≤ C d ( p + φ ),. (4.72) (4.73). 123.

(27) J Sci Comput. where · is the corresponding discrete L 2 norm in the region and Cd depends on the size of the domain . The proof of Lemma 4.7 can be obtained by following the same processes in [29]. Thus, the proof is omitted. Now the proof of Lemma 4.6 is presented as follows: Proof of Lemma 4.6 Assume that two vectors Ui+1/2, j , 0 ≤ i ≤ M, −N ≤ j ≤ N + 1 and Vi, j+1/2 , 0 ≤ i ≤ M + 1, −N + 1 ≤ j ≤ N are defined in Lemma 4.7 satisfying (4.67)–(4.73). For the Stokes region, the velocity field (u, ˜ v) ˜ is defined as follows: u˜ i+1/2, j = Ui+1/2, j , 0 ≤ i ≤ M, 0 ≤ j ≤ N + 1, v˜i, j+1/2 = Vi, j+1/2 , 0 ≤ i ≤ M + 1, 0 ≤ j ≤ N . It is easy to check that (4.58)–(4.62) are obtained from the definition of (u, ˜ v) ˜ and Lemma 4.7. For the Darcy region, thanks to (4.72), the existence of a vector φ˜ satisfying (4.63)–(4.65) is equivalent to show that the vector φ˜ satisfies Dx− φ˜ 1, j = Dx− φ˜ M+1, j = 0, −N + 1 ≤ j ≤ 0, −˜ ˜ D− y φi,−N +1 = 0, −κ D y φi,1 = v˜i,1/2 , 1 ≤ i ≤ M, Ui+1/2, j =. −κ Dx− φ˜ i+1, j ,. Vi, j+1/2 =. (4.74) (4.75). ˜ −κ D − y φi, j+1 ,. 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0.. (4.76). From (4.74) and (4.76), φ˜ is determined by the choice of φ˜ 1, j , j = −N + 1, . . . , 0. From (4.75) and (4.76), φ˜ is determined by the choice of φ˜ i,0 , i = 1, . . . , M. Therefore, this vector φ˜ exists up to a constant. Note that, the velocity field in Darcy region are defined by (Ui+1/2, j , Vi, j+1/2 ) = (−κ Dx− φi+1, j , −κ D − y φi, j+1 ). To show (4.66), we use the fact M−1 i=1. =. u˜ i+1/2,1 + u˜ i+1/2,0 2. M−1 i=1. −. N +1 j=2. 2 x. D− y u˜ i+1/2, j. +. 0 . D− y Ui+1/2, j. j=−N +1. 2. ( y)2 x 4. Ly 2 D − ≤ y U . 2. (4.77). By applying (4.77)and discrete Poincare inequality, we have ˜ 2 + D − ˜ 2 + Dx− v ˜ 2 + D − ˜ 2 Dx− u y u y v ˜ 2 + D − ˜ 2 +u ˜ 2 + v ˜ 2 + Dx− φ y φ +. M−1 i=1. u˜ i+1/2,1 + u˜ i+1/2,0 2. Ly 2 − 2 − 2 + 1+ ≤ D − y U + Dx V + D y V 2 . 1 + 1 + 2 U 2 + V 2 κ Dx− U 2. 123. 2 x.

(28) J Sci Comput. Ly 1 2 L x + L 2y ≤ 1+ + 1+ 2 Dx− U 2 2 κ 2 − 2 − 2 +D − y U + Dx V + D y V .. (4.78)

(29). Then, (4.66) is derived by using Lemma 4.7 and (4.78).. Now, we are ready to state and show the boundedness of the pressure p and φ as follows: Theorem 4.2 Let the mesh widths satisfy (4.35) and y = x, we have p2 + φ2. να22 2ν ˜ 2ν u 2 v 2 ≤ C p = max 2ν, κ, 1, K1 + K1 + f + f . Cd 16C f + α2 α2 Ly . Proof By multiplying (3.17), (3.18), (4.65), (4.19) and (4.20) by u˜ i+1/2, j x y , v˜i, j+1/2 x y, φi, j x y, u˜ i+1/2, j x y and v˜i, j+1/2 x y respectively, summing for all i and j and adding the resulting equations, we obtain φ2 + p2 = 2ν. N M . Dx− u˜ i+1/2, j Dx− u i+1/2, j x y. i=1 j=1. +ν. −1 M−1 N. − D− y u˜ i+1/2, j+1 +Dx v˜i+1, j+1/2. . − D− y u i+1/2, j+1 +Dx vi+1, j+1/2. x y. i=1 j=1. M−1 ν − − − − D y u˜ i+1/2,N +1 D y u i+1/2,N +1 + D y u˜ i+1/2,1 D y u i+1/2,1 x y + 2 i=1. +ν. M−1 i=1. +. M−1 u˜ i+1/2,1 + u˜ i+1/2,0 − D y u i+1/2,1 x + ν u˜ i+1/2,1 Dx− vi+1,1/2 x 2. N −1 . ν 2. i=1. Dx− v˜ M+1, j+1/2 Dx− v M+1, j+1/2 + Dx− v˜1, j+1/2 Dx− v1, j+1/2 x y. j=1. +2ν. N M . − D− y v˜i, j+1/2 D y vi, j+1/2 x y +. i=1 j=1. +κ. M . Dx− φi, j Dx− φ˜ i, j x y + κ. i=2 j=−N +1. −. v˜i,1/2 φi,0 − pi,1 + 2ν D − v x i,3/2 y. i=1. 0 . N M−1 . M . 0 . −˜ D− y φi, j D y φi, j x y. i=1 j=−N +2. u u˜ i+1/2, j f i+1/2, j x y −. i=1 j=1. M . −1 M N . v˜i, j+1/2 f i,v j+1/2 x y.. (4.79). i=1 j=1. By using the fact ν. M−1 i=1. +. M−1 u˜ i+1/2,1 + u˜ i+1/2,0 − u˜ i+1/2,1 Dx− vi+1,1/2 x D y u i+1/2,1 x + ν 2 i=1. ν 2. M−1 . − D− y u˜ i+1/2,1 D y u i+1/2,1 x y. i=1. 123.

(30) J Sci Comput. M−1 u i+1/2,1 + u i+1/2,0 2ν u˜ i+1/2,1 + u˜ i+1/2,0 x = α2 2 2 i=1. +. ν 2. M−1 . − − D− y u˜ i+1/2,1 Dx vi+1,1/2 + D y u i+1/2,1 x y,. i=1. and applying Cauchy–Schwarz inequality, Theorem 4.1, Lemmas 4.1 and 4.5, it is inferred from (4.79) that ε1 p2 + φ2 ≤ 2ν Dx− u ˜ 2 + D − ˜ 2 + Dx− v ˜ 2 + D − ˜ 2 y u y v 2 2 2 ˜ 2 + κD − ˜ 2 +u ˜ + v ˜ + κDx− φ y φ. M−1 να 2 2ν u˜ i+1/2,1 + u˜ i+1/2,0 2 1 2ν 16C f + x + K1 + 2 K1 α2 2 2ε1 α2 Ly i=1 + f u 2 + f v 2 , (4.80) +. where ε1 is a positive parameter which will be determined later. Then, by applying Lemma 4.6, we can deduce from (4.80) that ε1 2ν ˜ p2 + φ2 ≤ max 2ν, κ, 1, Cd ( p2 + φ2 ) 2 α2 να 2 2ν 1 16C f + K 1 + 2 K 1 + f u 2 + f v 2 . + 2ε1 α2 Ly. (4.81). By choosing ε1 to satisfy ε1 =. 1 , ˜ max 2ν, κ, 1, 2ν α2 C d

(31). the proof of Theorem 4.2 is complete.. 5 Error Analysis In this section, the error estimates for the scheme (3.1)–(3.14) will be derived. We begin with introducing the following notations for the errors: u ei+1/2, j = u(x i+1/2 , y j ) − u i+1/2, j ,. ei,v j+1/2 = v(xi , y j+1/2 ) − vi, j+1/2 , φ. p. ei, j = p(xi , y j ) − pi, j , ei, j = φ(xi , y j ) − φi, j . Thus, we can derive the following equations for the errors: u + − u − u −ν Dx+ Dx− ei+1/2, j − ν D y D y ei+1/2, j + Dx ei+1, j = Ri+1/2, j , p. 1 ≤ i ≤ M − 1, 1 ≤ j ≤ N , −ν Dx+ Dx− ei,v j+1/2. − v − ν D+ y D y ei, j+1/2. 1 ≤ i ≤ M, 1 ≤ j ≤ N − 1,. 123. (5.1) +. p D− y ei, j+1. =. Ri,v j+1/2 , (5.2).

(32) J Sci Comput u − v d Dx− ei+1/2, j + D y ei, j+1/2 = Ri, j , 1 ≤ i ≤ M, 1 ≤ j ≤ N , φ −κ Dx+ Dx− ei, j. − φ − κ D+ y D y ei, j. =. φ Ri, j , 1. (5.3). ≤ i ≤ M, −N + 1 ≤ j ≤ 0,. (5.4). φ. u d v where the terms Ri+1/2, j , Ri, j+1/2 , Ri, j and Ri, j are defined as u + − + − − Ri+1/2, j = −ν Dx Dx u(x i+1/2 , y j ) − ν D y D y u(x i+1/2 , y j ) + Dx p(x i+1 , y j ) u − f i+1/2, j , 1 ≤ i ≤ M − 1, 1 ≤ j ≤ N ,. (5.5) (5.6). Ri,d j =. − − −ν Dx+ Dx− v(xi , y j+1/2 ) − ν D + y D y v(x i , y j+1/2 ) + D y p(x i , y j+1 ) − f i,v j+1/2 , 1 ≤ i ≤ M, 1 ≤ j ≤ N − 1, Dx− u(xi+1/2 , y j ) + D − y v(x i , y j+1/2 ), 1 ≤ i ≤ M, 1 ≤ j ≤ N ,. φ Ri, j. −κ Dx+ Dx− φ(xi ,. Ri,v j+1/2. =. =. − y j ) − κ D+ y D y φ(x i ,. (5.7). y j ) − ( f 2 )i, j ,. 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0.. (5.8). For boundary conditions, we have u u v e1/2, j = e M+1/2, j = 0, 1 ≤ j ≤ N ; ei,N +1/2 = 0, 1 ≤ i ≤ M.. (5.9). As for the extrapolation of the ghost points on the boundary f ∪ p , we define u(xi+1/2 , y N +1 ) = −u(xi+1/2 , y N ), 0 ≤ i ≤ M, v(x0 , y j+1/2 ) = −v(x1 , y j+1/2 ), v(x M+1 , y j+1/2 ) = −v(x M , y j+1/2 ), 0 ≤ j ≤ N , φ(x0 , y j ) = φ(x1 , y j ), φ(x M+1 , y j ) = φ(x M , y j ), −N + 1 ≤ j ≤ 0, ∂φ φ(xi , y−N ) = φ(xi , y−N +1 ), φ(xi , y1 ) = φ(xi , y0 ) + y (xi , y1/2 ), ∂y 1 ≤ i ≤ M, which imply u u ei+1/2,N +1 = −ei+1/2,N , 0 ≤ i ≤ M, v e0, j+1/2 φ e0, j φ ei,−N. = = =. v v −e1, j+1/2 , e M+1, j1/2 φ φ φ e1, j , e M+1, j = e M, j , φ ei,−N +1 , 1 ≤ i ≤ M.. =. (5.10) −evM, j+1/2 ,. 0 ≤ j ≤ N,. (5.11). −N + 1 ≤ j ≤ 0,. (5.12) (5.13). Regarding the interface conditions on , we have φ. v = −κ D − ei,1/2 y ei,1 , 1 ≤ i ≤ M, p φ ei,1 − ei,0 u u + ei+1/2,0 ei+1/2,1. 2. = =. (5.14). p v 2ν D − y ei,3/2 + (R )i,1/2 ,. k˜. 1 ≤ i ≤ M, u − v u D− e + D e y i+1/2,1 x i+1,1/2 + (R )i+1/2,1/2 ,. α1 1 ≤ i ≤ M − 1.. (5.15). (5.16). p. where (R )i,1/2 and (Ru )i+1/2,1/2 are defined as (R )i,1/2 = p(xi , y1 ) − φ(xi , y0 ) − 2ν D − y v(x i , y3/2 ), 1 ≤ i ≤ M, p. (Ru )i+1/2,1/2. (5.17). u(xi+1/2 , y1 ) + u(xi+1/2 , y0 ) = 2. 123.

(33) J Sci Comput. −. k˜. D− y u(x i+1/2 , y1 ) +. α1 1 ≤ i ≤ M − 1, and due to that. u ei+1/2,0. . Dx− v(xi+1 , y1/2 ). , (5.18). is evaluated at the ghost point, it is defined as follows. u ei+1/2,0 = u(xi+1/2 , y0 ) − u i+1/2,0 ,. u(xi+1/2 , y0 ) =. (5.19). −( y − α2 ) α2 y u(xi+1/2 , y1 ) + D − v(xi+1 , y1/2 ). y + α2 y + α2 x. (5.20). About the condition of the uniqueness, we have N M . p. ei, j x y +. M . 0 . φ. ei, j x y = R p ,. (5.21). i=1 j=−N +1. i=1 j=1. where Rp =. N M . p(xi , y j ) x y +. M . 0 . φ(xi , y j ) x y.. (5.22). i=1 j=−N +1. i=1 j=1. Notice that the definitions of the discrete norms of the errors eu , ev , e p and eφ are the same as u, v, p and φ respectively. Before establishing the error estimates, the following lemmas are needed in the sequel. Lemma 5.1 eu 2 ≤ L 2x Dx− eu 2 , y v 2 v 2 L y D − e ≤ L y + ye . 2. (5.23) (5.24). The proof of Lemma 5.1 is similar to the proof in Lemma 4.1. Thus, the proof is omitted. Lemma 5.2 u + − v + d −ν Dx+ Dx− ei+1/2, j − ν D y D y ei, j+1/2 = −ν Dx Ri, j ,. 1 ≤ i ≤ M − 1, 1 ≤ j ≤ N , u −ν Dx+ Dx− ei+1/2, j. (5.25). − v − ν D+ y D y ei, j+1/2. =. d −ν D − y Ri, j+1 ,. 1 ≤ i ≤ M, 1 ≤ j ≤ N − 1.. (5.26). The proof of Lemma 5.2 is established by directly applying (5.3). Lemma 5.3 N M−1 . u − ei+1/2, j Dx ei+1, j x y + p. i=1 j=1. =−. M . v ei,1/2 ei,1 x − p. p. M . N M . p. Ri,d j ei, j x y,. i=1 j=1 0 . i=1 j=−N +1. 123. ei,v j+1/2 D − y ei, j+1 x y. i=1 j=1. i=1. −κ. −1 M N . . φ φ − φ Dx+ Dx− ei, j + D + D e y y i, j ei, j x y,. (5.27).

(34) J Sci Comput. =κ. M−1 . . 0 . φ. Dx− ei+1, j. 2. x y + κ. i=1 j=−N +1. −κ. φ φ M ei,1 − ei,0 i=1. y. M . . 0 . φ. D− y ei, j. 2. x y. i=1 j=−N +2 φ. ei,0 x.. (5.28). Lemma 5.4 −ν. N M−1 . u + − u − u 2 ei+1/2, j Dx Dx ei+1/2, j x y = νDx e ,. (5.29). i=1 j=1. −ν. N M−1 . u + − u ei+1/2, j D y D y ei+1/2, j x y. i=1 j=1 u 2 = νD − y e +ν. u M−1 ei+1/2,1 i=1. −ν. −1 M N . u + ei+1/2,0. 2. u D− y ei+1/2,1 x,. (5.30). ei,v j+1/2 Dx+ Dx− ei,v j+1/2 x y. i=1 j=1. =ν. −1 M−1 N. v 2 |Dx− ei+1, j+1/2 | x y. i=1 j=1. +. N −1 ν v 2 | |Dx− evM+1, j+1/2 |2 + |Dx− e1, x y, j+1/2 2. (5.31). j=1. −ν. −1 M N . − v − v 2 ei,v j+1/2 D + y D y ei, j+1/2 x y = νD y e + ν. i=1 j=1. M . v v ei,1/2 D− y ei,3/2 x,. i=1. (5.32) −ν. N M−1 . u + − v ei+1/2, j Dx D y ei, j+1/2 x y. i=1 j=1. =ν. −1 M−1 N. u − v D− y ei+1/2, j+1 Dx ei+1, j+1/2 x y + ν. i=1 j=1. M−1 . u v ei+1/2,1 Dx− ei+1,1/2 x,. i=1. (5.33) −ν. −1 M N . − u ei,v j+1/2 D + y Dx ei+1/2, j x y. i=1 j=1. =ν. −1 M−1 N. u − v D− y ei+1/2, j+1 Dx ei+1, j+1/2 x y.. (5.34). i=1 j=1. The proofs of Lemma 5.3 and 5.4 are similar to the proofs in Lemma 4.3 and 4.4. Thus, the proofs are omitted, too. For the convenience of the notation, we denote the maximum norm of the r -th derivatives r of any function u as |Dr u|∞ = | ∂ m∂x∂un y |∞ where r = m + n and m and n are nonnegative integers.. 123.

(35) J Sci Comput. To perform the error estimates of the finite difference scheme, we need the following p u d u v lemma for the estimates of the terms Ri+1/2, j , Ri, j+1/2 , Ri, j , (R )i,1/2 , (R )i+1/2,1/2 and p R : Lemma 5.5 u 2 4 3 2 4 Ri+1/2, j = O(( x) )(|D u|∞ + |D p|∞ ) + O(( y) )|D u|∞ + O( x y), j = 1, N u Ri+1/2,1 = O(( x)2 )(|D 4 u|∞ + |D 3 p|∞ ) + O( y)|D 3 u|∞ + O( x y), ν u = u yy (xi+1/2 , y N ) + O(( x)2 )(|D 4 u|∞ + |D 3 p|∞ ) Ri+1/2,N 4 +O( y)|D 3 u|∞ + O( x y),. Ri,v j+1/2 = O(( x)2 )|D 4 v|∞ + O(( y)2 )(|D 4 u|∞ + |D 3 p|∞ ) + O( x y), i = 1, M ν Ri,v j+1/2 = vx x (xi , y j+1/2 ) + O( x)|D 3 v|∞ + O(( y)2 )(|D 4 u|∞ + |D 3 p|∞ ) 4 +O( x y), i = 1, M Ri,d j = O(( x)2 )|D 3 u|∞ + O(( y)2 )|D 3 v|∞ , Dx+ Ri,d j = O(( x)2 )|D 4 u|∞ + O(( x)2 + ( y)2 )|D 4 v|∞ , d 2 2 4 2 4 D− y Ri, j+1 = O(( x) + ( y) )|D u|∞ + O(( y) )|D v|∞ , φ. Ri, j = (O(( x)2 ) + O(( y)2 ))(|D 3 φ|∞ + |D 2 f 2 |∞ ), i = 1, M and j = 0, −N + 1, φ. Ri, j = (O( x) + O(( y)2 ))(|D 3 φ|∞ + |D 2 f 2 |∞ ), i = 1, M and j = 0, −N + 1, φ. Ri, j = (O(( x)2 ) + O( y))(|D 3 φ|∞ + |D 2 f 2 |∞ ), i = 1, M and j = 0, −N + 1, φ. Ri, j = (O( x) + O( y))(|D 3 φ|∞ + |D 2 f 2 |∞ ), i = 1, M and j = 0, −N + 1, p. (R )i,1/2 = O( y)(|Dp|∞ + |Dφ|∞ + |D 2 v|∞ ), (Ru )i+1/2,1/2 = O(( x)2 )|D 3 v|∞ + O(( y)2 )|D 3 u|∞ , R p = (O(( x)2 ) + O(( y)2 ))(|D 2 φ|∞ + |D 2 p|∞ ). The proof of Lemma 5.5 is based on Taylor’s expansion and the assumption on the regularity of the solutions and the forcing. The proof of Lemma 5.5 is skipped. Remark that in the fluid region, the truncation errors inside the domain are of second order but the one near the boundary is some O(1) term plus first-order terms; in the porous region, the truncation errors inside the domain are of second order but the one near the boundary is of the first order. Near the interface, the truncation errors are both of the first order. The following lemma is also needed to control the boundedness of the error estimates of the pressure and related boundary terms: Lemma 5.6 Let the mesh widths satisfy (4.35) and y = x, we have u M−1 ei+1/2,1 i=1 p 2. u + ei+1/2,0. 2. 2. x ≤ K 2 = 2 K 1 + L x |u|2∞ ,. e + e ≤ K 3 = 2 C p + L x L y | p|2∞ + L x L y |φ|2∞ , φ 2. where K i , i = 2, 3 are constants depending only on the size of the domain, the forcing and the solutions.. 123.

(36) J Sci Comput. The proof of Lemma 5.6 can be proven by using the triangle inequality, Lemma 4.5 and Theorem 4.2. Now, we state and prove our main theorem of the error estimates for the scheme (3.1)– (3.14) as follows: Theorem 5.1 Let the mesh widths satisfy (4.35) and y = x, we have ν ν ν ν ev 2 + Dx− eu 2 + D − eu 2 + e v 2 2L 2x 8L 2y + 2νκ 2 16 y +ν. −1 M−1 N. u − v D− y ei+1/2, j+1 + Dx ei+1, j+1/2. 2. x y. i=1 j=1. +. M−1 x y ν − u |D y ei+1/2,N +1 |2 2 2 i=1. +. ν 2. N −1 . v 2 − v 2 |Dx− e1, j+1/2 | + |Dx e M+1, j+1/2 |. x y 2. j=1. φ 2 +κDx− eφ 2 + κD − ye . = O ( x)2 + ( y)2 .. (5.35). u v Proof By multiplying (5.1), (5.2), (5.4), (5.25) and (5.26) by ei+1/2, j x y, ei, j+1/2 x y, φ. u v ei, j x y, ei+1/2, j x y and ei, j+1/2 x y respectively, summing for all i and j and adding all the resulting equations, we obtain. 2νDx− eu 2 + ν. −1 M−1 N. u − v D− y ei+1/2, j+1 + Dx ei+1, j+1/2. 2. x y. i=1 j=1. +ν. M−1 . u 2 − u 2 |D − y ei+1/2,N +1 | + |D y ei+1/2,1 |. x y 2. i=1. +. N −1 ν − v |Dx e1, j+1/2 |2 + |Dx− evM+1, j+1/2 |2 x y 2 j=1. v 2 − φ 2 − φ 2 +2νD − y e + κDx e + κD y e M p φ v v v v ei,1/2 = ei,1 − ei,1/2 ei,0 − 2νei,1/2 D− e y i,3/2 x i=1. −ν. u M−1 ei+1/2,1. u + ei+1/2,0. 2. i=1. +κ. M i=1. +. φ. 2 |D − y ei,1 |. N M−1 i=1 j=1. u D− y ei+1/2,1 x − ν. M−1 . u v ei+1/2,1 Dx− ei+1,1/2 x. i=1. x y 2. u u ei+1/2, j Ri+1/2, j x y +. −1 M N . ei,v j+1/2 Ri,v j+1/2 x y. i=1 j=1. 123.

(37) J Sci Comput. +. M . 0 . φ. φ. Ri, j ei, j x y +. i=1 j=−N +1. −ν. N M−1 . N M . p. Ri,d j ei, j x y. i=1 j=1. u Dx+ Ri,d j ei+1/2, j x y − ν. i=1 j=1. −1 M N . d v D− y Ri, j+1 ei, j+1/2 x y. i=1 j=1. := I6 + I7 + I8 + I9 + I10 + I11 .. (5.36). Now, we estimate the terms Ii , i = 6, . . . , 11. For the term I6 , by applying the interface conditions (5.14) and (5.15), we obtain I6 =. M . v ei,1/2. p φ v − v ei,1 − ei,0 − 2νei,1/2 D y ei,3/2 x. i=1. =. M . v ei,1/2 (R )i,1/2 x p. i=1. =. N M . v −D − y ei, j+1/2 y(R )i,1/2 x p. i=1 j=1. ≤. M N i=1. v |D − y ei, j+1/2 | y. 2 1/2 M. j=1. v 2 ≤ ε2 νD − ye +. 1/2 p. |(R )i,1/2 |2. x. i=1 M Ly p |(R )i,1/2 |2 x, 4ε2 ν. (5.37). i=1. where ε2 is a positive parameter which is determined later. To estimate the term I7 , we have I7 = −ν. u M−1 ei+1/2,1 i=1. −ν. u M−1 ei+1/2,1 i=1. =−. u + ei+1/2,0. 2 u − ei+1/2,0. 2. u − v (D − y ei+1/2,1 + Dx ei+1,1/2 ) x. v x Dx− ei+1,1/2. 2 u M−1 u 2ν ei+1/2,1 + ei+1/2,0 x α2 2 i=1. +. 2ν α2. −ν. u M−1 ei+1/2,1. 2. i=1. u M−1 ei+1/2,1 i=1. u + ei+1/2,0. u − ei+1/2,0. 2. (Ru )i+1/2,1/2 x. v x, Dx− ei+1,1/2. (5.38). where the interface conditions (5.16) are applied. To estimate the second term on the righthand side of (5.38), by applying Cauchy–Schwarz inequality and Lemma 5.6, we have u M−1 u 2ν ei+1/2,1 + ei+1/2,0 u (R )i+1/2,1/2 x α2 2 i=1. 123.

(38) J Sci Comput. 2ν ≤ α2 =. M−1 . u u + ei+1/2,0 ei+1/2,1. 2. i=1. 1/2 M−1 2ν K 2. α2. 2 x. 1/2 M−1 i=1. 1/2 |(Ru )i+1/2,1/2 |2 x. 1/2 |(Ru )i+1/2,1/2 |2 x. .. (5.39). i=1. To estimate the third term on the right-hand side of (5.38), we observe u = ei+1/2,0. −( y − α2 ) u α2 y 2 y v ei+1/2,1 + Dx− ei+1,1/2 + (R u )i+1/2,1/2 . y + α2 y + α2 y + α2 (5.40). This leads to u u − ei+1/2,0 ei+1/2,1. 2. =. α2 y y eu − D − ev y + α2 i+1/2,1 2( y + α2 ) x i+1,1/2 y − (R u )i+1/2,1/2 . y + α2 . (5.41). Plugging (5.41) into the third term on the right-hand side of (5.38), we obtain −ν. u M−1 ei+1/2,1. 2. i=1. = −ν. M−1 i=1. −. y y + α2. = −ν. M−1 i=1. +ν. M−1 i=1. +ν. u − ei+1/2,0. M−1 i=1. v x Dx− ei+1,1/2. y α2 y eu − D − ev y + α2 i+1/2,1 2( y + α2 ) x i+1,1/2 v (Ru )i+1/2,1/2 Dx− ei+1,1/2 x. y eu D − ev x y + α2 i+1/2,1 x i+1,1/2 α2 y D − ev x D − ev 2( y + α2 ) x i+1,1/2 x i+1,1/2 y v (R u )i+1/2,1/2 Dx− ei+1,1/2 x y + α2 . = J3 + J4 + J5 .. (5.42). To estimate the terms J3 and J4 , we apply the similar procedure in estimating J1 and J2 respectively in Theorem 4.1 and have L y y v 2 Dx− eu 2 + νD − ye , 4α22 2 M−1 u u 2 y 2 ei+1/2,1 + ei+1/2,0 x J4 ≤ ν 1+ ε3 2 ( y + α2 ) α2 J3 ≤ ν. (5.43). i=1. M−1 ε3 − u x y α2 1+ +ν |D y ei+1/2,1 |2 , y + α2 2 2. (5.44). i=1. 123.

(39) J Sci Comput. where ε3 is a positive parameter which will be determined later. As for the term J5 , by applying Cauchy–Schwarz inequality and Lemma 5.6, we have M−1 . ν. i=1. ≤ ≤ ≤. y v (R u )i+1/2,1/2 Dx− ei+1,1/2 x y + α2 . ν α2. M−1 . v Dx− ei+1,1/2. i=1. 1/2 2ν L y. α2. v D − ye . M−1 . 2. x( y)2. 1/2 M−1 . 1/2 |(Ru )i+1/2,1/2 |2 x. i=1. 1/2. |(Ru )i+1/2,1/2 |2 x. i=1. M−1 128L y ν ν v 2 e + |(Ru )i+1/2,1/2 |2 x, D − 16 y α22 i=1. (5.45). where the following inequality which is obtained by applying the same idea in (4.53) is used M−1 . v Dx− ei+1,1/2. 2. 2 x( y)2 ≤ 4L y D − y .. i=1. Combining (5.38), (5.39), (5.43), (5.44) and (5.45), we obtain 2 M−1 u u 2 ei+1/2,1 + ei+1/2,0 y 2ν I7 ≤ − 1+ 1− x α2 ( y + α2 ) ε3 2 i=1. +ν. α2 y + α2. . M−1 ε3 − u x y 1+ |D y ei+1/2,1 |2 2 2 i=1. L y y 17ν v 2 Dx− eu 2 + +ν D − ye 16 4α22 1/2 M−1 1/2 M−1 2ν K 2 128L y ν u 2 + |(R )i+1/2,1/2 | x + |(Ru )i+1/2,1/2 |2 x. 2 α2 α 2 i=1 i=1 (5.46) Now, by taking ε3 = I7 ≤ ν. M−1 i=1. +. 2 y α2 ,. we can infer from (5.46) that. u 2 |D − y ei+1/2,1 | 1/2 M−1 . 2ν K 2 α2. L y y 17ν x y v 2 +ν D − Dx− eu 2 + ye 2 16 4α22 1/2. |(Ru )i+1/2,1/2 |2 x. i=1. +. M−1 128L y ν |(Ru )i+1/2,1/2 |2 x. α22 i=1. (5.47) To estimate the terms I8 , I9 , I10 and I11 , we have I8 ≤. L y y v 2 D − ye , 2κ. I9 ≤. ε2 νDx− eu 2. 123. (5.48) . ε2 ν L 2x ˜ u 2 L y L y + v 2 + D − R + ye + 2 ε2 ν 2ε2 ν. y 2. R˜ v 2.

(40) J Sci Comput. +. M−1 ν − u x y |D y ei+1/2,N +1 |2 2 2 i=1. +. +. ν 2. N −1 . v 2 − v 2 |Dx− e1, j+1/2 | + |Dx e M+1, j+1/2 |. x y. j=1. 2. M−1. 2 x y ν ( y)2 u yy xi+1/2 , y N 32 2 i=1. −1 2 N . 2 . 2 x y ν ( x) vx x x1 , y j+1/2 , (5.49) + vx x x M+1 , y j+1/2 32 2 j=1 1/2 R φ + R d , ≤ K3 (5.50) L y L y + y 2 ε2 ν L 2x v 2 + u 2 v 2 ≤ ε2 νDx− eu 2 + e + R + D − D D − y x y R , 2 ε2 ν 2ε2 ν (5.51) +. I10 I11. where R˜ u and R˜ v are defined as. u Ri+1/2, u j, R˜ i+1/2, j = u Ri+1/2, j − ν4 u yy (xi+1/2 , y N ),. if j = N , if j = N ,. and. R˜ i,v j+1/2,. =. Ri,v j+1/2 , Ri,v j+1/2 − ν4 vx x (xi , y j+1/2 ),. if i = 1, M, if i = 1, M.. Combining the estimates of Ii , i = 6, . . . , 11 and taking ε2 = 41 , we can infer from (5.36) that . −1 M−1 2 N 3 L y y u − v D− − νDx− eu 2 + ν x y y ei+1/2, j+1 + Dx ei+1, j+1/2 2 4α2 i=1 j=1. +. ν 2. M−1 . u 2 |D − y ei+1/2,N +1 |. i=1. x y 2. ν v 2 − v 2 x y |Dx− e1, j+1/2 | + |Dx e M+1, j+1/2 | 2 2 j=1 L y y 7ν v 2 − 2 − 2 D − − + y e + κDx φ + κD y φ 16 2κ 1/2 M 1/2 M−1 2ν K 2 Ly p ≤ |(R )i,1/2 |2 x + |(Ru )i+1/2,1/2 |2 x ν α2 i=1 i=1 y M−1 2 L + 2L y y 128L y ν 2 4L x ˜ u 2 + |(Ru )i+1/2,1/2 |2 x + R + R˜ v 2 ν ν α22 +. N −1 . i=1. 123.

(41) J Sci Comput. +. M−1 ν( y)2 x y (u yy (xi+1/2 , y N ))2 32 2 i=1. N −1 ν( x)2 . x y. ((vx x (x1 , y j+1/2 ))2 + vx x (x M+1 , y j+1/2 ))2 2 j=1 2L y L y + y 2 4L 2x 1/2 v 2 +K 3 (R φ + R d ) + (5.52) Dx+ R u 2 + D − y RS . ν ν By taking y to satisfy the condition (4.35) and applying Lemma 5.5, the proof of Theorem 5.1 is complete.

(42). +. 32. To show the error estimates of the pressure p and φ, we need the following lemmas: Lemma 5.7 Let the mesh widths satisfy (4.35) and y = x, we have 2 u u M−1 ei+1/2,1 + ei+1/2,0 x = O(( x)2 + ( y)2 ). 2 i=1. The proof of Lemma 5.7 follows Lemma 4.7 and Theorem 5.1. Like Lemmas 4.7 and 4.6, we also need the following lemma related to the finite difference method for the divergence problem in the whole domain : φ. p. Lemma 5.8 For any given ei, j , 1 ≤ i ≤ M, 1 ≤ j ≤ N and ei, j , 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0 satisfying (5.21), there exist two vectors (eu˜ , ev˜ ) satisfies the following properties: u˜ u˜ e1/2, j = e M+1/2, j = 0, 0 ≤ j ≤ N + 1, u˜ u˜ ei+1/2,N +1 = −ei+1/2,N , 0 ≤ i ≤ M, v˜ v˜ v˜ v˜ e0, j+1/2 = −e1, j+1/2 , e M+1, j+1/2 = −e M, j+1/2 , 0 ≤ j ≤ N , u˜ − v˜ Dx− ei+1/2, j + D y ei, j+1/2 = ei, j − p. Rp , 1 ≤ i ≤ M, 1 ≤ j ≤ N , 2L x L y. φ˜. and there exists a vector ei, j , 0 ≤ i ≤ M + 1, −N ≤ j ≤ 1 satisfying φ˜. φ˜. Dx− e1, j = Dx− e M+1, j = 0, −N + 1 ≤ j ≤ 0, φ˜. φ˜. − v˜ D− y ei,−N +1 = 0, −κ D y ei,1 = ei,1/2 , 1 ≤ i ≤ M, φ˜. −κ. φ˜. Dx− ei+1, j − Dx− ei, j. φ˜. −κ. φ˜. − D− y ei, j+1 − D y ei, j. x 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0.. y. φ. = ei, j −. Rp , 2L x L y. Moreover, we have u˜ 2 − v˜ 2 − v˜ 2 Dx− eu˜ 2 + D − y e + Dx e + D y e u˜ 2. v˜ 2. +e + e . ˜ + Dx− eφ 2. φ˜ 2 + D − ye . +. M−1 . . u˜ u˜ + ei+1/2,0 ei+1/2,1. i=1. ≤ C˜ e (e p 2 + eφ 2 + (R p )2 ), where C˜ e is a constant independent of the mesh widths x and y.. 123. 2. 2 x.

(43) J Sci Comput. The proof of Lemma 5.8 can be proven similarly by following the process of showing Lemma 4.6, which needs the following lemma: φ. p. Lemma 5.9 For any given ei, j , 1 ≤ i ≤ M, 1 ≤ j ≤ N and ei, j , 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0 satisfying (5.21), there exist a positive constant Ce independent of the mesh widths U V x and y and two vectors ei+1/2, j , 0 ≤ i ≤ M, −N ≤ j ≤ N + 1 and ei, j+1/2 , 0 ≤ i ≤ M + 1, −N + 1 ≤ l ≤ N satisfying U = eU e1/2,l M+1/2,l = 0, −N ≤ l ≤ N + 1,. (5.53). U U U U (5.54) ei+1/2,N +1 = −ei+1/2,N , ek+1/2,−N = −ei+1/2,−N +1 , 0 ≤ i ≤ M, V V (5.55) ei,N +1/2 = ei,−N +1/2 = 0, 0 ≤ i ≤ N + 1, V V V V e0, (5.56) j+1/2 = −e1, j+1/2 , e M+1, j+1/2 = −e M, j+1/2 , −N + 1 ≤ j ≤ N , φ R p U − V , 1 ≤ i ≤ M, 1 ≤ j ≤ N , (5.57) Dx− ei+1/2, j + D y ei, j+1/2 = ei, j − 2L x L y Rφ φ U − V , 1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0, (5.58) Dx− ei+1/2, j + D y ei, j+1/2 = ei, j − 2L x L y U 2 − V 2 − V 2 p 2 φ 2 p 2 Dx− eU 2 + D − y e + Dx e + D y e ≤ C e (e + e + (R ) ).. (5.59) The proof of Lemma 5.9 can be proven similarly in [29]. Thus, the proofs of Lemma 5.8 and 5.9 are omitted. Now, we are ready to state the error estimates of the pressure p and φ as follows: Theorem 5.2 Let the mesh widths satisfy (4.35) and y = x, we have e p 2 + eφ 2 = O(( x)2 + ( y)2 ). The proof of Theorem 5.2 can be shown by following the similar procedure in Theorem 4.2 with the aid of Lemmas 5.5, 5.7, 5.9, 5.8 and triangle inequality.. 6 Numerical Results In this section, we carry out some numerical tests for the present MAC scheme of the Stokes/Darcy coupling problem. The computational domain f = [0, 1] × [1, 2], p = [0, 1] × [0, 1] and the interface is located at y = 1 in the domain. For simplicity, all physical ˜ α1 are all equal to 1. Throughout the paper, we choose the mesh widths constants ν, κ, k, x = y so the grid size M = N . The discrete L 2 error norms for all variables u, v, p, φ are computed based on the formulas (4.1)–(4.4). Example 1 The first analytic solution is taken from the initial condition of the test example constructed in [11] written by 1 u = − e y sin(π x), π. v = e y − e cos (π x) , p = 2e y cos(π x),. φ = e y − ye cos(π x).. 123.

(44) J Sci Comput Table 1 Grid refinement analysis of L 2 errors for the solutions u and v in Example 1. Table 2 Grid refinement analysis of L 2 errors for the solutions p and φ in Example 1. M × 2N. u exact − u. Rate. vexact − v. Rate. 16 × 32. 2.698 × 10−3. –. 3.549 × 10−3. –. 32 × 64. 6.884 × 10−4. 1.97. 9.068 × 10−4. 1.97. 64 × 128. 1.734 × 10−4. 1.99. 2.287 × 10−4. 1.99. 128 × 256. 4.349 × 10−5. 2.00. 5.734 × 10−5. 2.00. M × 2N. pexact − p. Rate. φexact − φ. Rate. 16 × 32. 7.310 × 10−2. –. 2.083 × 10−3. –. 32 × 64. 2.373 × 10−2. 1.62. 5.248 × 10−4. 1.99. 64 × 128. 7.233 × 10−3. 1.71. 1.315 × 10−4. 2.00. 128 × 256. 2.112 × 10−3. 1.78. 3.291 × 10−5. 2.00. One can easily check that this solution satisfies all three interface conditions (2.10)–(2.12). In addition, the solution satisfies zero normal velocity (v = 0) across the interface y = 1. Tables 1 and 2 show the grid refinement analysis of L 2 errors for the velocity u, v and the pressure p, φ, respectively. One can see both velocity components u, v in Stokes region and the pressure φ in Darcy’s region behave like second-order convergent. However, the pressure p in Stokes region is better than first-order but not exactly second-order convergent. The numerical results show better convergence rates than the ones obtained from the present theoretical analysis. Example 2 The second test example is given by u = (y − 1)2 + x(y − 1) + 3x − 1, v = x(x − 1) − 0.5(y − 1)2 − 3y + 1, p = 2x + y − 1, (y − 1)3 + 2x + 2y + 4. 3 Again, the solution satisfies all three interface conditions (2.10)–(2.12) but unlike Example 1, here, the normal velocity v across the interface y = 1 is nonzero. Tables 3 and 4 show the grid refinement analysis of L 2 errors for the velocity u, v and the pressure p, φ, respectively. As in Example 1, both velocity components u, v in Stokes region and the pressure φ in Darcy’s region behave like second-order convergent. However, the pressure p in Stokes region behaves exactly first-order convergent. Again, the numerical results show better convergence rates than the ones obtained from the present theoretical analysis. φ = x(1 − x)(y − 1) +. 7 Concluding Remarks The stability and error estimates for both velocity and pressure have been established for the MAC scheme of stationary Stokes/Darcy coupling problem based on finite difference methods. The stability of the velocity in both Stokes and Darcy regions is derived by performing. 123.

(45) J Sci Comput Table 3 Grid refinement analysis of L 2 errors for the solutions u and v in Example 2. Table 4 Grid refinement analysis of L 2 errors for the solutions p and φ in Example 2. M × 2N. u exact − u. Rate. vexact − v. Rate. 16 × 32. 2.308 × 10−4. –. 3.748 × 10−4. –. 32 × 64. 6.069 × 10−5. 1.93. 9.343 × 10−5. 2.00. 64 × 128. 1.618 × 10−5. 1.91. 2.242 × 10−5. 2.06. 128 × 256. 4.759 × 10−6. 1.77. 5.606 × 10−6. 2.00. M × 2N. pexact − p. Rate. φexact − φ. Rate. 16 × 32. 1.686 × 10−1. –. 1.815 × 10−4. –. 32 × 64. 8.277 × 10−2. 1.03. 4.637 × 10−5. 1.97. 64 × 128. 4.098 × 10−2. 1.01. 1.171 × 10−5. 1.99. 128 × 256. 2.039 × 10−2. 1.01. 2.951 × 10−6. 1.99. careful estimates and the stability of the pressure in both regions is established by considering an analogue of discrete divergence problem . We remark that the mesh width y needs to be below a threshold to make sure the stability. However, there is no limitation on the mesh size x. This is due to the control of the related estimates on the interface. Following the similar analysis on stability, the error estimates for the velocity and the pressure in both regions are performed. The theoretical results show first-order convergence of the scheme in discrete L 2 norm for both velocity and the pressure. The main problematic term comes form the estimate of the first order discrete approximation (3.13) of the interface condition (2.11). Moreover, in fluid region, the first-order convergence for the x-derivative of velocity component u and the y-derivative of velocity component v is also obtained in discrete L 2 norms. However, numerical tests presented in Sect. 5 show one order better for the velocity in Stokes region and the pressure in Darcy region. This is a gap between theoretical results and numerical evidences which will be studied elsewhere. Acknowledgements The work of M.-C. Lai was supported in part by Ministry of Science of Technology of Taiwan under research grant MOST-104-2115-M-009-014-MY3 while M.-C. Shiue was supported in part by the Grant MOST-104-2115-M-009-012-MY2. K.C. Ong was supported in part by NCTU Taiwan Elite Internship Program at National Chiao Tung University and NCTS.. References 1. Arbogast, T., Brunson, D.S.: A computational method for approximating a Darcy–Stokes system governing a vuggy porous medium. Comput. Geosci. 11, 207–218 (2007) 2. Tlupova, S., Cortez, R.: Boundary integral solutions of coupled Stokes and Darcy flows. J. Comput. Phys. 228, 158–179 (2009) 3. Chen, Q.: Stable and convergent approximation of two-dimensional vector fields on unstructured meshes. J. Comput. Appl. Math. 307, 284–306 (2016) 4. Chou, S.H.: Analysis and convergence of a covolume method for the generalized Stokes problem. Math. Comp. 66, 85–104 (1997) 5. Cai, M., Mu, M., Xu, J.: Preconditioning techniques for a mixed Stokes/Darcy model in porous media applications. J. Comput. Appl. Math. 233, 346–355 (2009) 6. Cao, Y., Gunzburger, M., Hu, X., Hua, F., Wang, X., Zhao, W.: Finite element approximations for Stokes– Darcy flow with Beavers–Joseph interface conditions. SIAM J. Numer. Anal. 47, 4239–4256 (2010). 123.