Efficient and Long-Time Accurate Second-Order Methods for Stokes-Darcy Systems

# Efficient and Long-Time Accurate Second-Order Methods for Stokes-Darcy Systems

Wenbin Chen School of Mathematical Sciences, Fudan University (wbchen@fudan.edu.cn).    Max Gunzburger Department of Scientific Computing, Florida State University (gunzburg@fsu.edu).    Dong Sun Department of Mathematics, Florida State University (dsun@math.fsu.edu).    Xiaoming Wang Department of Mathematics, Florida State University (wxm@math.fsu.edu).
###### Abstract

We propose and study two second-order in time implicit-explicit (IMEX) methods for the coupled Stokes-Darcy system that governs flows in karst aquifers. The first is a combination of a second-order backward differentiation formula and the second-order Gear’s extrapolation approach. The second is a combination of the second-order Adams-Moulton and second-order Adams-Bashforth 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 fully-discretized schemes using finite element spatial discretizations are derived. Numerical examples are provided that illustrate the accuracy, efficiency, and long-time stability of the two schemes.

Key words. Stokes-Darcy systems, backward differentiation formulas, Gear’s extrapolation, Adams-Moulton and Adams-Bashforth methods, unconditional stability, long-time 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 high-rain seasons, the water pressure in the conduits is larger than that in the ambient matrix so that conduit-borne 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 well-known 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, lab-size 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 Stokes-Darcy model that includes the classical Beavers-Joseph [6] matrix-conduit 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 Navier-Stokes/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 Stokes-Darcy 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, long-time 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 long-time 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 Stokes-Darcy system. We discretize the system in time via either a combination of second-order BDF and and Gear extrapolation methods or a combination of second-order Adams-Moulton and Adams-Bashforth methods. These algorithms are special cases of the implicit-explicit (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 long-time 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 long-time 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 time-dependent non-iterative version of the steady-state domain decomposition work in [13, 16] and as a generalization of the first-order schemes in [10, 37] that achieve the desirable second-order accuracy withtout increasing the complexity. The backward differentiation-based algorithm can be viewed as an infinite-dimensional version of the scheme presented in [33], but with the additional important result on time-uniform error estimates. The Adams-Moulton/Bashford based algorithm is new so far as the Stokes-Darcy problem is concerned. To the best of our knowledge, our uniform in time error estimates are the first of their kind for Stokes-Darcy and related systems.

The rest of the paper is organized as follows. In §LABEL:sec:2, we introduce the coupled Stokes-Darcy system and the associated weak formulation as well as the two second-order in time schemes. The unconditional and long-time 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 first-order 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 long-time stability of our our algorithms are given in §LABEL:sec:6. We close by providing some concluding remarks in §LABEL:sec:7.

## 2 The Stokes-Darcy system and two types of IMEX methods

### 2.1 The Stokes-Darcy 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 Stokes-Darcy system governing fluid flow in the karst system is given by [7, 16]

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩S∂ϕ∂t−∇⋅(K∇ϕ)=fin Ωp∂uf∂t−∇⋅T(uf,p)=fand∇⋅uf=0in Ωf, \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 Beavers-Joseph-Saffman-Jones interface conditions [6, 26, 39]

 ⎧⎪⎨⎪⎩uf⋅nf=up⋅nf=−(K∇ϕ)⋅nf−τj⋅(T(uf,pf)⋅nf)=αBJSJτj⋅uf,j=1,…,d−1−nf⋅(T(uf,pf)⋅nf)=gϕ, \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 Beavers-Joseph-Saffman-Jones 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

 Hf ={v∈(H1(Ωf))d∣v=0 on ∂Ωf∖Γ} Hp ={ψ∈H1(Ωp)∣ψ=0 % on ∂Ωp∖Γ} Q =L2(Ωf),W=Hf×Hp.

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 Stokes-Darcy 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 second-derivative 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

 ⟨⟨→ut,→v⟩⟩+a(→u,→v)+b(v,p)+aΓ(→u,→v) =⟨⟨⟨→f,→v⟩⟩⟩ \hb@xt@.01(2.3) b(u,q) =0,

where , , and and where . In (LABEL:weakform), we have that

 ⟨⟨→ut,→v⟩⟩=⟨ut,v⟩Ωf+gS⟨ϕt,ψ⟩Ωp,b(v,q)=−(q,∇⋅v)Ωf \hb@xt@.01(2.4) a(→u,→v)=af(u,v)+ap(ϕ,ψ)+aBJSJ(u,v) aΓ(→u,→v)=g(ϕ,v⋅nf)Γ−g(u⋅nf,ψ)Γ,⟨⟨⟨→f,→v⟩⟩⟩=⟨f,v⟩Ωf+g⟨f,ψ⟩Ωp,

where

 af(u,v)=ν( ∇u,∇v)Ωf,ap(ϕ,ψ)=g(K∇ϕ,∇ψ)Ωp aBJSJ(u,v)=αBJ(u⋅→τ,v⋅→τ)Γ.

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]

 a(→u,→u)≥(ν∥∇u∥2+gKmin∥∇ϕ∥2+αBJ∥u⋅→τ∥2Γ)≥Ca∥∇→u∥2, \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

 cs∥→v∥S≤∥→v∥≤CS∥→v∥S. \hb@xt@.01(2.6)

