Efficient and LongTime Accurate SecondOrder Methods for StokesDarcy Systems
Abstract
We propose and study two secondorder in time implicitexplicit (IMEX) methods for the coupled StokesDarcy system that governs flows in karst aquifers. The first is a combination of a secondorder backward differentiation formula and the secondorder Gear’s extrapolation approach. The second is a combination of the secondorder AdamsMoulton and secondorder AdamsBashforth methods. Both algorithms only require the solution of two decoupled problems at each time step, one Stokes and the other Darcy. Hence, these schemes are very efficient and can be easily implemented using legacy codes. We establish the unconditional and uniform in time stability for both schemes. The uniform in time stability leads to uniform in time control of the error which is highly desirable for modeling physical processes, e.g., contaminant sequestration and release, that occur over very long time scales. Error estimates for fullydiscretized schemes using finite element spatial discretizations are derived. Numerical examples are provided that illustrate the accuracy, efficiency, and longtime stability of the two schemes.
Key words. StokesDarcy systems, backward differentiation formulas, Gear’s extrapolation, AdamsMoulton and AdamsBashforth methods, unconditional stability, longtime stability, uniform in time error estimates, finite element methods, karst aquifers
AMS subject classifications. 35M13, 35Q35, 65N30, 65N55, 76D07, 76S05
1 Introduction
Karst is a common type of landscape formed by the dissolution of layers of soluble bedrock, usually including carbonate rock, limestone, and dolomite. Karst regions often contain karst aquifers which are important sources of potable water. For example, about 90% of the fresh water used in Florida comes from karst aquifers [27]. Clearly, the study of karst aquifers is of great importance, especially because they are seriously threatened by contamination [29].
A karst aquifer, in addition to a porous limestone or dolomite matrix, typically has large cavernous conduits that are known to have great impact on groundwater flow and contaminant transport within the aquifer. During highrain seasons, the water pressure in the conduits is larger than that in the ambient matrix so that conduitborne contaminants can be driven into the matrix. During dry seasons, the pressure differential reverses and contaminants long sequestered in the matrix can be released into the free flow in the conduits and exit through, e.g., springs and wells, into surface water systems. Therefore, the understanding of the interaction between the free flow in the conduits and the Darcy flow in the matrix is crucial to the study of groundwater flows and contaminant transport in karst region.
The mathematical study of flows in karst aquifers is a wellknown challenge due to the coupling of the flow in the conduits and the flow in the surrounding matrix, the complex geometry of the network of conduits, the vastly disparate spatial and temporal scales, the strong heterogeneity of the physical parameters, and the huge associated uncertainties in the data. Even for a small, labsize conceptual model with only one conduit (pipe) imbedded in a homogenous porous media (matrix), significant mathematically rigorous progress has only recently been achieved. For the coupled StokesDarcy model that includes the classical BeaversJoseph [6] matrixconduit interface boundary condition, see [7, 8, 9]. For various simplified interface conditions, see, e.g., [7, 16, 30]. Nonlinear interface conditions have also been proposed for NavierStokes/Darcy modeling; see, e.g., [11, 18].
Due to the practical importance of the problem of flow and contaminant transport in karst aquifers, there has been a lot of attention recently paid to the development of accurate and efficient numerical methods for the coupled StokesDarcy system; see, e.g., [10, 16, 30, 36] among many others. The efficiency of the algorithms is a particularly important issue due to the large scale of field applications. Because of the disparity of governing equations and physics in the conduit and matrix, domain decomposition methods (also called partitioned methods by some authors) that only requires separate Stokes and Darcy solves seems natural; see, e.g., [12, 13, 16, 17, 28, 31, 32, 33, 37]. On the other hand, longtime accuracy of the schemes is also highly desirable because the physical phenomena of retention and release of contaminants takes place over a very long time scale. Therefore, there is a need to ensure the longtime accuracy of the discretization algorithms in addition to the standard notion of accuracy on an order one time scale.
The purpose of this work is to propose and investigate two types of numerical methods for the coupled StokesDarcy system. We discretize the system in time via either a combination of secondorder BDF and and Gear extrapolation methods or a combination of secondorder AdamsMoulton and AdamsBashforth methods. These algorithms are special cases of the implicitexplicit (IMEX) class of schemes. The coupling terms in the interface conditions are treated explicitly in our algorithm so that only two decoupled problems (one Stokes and one Darcy) are solved at each time step. Therefore, these schemes can be implemented very efficiently and, in particular, legacy codes for each of the two components can be utilized. Moreover, we show that our schemes are unconditionally stable and longtime stable in the sense that the solutions remain uniformly bounded in time. The uniform in time bound of the solution further leads to uniform in time error estimates. This is a highly desirable feature because one would want to have reliable numerical results over the longtime scale of contaminant sequestration and release. Uniform in time error estimates for fully discrete schemes using finite element spatial discretizations are also presented. Our numerical experiments illustrate our analytical results.
Our work can be viewed as a timedependent noniterative version of the steadystate domain decomposition work in [13, 16] and as a generalization of the firstorder schemes in [10, 37] that achieve the desirable secondorder accuracy withtout increasing the complexity. The backward differentiationbased algorithm can be viewed as an infinitedimensional version of the scheme presented in [33], but with the additional important result on timeuniform error estimates. The AdamsMoulton/Bashford based algorithm is new so far as the StokesDarcy problem is concerned. To the best of our knowledge, our uniform in time error estimates are the first of their kind for StokesDarcy and related systems.
The rest of the paper is organized as follows. In §LABEL:sec:2, we introduce the coupled StokesDarcy system and the associated weak formulation as well as the two secondorder in time schemes. The unconditional and longtime stability with respect to the norm are presented in §LABEL:sec:3. Section LABEL:sec:4 is devoted to the stability with respect to the norm. The estimates are important for the finite element analysis; this is another new feature of our work, even for firstorder schemes. In §LABEL:sec:5, we focus on the error analysis of the fully discretized scheme using finite element spatial discretizations. Numerical results that illustrate the accuracy, efficiency, and longtime stability of our our algorithms are given in §LABEL:sec:6. We close by providing some concluding remarks in §LABEL:sec:7.
2 The StokesDarcy system and two types of IMEX methods
2.1 The StokesDarcy system
For simplicity, we consider a conceptual domain for a karst aquifer that consists of a porous media (matrix), denoted by , and a conduit, denoted by , where denotes the spatial dimension. The interface between the matrix and the conduit is denoted . The remaining parts of the boundaries of and are denoted by and , respectively. See Fig. LABEL:fig:1.
The coupled StokesDarcy system governing fluid flow in the karst system is given by [7, 16]
\hb@xt@.01(2.1) 
where the unknowns are the fluid velocity and the kinematic pressure in the conduit and the hydraulic head in the matrix; the velocity in the matrix is recovered from . In (LABEL:SDsys), and denote external body forces acting on the domains and respectively, and denotes the stress tensor. The parameters appearing in (LABEL:SDsys) are the water storage coefficient , the hydraulic conductivity tensor , and the kinematic viscosity of the fluid .
For simplicity, we assume homogeneous Dirichlet boundary conditions for the hydraulic head and fluid velocity on the outer boundaries and , respectively. On the interface , we impose the BeaversJosephSaffmanJones interface conditions [6, 26, 39]
\hb@xt@.01(2.2) 
where denotes the outer unit normal vector for and , denotes a linearly independent set of vectors tangent to the interface . The additional parameters appearing in (LABEL:BJSJ) are the gravitational constant and the BeaversJosephSaffmanJones coefficient .
2.2 Weak formulation
We denote by and the standard inner product and norm, respectively, where may be , , or . We often suppress the subscript if there is no possibility of confusion. We define the spaces
Dual spaces are denoted by and duality parings between spaces and their duals induced by the inner product on the appropriate domain are denoted by .
A weak formulation of the StokesDarcy system (LABEL:SDsys) is derived by multiplying the three equations in that system by test functions , , and , respectively, then integrating over the corresponding domains, then integrating by parts the terms involving secondderivative operators, and then substituting the interface conditions (LABEL:BJSJ) in the appropriate terms. The resulting weak formulation is given as follows [7, 15]: given and , seek , , and , with and , satisfying
\hb@xt@.01(2.3)  
where , , and and where . In (LABEL:weakform), we have that
\hb@xt@.01(2.4)  
where
In (LABEL:weakform), , , and are the primary variables; as mentioned before, once the hydraulic head is known, one can recover , the velocity in the porous media, via the Darcy relation .
It can be shown that the bilinear form is coercive; indeed, we have that [7]
\hb@xt@.01(2.5) 
where and where denotes the smallest eigenvalue of . We define the norms and . We have that is equivalent to the norm, i.e., we have that
\hb@xt@.01(2.6) 
2.3 The secondorder backwarddifferentiation scheme (BDF2)
The first scheme we introduce discretizes in time via a secondorder BDF whereas the interface term is treated via a secondorder explicit Gear’s extrapolation formula. We propose the following algorithm: for any and ,
\hb@xt@.01(2.7)  
where the artificial stabilizing term is defined as
\hb@xt@.01(2.8) 
with parameters and is defined as
2.4 The secondorder AdamsMoultonBashforth method (AMB2)
For the second scheme, we combine the secondorder implicit AdamsMoulton treatment of the symmetric terms and the secondorder explicit AdamsBashforth treatment of the interface term to propose the following secondorder scheme: for any and ,
\hb@xt@.01(2.9)  
where denotes the difference operator that depends on a parameter and is defined by . The stabilizing term is defined as in (LABEL:ast).
2.5 Efficiency of the schemes
The implemented schemes are highly efficient because we can decouple the Stokes and Darcy subproblems:

given

set so that all terms involving drop out and we only need to use a fast Stokes solver to obtain ;

set so that all terms involving drop out and we only need a fast Poisson solver to obtain ;

Set and return to step 1.
Note that steps 2 and 3 can be solved independently. Moreover, legacy codes can be used in each of those steps.
3 Unconditional and longtime stability
The goal of this section is to demonstrate the unconditional and longtime stability, with respect to the norm, of the two secondorder schemes proposed in §LABEL:sec:2. We first recall a few basic facts and notations that are needed below.
Recall that the matrix associated with the classical secondorder BDF is given by
with the associated norm given by The following identity is wellknown (see, e.g., [23]): for any , ,
\hb@xt@.01(3.1) 
where and . We also apply the matrix to functions belonging to : for any , define Then, for any , ,
where and .
The norms are equivalent norms on in the sense that there exists such that
We next recall the following basic inequalities:

trace inequality: if , then
\hb@xt@.01(3.2) 
Poincaré inequality: if , then

Young inequality:
Other variants of Young’s inequality will also be used.
The following estimate follows from the basic inequalities.
Lemma 3.1
Let and be defined as in (LABEL:bilinearforms) and (LABEL:ast), respectively. Then, there exists a constant such that
Proof. By the definition (LABEL:ast) of , we have
\hb@xt@.01(3.3)  
where the triangle and CauchySchwarz inequalities are used and . Similarly, by the definition (LABEL:bilinearforms) of , we have
\hb@xt@.01(3.4) 
Note that so that
Then, combining (LABEL:astineq) and (LABEL:agamma), we obtain
The lemma is proved by setting .
For the sake of brevity, we introduce the BDF difference operator and the central difference operator .
3.1 Unconditional stability of the the BDF2 and AMB2 schemes
3.1.1 Unconditional stability of the BDF2 scheme
xxx
Theorem 3.2
Let be any fixed time. Then, the BDF2 scheme (LABEL:schemeBDF) is unconditionally stable on .
Proof. Setting in the BDF2 scheme (LABEL:schemeBDF), we have
From (LABEL:Gnormid) and the skewsymmetry of , we obtain
\hb@xt@.01(3.5)  
where . Also, from the definition of the bilinear form , Lemma LABEL:Lemma1, the trace inequality, and Young’s inequality, we have
\hb@xt@.01(3.6)  
For the forcing term, we have
\hb@xt@.01(3.7) 
After we discard the nonnegative terms and , noting that , and using (LABEL:BDF:interface1) and (LABEL:forcingterm), (LABEL:ineq1) becomes
Next, by adding to both sides of this inequality, we deduce
where
Thus, the unconditional stability of the BDF2 scheme is proved.
3.1.2 Unconditional stability of the AMB2 scheme
We introduce the parameters
\hb@xt@.01(3.8)  
Theorem 3.3
Let be any fixed time and let . Then, the AMB2 scheme (LABEL:schemeAMB) is unconditionally stable in .
Proof. Setting in (LABEL:schemeAMB), we deduce
\hb@xt@.01(3.9)  
Combining the two terms and using the basic equality , we have
\hb@xt@.01(3.10)  
Note that provided . Therefore, and hence when . By the CauchySchwarz inequality, we then have
\hb@xt@.01(3.11)  
Similarly as for (LABEL:BDF:interface1), for the interface coupling term, there exists a constant such that
\hb@xt@.01(3.12)  
For the forcing term, there exists a constant such that
\hb@xt@.01(3.13) 
Substituting (LABEL:goodtermAMB)–(LABEL:forcingAMB) into (LABEL:weakformAMB) yields
\hb@xt@.01(3.14)  
Define the energy