### 2.3 The second-order backward-differentiation scheme (BDF2)

The first scheme we introduce discretizes in time via a second-order BDF whereas the interface term is treated via a second-order explicit Gear’s extrapolation formula. We propose the following algorithm: for any and ,

 \hb@xt@.01(2.7) =⟨⟨⟨→fn+1,→v⟩⟩⟩−˜aΓ(2→un−→un−1,→v) b(un+1,q)=0,

where the artificial stabilizing term is defined as

 ast(→u,→v)=γf(u⋅nf,v⋅nf)Γ+γp(ϕ,ψ)Γ \hb@xt@.01(2.8)

with parameters and is defined as

 ˜aΓ(→u,→v)=aΓ(→u,→v)−ast(→u,→v).

### 2.4 The second-order Adams-Moulton-Bashforth method (AMB2)

For the second scheme, we combine the second-order implicit Adams-Moulton treatment of the symmetric terms and the second-order explicit Adams-Bashforth treatment of the interface term to propose the following second-order scheme: for any and ,

 ⟨⟨→un+1−→unΔt,→v⟩⟩+a(Dα→un+1,→v)+b(v,Dαpn+1)+ast(Dα→un+1,→v) \hb@xt@.01(2.9) =⟨⟨⟨→fn+12,→v⟩⟩⟩−˜aΓ(32→un−12→un−1,→v) b(Dαun+1,q)=0,

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:

1. given

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

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

Note that steps 2 and 3 can be solved independently. Moreover, legacy codes can be used in each of those steps.

## 3 Unconditional and long-time stability

The goal of this section is to demonstrate the unconditional and long-time stability, with respect to the norm, of the two second-order 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 second-order BDF is given by

 G=(12−1−152)

with the associated -norm given by The following identity is well-known (see, e.g., [23]): for any , ,

 (32v2−2v1+12v0,v2)=12(∥w1∥2G−∥w0∥2G)+∥v2−2v1+v0∥24, \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

 Cl∥w∥2G≤∥w∥2≤Cu∥w∥2GandCl∥w∥2G≤|w|2G≤Cu∥w∥2G.

We next recall the following basic inequalities:

• trace inequality: if , then

 ∥→v∥Γ≤Ctr√∥→v∥∥∇→v∥,∥→v∥Γ≤Ctr∥∇→v∥ \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

 |ast(→u,→v)|+|aΓ(→u,→v)|≤Cct∥→u∥Γ∥→v∥Γ∀→u,→v∈W.

Proof. By the definition (LABEL:ast) of , we have

 |ast(→u,→v)| ≤γf|(u⋅nf,v⋅nf)Γ|+γp|(ϕ,ψ)Γ| \hb@xt@.01(3.3) ≤γf∥u⋅nf∥Γ∥v⋅→nf∥Γ+γp∥ϕ∥Γ∥ψ∥Γ ≤γmax(∥u⋅nf∥Γ∥v⋅nf∥Γ+∥ϕ∥Γ∥ψ∥Γ),

where the triangle and Cauchy-Schwarz inequalities are used and . Similarly, by the definition (LABEL:bilinearforms) of , we have

 |aΓ(→u,→v)|≤g(∥ϕ∥Γ∥v⋅nf∥Γ+∥u⋅nf∥Γ∥ψ∥Γ). \hb@xt@.01(3.4)

Note that so that

 ∥→u∥2Γ=∥u∥2Γ+∥ϕ∥2Γ=∥u⋅nf∥2Γ+∥u⋅τ∥2Γ+∥ϕ∥2Γ.

Then, combining (LABEL:astineq) and (LABEL:agamma), we obtain

 |ast(→u ,→v)|+|aΓ(→u,→v)| ≤(γmax(∥u⋅→nf∥Γ∥v⋅→nf∥Γ+∥ϕ∥Γ∥ψ∥Γ)+g(∥ϕ∥Γ∥v⋅→nf∥Γ+∥u⋅→nf∥Γ∥ψ∥Γ)) ≤max{γmax,g}(∥u⋅→nf∥Γ+∥ϕ∥Γ)(∥v⋅→nf∥Γ+∥ψ∥Γ) ≤√2max{γmax,g}∥→u∥Γ∥→v∥Γ.

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

 1Δt⟨⟨D→un+1,→un+1⟩⟩ +a(→un+1,→un+1)+ast(δ→un+1,→un+1) =⟨⟨⟨→fn+1,→un+1⟩⟩⟩−aΓ(2→un−→un−1,→un+1).

From (LABEL:Gnormid) and the skew-symmetry of , we obtain

 12|→wn|2G− 12|→wn−1|2G+14∥δ→un+1∥2S+Δta(→un+1,→un+1)+Δtast(→un+1,→un+1) \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

 ˜aΓ(−2→un +→un−1,→un+1)≤Cct∥−2→un+→un−1∥Γ∥→un+1∥Γ \hb@xt@.01(3.6) ≤CctC2tr∥−2→un+→un−1∥12∥−2∇→un+∇→un−1∥12∥∇→un+1∥ ≤CctC2tr∥−2→un+→un−1∥12∥∇→un+1∥(√2∥∇→un∥12+∥∇→un−1∥12) ≤C12|→wn−1|2+Ca6∥∇→un+1∥2+Ca3∥∇→un∥2 +C22|→wn−1|2+Ca6∥∇→un+1∥2+Ca6∥∇→un−1∥2.

For the forcing term, we have

 ⟨⟨⟨→fn+1,→un+1⟩⟩⟩≤C32∥→fn+1∥2+Ca6C2P∥∥→un+1∥∥2. \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

 |→wn|2G +CaΔt∥∇→un+1∥2≤C3∥→fn+1∥2Δt +(1+(C1+C2)Δt)|→wn−1|2G+2CaΔt3∥∇→un∥2+CaΔt3∥∇→un−1∥2.

Next, by adding to both sides of this inequality, we deduce

 En+CaΔt2(C1+C2)(1+(C1+C2)Δt)∥∇→un+1∥2+CaΔt2(C1+C2)3(1+(C1+C2)Δt)∥∇→un∥2 ≤C3∥→fn+1∥2Δt+(1+(C1+C2)Δt)En−1 ≤e(C1+C2)TE0+C3C1+C2e(C1+C2)Tmaxn∥→fn+1∥2,

where

 En=|→wn|2G+CaΔt(1+(C1+C2)Δt)∥∇→un+1∥2+CaΔt3(1+(C1+C2)Δt)∥∇→un∥2.

Thus, the unconditional stability of the BDF2 scheme is proved.

#### 3.1.2 Unconditional stability of the AMB2 scheme

We introduce the parameters

 α1=∣∣32−2α∣∣,α2=∣∣α−12∣∣, \hb@xt@.01(3.8) β3=α1+α2,β1=2α−β3,β2=12(β1+β3).
###### 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

 ⟨⟨→un+1−→unΔt,→un+1⟩⟩+a(Dα→un+1,→un+1)+ast(Dα→un+1,→un+1) \hb@xt@.01(3.9)

Combining the two terms and using the basic equality , we have

 1Δt(∥→un+1∥2S−∥→un∥2S+∥→un+1−→un∥2S)+2a(Dα→un+1,→un+1) \hb@xt@.01(3.10) =2(→fn+12,→un+1)−2aΓ(32→un−12→un−1,→un+1)−2αast(δ→un+1,→un+1).

Note that provided . Therefore, and hence when . By the Cauchy-Schwarz inequality, we then have

 2a(Dα→un+1,→un+1)≥2αa(→un+1,→un+1)−α1(a(→un+1,→un+1)+a(→un,→un)) \hb@xt@.01(3.11) −α2(a(→un+1,→un+1)+a(→un−1,→un−1)) =β1a(→un+1,→un+1)−α1a(→un,→un)−α2a(→un−1,→un−1).

Similarly as for (LABEL:BDF:interface1), for the interface coupling term, there exists a constant such that

 −2aΓ(32→un−12→un−1,→un+1)−2αast(−2→un+→un−1,→un+1) \hb@xt@.01(3.12) ≤Ca(β1−β2)4∥∇→un∥2+2C5∥→un∥2S+Ca(β1−β2)4∥∇→un+1∥2 +Ca(β1−β2)8∥∇→un−1∥2+C5∥→un−1∥2S+Ca(β1−β2)8∥∇→un+1∥2.

For the forcing term, there exists a constant such that

 2Δt⟨⟨⟨→fn+12,→un+1⟩⟩⟩≤C6Δt∥→fn+12∥2+Ca(β1−β2)8Δt∥∇→un+1∥2. \hb@xt@.01(3.13)

Substituting (LABEL:goodtermAMB)–(LABEL:forcingAMB) into (LABEL:weakformAMB) yields

 ∥→un+1∥2S+Ca(β1−β2)2Δt∥∥∇→un+1∥∥2+Δtβ2∥→un+1∥2a \hb@xt@.01(3.14) +∥→un+1−→un∥2S+2αΔtast(→un+1,→un+1) ≤C6∥→fn+12∥2Δt+(1+2C5Δt)∥→un∥2S+C5Δt∥→un−1∥2S+Δtα1∥→un∥2a +Δtα2∥→un−1∥2a+Ca(β1−β2)4Δt∥∇→un∥2+Ca(β1−β2)8Δt∥∇→un−1∥2.

Define the energy

 En=∥→un∥2S+C5Δt1+3C5Δt∥→un−1∥2S+β3Δt1+3C5Δt∥→un∥2a+Δtα21+3C5Δt∥→un−1∥