A hybridizable discontinuous Galerkin method for the fully coupled time-dependent Stokes/Darcy-transport problem

We present a high-order hybridized discontinuous Galerkin (HDG) method for the fully coupled time-dependent Stokes-Darcy-transport problem where the fluid viscosity and source/sink terms depend on the concentration and the dispersion/diffusion tensor depends on the fluid velocity. This HDG method is such that the discrete flow equations are compatible with the discrete transport equation. Furthermore, the HDG method guarantees strong mass conservation in the $H^{\rm div}$ sense and naturally treats the interface conditions between the Stokes and Darcy regions via facet variables. We employ a linearizing decoupling strategy where the Stokes/Darcy and the transport equations are solved sequentially by time-lagging the concentration. We prove well-posedness and optimal a priori error estimates for the velocity and the concentration in the energy norm. We present numerical examples that respect compatibility of the flow and transport discretizations and demonstrate that the discrete solution is robust with respect to the problem parameters.


The Stokes/Darcy-transport system
Let Ω ⊂ R dim , dim = 2, 3, be a bounded polygonal domain. We denote its boundary by ∂Ω and the outward Figure 1. A depiction of a domain Ω ⊂ R 2 with its two sub-domains Ω s and Ω d . unit normal to ∂Ω by n. Domain Ω consists of two non-overlapping polygonal regions, a free flow region Ω s and a Darcy flow region Ω d , such that Ω = Ω s ∪ Ω d . The polygonal interface between Ω s and Ω d is denoted by Γ I and the external boundary of Ω j is denoted by Γ j := ∂Ω ∩ ∂Ω j , j = s, d. See Figure 1 for a depiction of a domain when dim = 2.
We denote the time interval of interest by J = [0, T ]. The fully coupled Stokes/Darcy-transport system for the velocity field u : Ω × J → R dim , fluid pressure p : Ω × J → R and concentration c : Ω × J → R is given by in Ω s × J, (1a) where ε(u) := ∇u + (∇u) T /2 is the strain rate tensor and χ d is the characteristic function that takes the value 1 in Ω d and 0 in Ω s . Here the fluid viscosity µ, the matrix K = κ µ , where κ is the permeability matrix of the porous medium, and the body force terms f s and f d are concentration dependent functions. The porosity φ of the medium in Ω d is a spatially varying function. In Ω s we set φ = 1. The functions g i and g p denote the source and sink terms related to injection and production wells and c I is the injected concentration. Furthermore, in the Stokes region D(u) = dI, where I is the dim × dim identity matrix and d is the diffusion coefficient. In the Darcy region D(u) = D(u), where D(u) denotes the diffusion dispersion tensor in Ω d .
We will denote the restriction of the velocity u, pressure p, and concentration c to sub-domain Ω j , j = s, d by, respectively, u j , p j , and c j . Then, on the interface Γ I , choosing the unit normal vector n to be pointing from Ω s to Ω d , we prescribe the following interface conditions that hold for t ∈ J: −2(ε(u s )n) · τ = γ u s · τ , = 1, . . . , dim − 1, (2c) where τ , = 1, . . . , dim − 1 denote the unit tangent vectors on Γ I . These conditions enforce the normal continuity of the velocity (2a), the normal continuity of the normal component of the stress (2b), continuity of the concentration (2d), and normal continuity of the concentration flux (2e). Equation (2c), where γ = α/ √ τ · κτ with α > 0 a constant, is the Beavers-Joseph-Saffman law which enforces a condition on the tangential component of the normal stress [5,46]. To close the model, we assume the following initial conditions: c(x, 0) = c 0 (x) in Ω.
The body force functions f s and f d are assumed to be Lipschitz continuous in c with Lipschitz constants L s f and L d f . Note that f s and f d depend on x and c, but they do not depend explicitly on t. We will further assume that 0 ≤ c I ≤ 1 a.e. in Ω d and that g i , g p ≥ 0, g i , g p ∈ L ∞ (J; L 2 (Ω d )) are such that A weak formulation of the problem defined by eqs. (1) to (3) was presented in [45]. The analysis for a weak solution of a more general version of this problem, in which the free fluid flow is governed by the Navier-Stokes equations, can be found in [14].

Preliminaries
We use the same notation that we used previously in [16,17]. Let T j := {K} be a shape-regular triangulation of Ω j , j = s, d, into non-overlapping elements (we only consider simplices) such that T s and T d match at the interface Γ I . We define the triangulation of the entire domain Ω as T := T s ∪ T d . The maximum diameter over all elements is h = max K∈T h K , where h K stands for the diameter of an element K. The boundary of an element K and its outward unit normal are denoted by ∂K and n, respectively. A facet of an element boundary is an interior facet if it is shared by two neighboring elements and it is a boundary facet if it is a part of ∂Ω. The set of all interior facets and all boundary facets inΩ j are denoted by F j i and F j b , j = s, d, respectively. We also collect the facets that lie on the interface Γ I in the set F I . The set of all facets that lie inΩ and inΩ j are denoted by F and F j , respectively. We point out that F j = F j i ∪ F j b ∪ F I , j = s, d. Furthermore, we define Γ 0 := ∪ F ∈F F and Γ j 0 := ∪ F ∈F j F , j = s, d. The finite element function spaces on Ω for the velocity and pressure are given by where P k (K) denotes the space of polynomials of degree at most k defined on the element K. The finite element spaces for the velocity and pressure traces are given bȳ Here P k (F ) denotes the space of polynomials of degree at most k defined on the facet F . Note that functions in V h are not defined on Γ d 0 \Γ I . The finite element function spaces for the concentration and its trace are defined as The semi-discrete and fully-discrete HDG methods for the flow and transport equations considered in this article are compatible when k c = k f − 1 [16]. For this reason we set k c = k f − 1.
To reduce the notational burden, we define Next, let us define the function spaces and set X := V × Q. As before, we use a superscript j to specify the restriction of these spaces to Ω j , j = s, d. The trace spaces of V restricted to Γ s 0 , Q j restricted to Γ j 0 , and C restricted to Γ 0 are denoted by, respectively, V ,Q j , andC. The trace operator γ V : V s →V restricts functions in V s to Γ s 0 , and similarly the trace operators γ Q j : Q j →Q j restrict functions in Q j to Γ j 0 , j = s, d. However, when it is clear from the context, we omit the subscript in the trace operator. Analogous to the discrete case, we introduce V := V ×V , Q := Q ×Q s ×Q d , and C := C ×C. We then define extended function spaces as We close this section by listing various norms on the spaces described above. We refer the reader to [1] for the definitions of the standard Sobolev spaces W m,p (D) and their corresponding norms · W m,p (D) . For ease of notation, we write · m,p,D instead of · W m,p (D) with the following simplifications. When m = 0, W 0,p (D) coincides with L p (D) and when p = 2, H m (D) = W m,p (D). For p = 2, we write · m,D to denote · W m,2 (D) and for m = 0, p = 2, we write · D instead of · 0,D .
On V s (h) we define the standard HDG-norm and its strengthened version as follows: On V (h) we then introduce the norms and note that |||v||| v and |||v||| v are equivalent on V h due to the fact that |||·||| v,s and |||·||| v ,s are equivalent on V s h (see, for example, [52, eq. (5.5)]).
On the pressure spaces Q j (h), j = s, d and Q(h), we define, respectively, Finally, for w h ∈ C(h), we define the following semi-norm:

Semi-discrete HDG scheme
The semi-discrete method we propose for the Stokes/Darcy-transport system in eqs. (1) and (2) is as follows: for all v h , q h ∈ X h and w h ∈ C h . The form B sd h in eq. (12a) collects the discretization terms for the Stokes/Darcy momentum and mass conservation equations as follows: Here a h (·, ·) is defined as where and β s > 0 is a penalty parameter. The bilinear forms b j h (·, ·) and b I,j h (·, ·) in eq. (13), j = s, d are defined as Before defining the terms related to the transport equation, we point out that eqs. (13) to (15) are the same as in [17] when the viscosity and κ are both constants. The form B tr h (u; c(t), w) in eq. (12b) discretizes the advective and diffusive parts of the transport equation: The advective part is defined as (17) where ∂K in denotes the inflow portion of the boundary on which u h · n < 0, and the diffusive part is defined as where β tr > 0 is a penalty parameter. Here we pause again to mention that eqs. (16) to (18) are the same as in [16], and a standard extension of the discretization analyzed in [52].
To complete the discretization, we project the initial conditions u 0 and c 0 eq. (3) into V h and C h , respectively.

Properties of the numerical scheme
The semi-discrete HDG scheme presented in Section 3.2 has various attractive features. Besides local momentum conservation, a property of all HDG methods, this particular HDG method also conserves mass strongly, according to the definition defined in [35].
To be specific, the discrete velocity enjoys the following properties: where · is the usual jump operator and n is the unit normal vector on F . Note that eqs. (19b) and (19c) imply that u h is H(div)-conforming on the whole domain. More details on eq. (19) can be found in [17,Section 3.3]. Additionally, the scheme is consistent, that is, the solution to eqs. (1) to (3) satisfies eq. (12), as we discuss next.
Lemma 3.1 (Consistency). Suppose that the solution (u, p, c) to the Stokes/Darcy-transport system eqs.
Proof. The proof is as in [17,Lemma 1] and [16,Lemma 6] with minor modifications. We do not repeat the proof here, but mention that it is based on integration by parts in a s h (c; u, v), j=s,d b j h (p j , v), B a h (u, c, w), and B d h (u, c, w), using γ(u) = u on Γ s 0 , γ(p j ) = p j on Γ j 0 , j = s, d, γ(c) = c on Γ 0 , the smoothness of the solution, the continuity of µ and D, and eqs. (1) and (2).

Continuity, coercivity, and an inf-sup condition
Let us recollect various known inequalities. Throughout this article we denote by C > 0 a generic constant that is independent of the mesh size and the time step. From [23, Lemma 1.46, Remark 1.47], for any K ∈ T , we have v ∂K ≤ Ch We will also use the following versions of the continuous trace inequality : where C in eq. (22) Furthermore, by [31,Theorem 4.4], for any Similarly, for any w h ∈ C h , for k c ≥ 1, The following Poincaré-type inequality follows from [ v Ω s ≤ C|||v||| v,s ∀v : The following version of Korn's first inequality is a consequence of [7, (1.19)], [31,Proposition 4.7], [42, p.110]: Continuity and coercivity a h (·, ·), follow from [17, Lemma 2, Lemma 3] keeping in mind that µ satisfies eq. (5b) and K satisfies eq. (7). They can be stated as follows: Lemma 4.1 (Continuity and coercivity of a h ). There exists a constant C > 0, independent of h, such that for all u, v ∈ V (h) and c ∈ C(h), In addition, there exists a constant C a > 0, independent of h but dependent on κ * , µ * , and µ * , and a constant β s The inf-sup condition on the discrete spaces V h and Q h was proved in [17] in the case of a continuous discrete velocity trace It is straightforward to show that the inf-sup condition also holds when the discrete velocity trace space is the larger discontinuousV h space.
There exists a constant c inf > 0, independent of h, such that for any q h ∈ Q h , Note that the proof of this theorem as well as the error analysis requires appropriate interpolation operators where we remark that F is an edge if dim = 2 and a face if dim = 3. Furthermore, for u ∈ W 1,∞ (K), [32, (2.33)], The interpolant onto the trace spaceV h is defined byΠ V : where PV is the L 2 -projection ontoV h . It is straightforward to deduce the following estimates using eq. (21) and the fact that We finish this section by noting that by eqs. (34a) to (34c) the solution u of eqs. (1) and (2) under the assumption that u ∈ [H k+1 (K)] dim , for all K ∈ T , satisfies

Fully discrete numerical scheme
Let us now describe the fully discrete HDG method and decoupling strategy used to solve the Stokes/Darcy and transport problems sequentially. For the time discretization, we partition the time interval J as: 0 = t 0 < t 1 < . . . < t N = T . For simplicity, we assume a uniform partition with t n+1 − t n = ∆t for 0 ≤ n ≤ N − 1. We denote a function h(t) at time level t n by h n := h(t n ) and for a sequence {u n } n≥1 we denote by d t u n = (u n − u n−1 )/∆t a first order difference operator.
In the first step of our sequential algorithm, given an initial velocity u 0 h in the Stokes domain and an initial concentration c 0 h , we solve the Stokes/Darcy problem and obtain a velocity in the entire region. This velocity, with properties given by eq. (19), is then substituted into the concentration problem. This approach is repeated for all time steps with the initial velocity and concentration being replaced by the last computed velocity and concentration solutions. We summarize the fully discrete problem in Algorithm 1.

Algorithm 1 Sequential algorithm
end for Remark 5.1. In Algorithm 1, Π C denotes the L 2 -projection onto C h and Π V u 0 is understood as Π V applied to the extension of u 0 to Ω by zero assuming u 0 ∈ H 1 0 (Ω s ). We note that this choice of u 0 h satisfies normal continuity across the interfaces in Γ s 0 and has zero divergence in Ω s under the additional assumption that ∇ · u 0 = 0 in Ω s . We further remark that the properties in eq. (19) hold for u n h for each time step n = 0, . . . , N .
We conclude this section by stating some preliminary results obtained by Taylor's theorem [15,Lemma 3.2]. For a function z defined on D × [0, T ], assuming enough regularity, we have the following results: Note that the inequalities in eq. (41) for = 0 were presented in [15,Lemma 3.2] and that it is straightforward to extend eq. (41b) to = 1.

Main results
In this section we present our main results. For the error estimates we will make use of the following definition of the discrete in time norm: Before proving a priori error estimates for the discrete velocity, pressure, and concentration, we first state the well-posedness of the discrete Stokes/Darcy problem eq. (39). Well-posedness of the discrete transport problem eq. (40) is proven in Section 6.3 as it depends on results obtained in Section 6.1.
s be as in Lemma 4.1 and n ≥ 1. Then given Proof. The result follows by applying the abstract theory for saddle point problems [27, Theorem 2.34] together with Theorem 4.2 and eq. (32).

Error estimates for the discrete velocity
In this section we derive estimates for the error u n h − u n , for each n ≥ 0, given error estimates for the discrete concentration in previous time steps. To do so, we define the following: For the case n = 0, Π V u 0 is understood as Π V applied to the extension of u 0 to Ω by zero assuming u 0 ∈ H 1 0 (Ω s ). Therefore, ζ 0 u = 0. Furthermore, note that the following identities hold: To be consistent with the notation used in previous sections, we set n u = ( n u ,¯ n u ), n p := ( n p ,¯ sn p ,¯ dn p ), and jn p := ( n p ,¯ jn p ), for = ξ, ζ and j = s, d. Here we recall the following results on the interpolation errors [17, Lemma 7, Lemma 8] Lemma 6.2. Let (u, p) be the velocity solution of eqs.
(1) to (3),ū = γ(u s ). Then for any n ≥ 1, Proof. The proof is based on the properties of the numerical scheme eq. (19), the properties of the BDM projection Π V in eq. (34), and the properties of the (47) is exactly the same as eq. (38), evaluated at t = t n , and rewritten by using the definitions of ξ n u andξ n u .
Theorem 6.3 (Error equation for eq. (39)). There holds Subtracting this equation from eq. (39), we obtain for all (v h , q h ) ∈ X h . Then, by eq. (43), the B sd h terms can be rewritten as where we applied eqs. (46) and (47). Combining this with eq. (49) yields the result.
The velocity error at each time step depends on the error in concentration from the previous time step. Therefore, for the velocity error estimates, we will need some auxiliary results related to the concentration error. To estimate the error of the concentration, we use the continuous interpolant Ic ∈ C h ∩ C 0 (Ω) of c [8], and we setĪc(t) = Ic| Γ 0 (t) ∈C h . Denoting the restriction of c to Γ 0 byc, we define Note that: Furthermore, we have the following interpolation estimate [8, Section 4.4] for r = 0, 1 and 1 ≤ ≤ k c : such that ∂ t c s ∈ L 2 (0, T, H 1 (Ω s )), andc = γ(c) on Γ 0 . Then we have the following estimates: Proof. We first prove eq. (54a). By the triangle inequality and the definition of d t , Multiplying this by ∆t, summing from m = 1 to n, and using eq. (41b), we obtain eq. (54a). We now prove eq. (54b). We have (24), the first term on the right side of eq. (55) is bounded as follows: Splitting the second term on the right side of eq. (55) by using Ic =Īc for any F ∈ F I gives: The first term on the right hand side of eq. (57) is bounded by eq. (21) and eq. (53) as follows: Using eq. (26) and the definition of |||·||| c , Collecting eqs. (55) to (59), we obtain Equation (54b) now follows after multiplying the above inequality by ∆t, summing from 1 to n, and using eq. (41b). We next prove eq. (54c). By the triangle inequality, The results follows by multiplying eq. (60) by ∆t, summing from m = 1 to n, and using eq. (41b) as before.
Now that the auxiliary result is established, we proceed with proving error estimates the velocity.
Theorem 6.5. Let c and u be the solutions of eqs.
Proof. Proof of eq. (62a): , and the coercivity of a h eq. (32) yields Using eq. (29) and employing Young's inequality for some > 0, It follows from eq. (31) and Young's inequality that Consider now I 3 : The first term on the right side of eq. (66) can be bounded as follows using Lipschitz continuity of µ, the generalized Hölder's inequality for integrals and sums, and Young's inequality: Next we bound I 32 . By Lipschitz continuity of µ, eq. (21), eq. (34d), generalized Hölder's inequality, and Young's inequality, We bound I 33 by using the assumption on K given in eq. (6): Again by the Lipschitz property of µ and Hölder's inequality, Combining eqs. (66) to (70), Since f s is Lipschitz continuous in c, with Lipschitz constant L s f , and recalling eq. (29), Since f d and µ are Lipschitz continuous in c, with Lipschitz continuity constants L d f and µ L , respectively, Combining the above bounds for I 1 to I 5 with eq. (63), letting = C a 16 (C a is the coercivity constant), multiplying by 2∆t, summing from 1 to n, noting that ζ 0 u = 0, and applying eq. (41), we obtain: Next, using eq. (34d) and eq. (45b), where the constant C > 0 depends on µ L , κ * , K * , γ j , L s f but is independent of h and ∆t. Equation (62a) follows by eqs. (54a) to (54c), and assumption eq. (61).

Proof of eq. (62b):
Let v s h = 0,v h = 0, and q h = 0 in eq. (48). Then On the other hand, letting v h = 0 and q s h = 0 in eq. (48), we have From [30, Lemma 3.2], since u sn − u sn h ∈ H div (Ω s ), there exists w ∈ H div (Ω d ) such that ∇ · w = 0 in Ω d , w · n = 0 on Γ d and w · n d = (γ(u sn ) − u sn h ) · n d on Γ I . With this choice of w, and using eqs. (34a) to (34c) we observe that Adding eqs. (75) and (76), and recalling eqs. (19b) and (19c), the definition ofΠ V , and eq. (34b), we obtain This leads us to consider v h = ζ n u + Π V w as test function in eq. (74). Using eqs. (7) and (77) we find that We will bound A 1 to A 4 by a series of Cauchy-Schwarz, Hölder's, triangle, and Young's inequalities together with the properties of µ and κ. First, Following the proof of eq. (69), while as the proof of eq. (73), Combining the bounds of A 1 to A 4 with eq. (78), choosing = 1/(8K * ), and applying another set of Young's inequalities, we obtain where C depends on the problem parameters κ * , K * , K * , µ L , L d f , and the regularity of u, p and c. Noting that from eq. (34d) and [30,Theorem 3.3], [17,Lemma 10], Then eqs. (34d), (60), (79) and (80) imply Equation (62b) is now a consequence of the assumptions on the regularity of the exact solution and eq. (61).
The following is a straightforward consequence of Theorem 6.5.
Corollary 6.6. Let u n and u n h be as defined in Theorem 6.5. Then for all n ≥ 1, Before moving on to the next section, we note another consequence of eq. (62) that will prove useful in analysis later on.
γ j , and the regularity of u 0 , u, c 0 , and c, but is independent of h, n and ∆t.
Proof. The proof is the same as in [16, Lemma 1]. Using eqs. (62a) and (62b), we have Combining the above estimate with eqs. (28) and (34d), for each K ∈ T , we have The result eq. (82) follows by using the assumption on ∆t and by taking the maximum over K ∈ T .
Remark 6.8. The restrictions on the polynomial degree and the time step are not necessary if we assume that |D(u)| ≤D for some positive constantD as in [43, 2.12]. It is also possible to avoid these restrictions by using an approach involving a cutoff operator on the velocity solution, as in [49,50], if one is interested in lower order approximations.
Remark 6.9. Compatibility, as defined in [22], can be achieved by choosing k c = k f − 1 [16]. However, the requirement that k c ≥ dim − 1 implies k f ≥ dim. Therefore, when dim = 2, our theory supports compatibility only for k f ≥ 2 and for k f ≥ 3 when dim = 3.

Error estimate for the pressure
In this section, we briefly discuss the a priori error estimate for the pressure approximation.
Proof. Setting q h = 0 in the error equation in Theorem 6.3, we obtain: By eq. (31) and Young's inequality, and using that |||·||| v and |||·||| v are equivalent on V h , we have Following the proof of eq. (71), we can show that As in eqs. (72) and (73), we find Using Cauchy-Schwarz and triangle inequalities, and eq. (34d), Therefore, combining (85)-(90), dividing both sides by |||v h ||| v , taking the supremum over v h ∈ V h , and using Theorem 4.2, we obtain Squaring both sides, multiplying by (c * inf ) −2 ∆t, summing from 1 to n, using eqs. (41a) and (45b), stability of Π V , and the regularity assumptions on u s , u d , and ∇p d yields Therefore, the result follows by eqs. (54a) to (54c) under the assumptions on the exact solution given in Theorem 6.5.
Remark 6.11. An immediate consequence Lemma 6.10 and Theorem 6.5 is This loss of ∆t is due to ∆t n m=1 d t ζ m u 2 Ω s in eq. (84). However, an improved estimate can be obtained by bounding this term as follows: testing eq. (48) with v h = d t ζ m u , multiplying by ∆t, summing from m = 1 to n, employing a summation-by-parts formula on the terms on the right hand side that are contained in a h (·; ·, d t ζ m u ) to transfer the discrete time derivative on d t ζ m u to the other terms, and assuming that the exact solution is sufficiently smooth in time, leads to We do not provide the details of this proof here, but instead refer to [18, p.42].

Existence and uniqueness of the concentration solution
In this section, we will prove existence and uniqueness of the discrete concentration solution c n h ∈ C h to eq. (40). First observe that assumption eq. (4b) on D implies that for and together with eq. (22) that where C depends on v 1,∞,K . Therefore by eqs. (82), (91) and (92), there exists a constant D max > 0 that depends on d and the upper bound in eq. (82) such that With eq. (93), the following coercivity result can be proved following the same steps as the proofs of [16, Lemmas 2 and 3].
Theorem 6.12 (Coercivity of B tr h (u n h ; w h , w h )). There exists a constant β 0 > 0 such that if β tr > β tr 0 , then for all w h ∈ C h , where C tr > 0 is a constant that depends on d, D min , and the upper bound in eq. (82). Now that we have coercivity, we proceed with the existence and stability proof for the discrete concentration.
Proof. Let w h = c n h in eq. (40). From the algebraic inequality (a−b)a ≥ 1 2 (a 2 −b 2 ), eq. (94), and the assumption that 0 ≤ c I ≤ 1 a.e., we have Then for sufficiently small ∆t, where C depends on φ * , φ * , d, D and the regularity of the solution but is independent of h and ∆t.
Proof. Setting w h = ζ n c in Lemma 6.14, using the inequality a(a − b) ≥ a 2 −b 2 2 , and Theorem 6.12, Using eq. (5a), the Cauchy-Schwarz inequality, Young's inequality with constant γ, and eq. (53), Again by eq. (5a), this time using Taylor's theorem in integral form, and applying Young's inequality, The following series of inequalities is dedicated to finding an upper bound for I 3 . By definition of B tr h , We will bound B a h and B d h separately, starting with B a h . Noting that ξ n c −ξ n c vanishes on facets, we have by eq. (17), I 311 ≤ C u n h 0,∞,Ω ξ n c Ω ∇ζ n c Ω ≤ Ch kc c n kc,Ω |||ζ n c ||| c .
By the triangle inequality and eq. (53), we immediately have
Note that this solution satisfies all the interface conditions and that ∇ · u s = 0 in Ω s . We present our numerical results for a wide range of values for κ and µ: κ = µ = 1; κ = 10 3 , µ = 10 −6 ; κ = 1, µ = 10 −6 ; and κ = 10 −3 , µ = 10 −6 . Since we are primarily interested in the spatial error, to minimize the temporal error as much as possible, we use the third order backward differentiation formulae (BDF3) as time stepping method even though the sequential algorithm 1 is only first order accurate in time. We choose ∆t = 0.1h k f /(k f + 1) and present errors and rates of convergence using k f = 2, k c = 1 in Tables 1 to 3 and using k f = 3, k c = 2 in Tables 4 to 6. Tables 1 and 2 for k f = 2, k c = 1, and Tables 4 and 5 for k f = 3, k c = 2 show that in the Stokes and Darcy regions u h and p h both converge optimally in the L 2 -norm with orders k f + 1 and k f , respectively. This is Table 2. Errors and rates of convergence at final time T = 0.1 for u h and p h in the Darcy region Ω d for the test case in Section 7.1 using k f = 2, k c = k f − 1, and BDF3 time stepping with ∆t = 0.1h 2 /3.  consistent with our theoretical convergence rate in Corollary 6.6 that predicts at least suboptimal rates for the velocity. The right most columns in these tables demonstrate pointwise mass conservation. Furthermore, even though the magnitude of the pressure error changes dramatically as we change the values of κ and µ, there is no significant change in the velocity errors. This is more pronounced in the case where both the permeability and the viscosity are small (10 −3 and 10 −6 ). This confirms that the velocity error bounds in Theorem 6.5 and Corollary 6.6 are independent of the pressure error.
We furthermore observe from Table 3 and Table 6 that c h converges optimally in the L 2 -norm with order k c + 1. This supports our result eq. (117) that shows at least suboptimal convergence.

Example 2
We now consider the fully coupled problem by incorporating the influence of the velocity solution on the dispersion/diffusion tensor and the dependence of the viscosity on the concentration solution. The source terms and boundary conditions for the Stokes/Darcy-transport problem (1) are chosen such that the exact solution is given by eq. (118). We define the diffusion dispersion tensor in Ω and the viscosity according to where we remark the the viscosity is defined as the quarter-power mixing rule [38] with where µ 0 = 0.9, µ 1 = 1.3.
We use k f = 3, k c = 2, and BDF3 time stepping with ∆t = 0.1h 3 /4. We present numerical results for κ = 10 3 , 1, 10 −3 . We observe from Tables 7 and 8 that when µ is changed from a constant to a concentration dependent function the rate of convergence reduces from k f + 1 to a value between k f and k f + 1. This is due to our choice k c = k f − 1 to achieve compatibility and is consistent with our a priori estimates eqs. (62a) and (62b) in the energy norm which imply that the rate of convergence of the velocity approximation is polluted by the concentration approximation. Indeed, from Table 9 we observe that c h converges in the L 2 -norm with order k f . Therefore, for the velocity we expect an order of at least k f − 1 in the energy norm and k f in the L 2 -norm. From the right most columns in Tables 7 and 8 we observe that the discretization is exactly mass conserving.

Example 3
In this final example, we simulate a more realistic problem similar to [17,Section 6.2] in which the permeability field in the Darcy region is highly heterogeneous. For this, let the boundary of the Stokes region be partitioned We impose the following boundary conditions: u = (x 2 (3/2 − x 2 )/5, 0) on Γ s 1 , (−2µε(u) + pI)n = 0 on Γ s 2 , u · n = 0 and (−2µε(u) + pI) t = 0 on Γ s 3 , u · n = 0 on Γ d 1 , p = −0.05 on Γ d 2 .  The first boundary condition on the left boundary Γ s 1 of Ω s imposes a parabolic velocity profile. We set the permeability to κ = 700(1 + 0.5(sin(10πx 1 ) cos(20πx 2 2 ) + cos 2 (6.4πx 1 ) sin(9.2πx 2 ))) + 100, a plot of which is given in Figure 2. The viscosity is defined by the quarter-power mixing rule as in eq. (119). The other parameters are set as µ = 0.1, α = 0.5, k f = 3, h = 1/80, ∆t = 10 −3 , T = 15, and the source/sink terms are set to zero. In the Darcy region Ω d the porosity is set to φ = 0.4. The dispersion/diffusion tensor is  We compute the solution using BDF3 time stepping. Figure 3 shows the computed velocity field after one time step and at the final time. In the Darcy region Ω d , the flow field avoids areas with low permeability as expected. Figure 4 shows the pressure contours at various times which demonstrates the effect of the concentration on the pressure around the plume of contaminants, especially in the Stokes region. Figure 5 presents the plume of contaminant spreading through the surface water region and infiltrating the porous medium. We plot the solution 6 different instances in time. The contaminant plume stays compact while in the surface water region. Once it reaches the subsurface region it spreads out following a path dictated by the heterogeneous permeability structure of the porous medium.

Conclusions
In this paper, we introduced and analyzed a fully discrete sequential method for the fully coupled Stokes/Darcytransport problem. The spatial discretization uses the HDG method which is higher-order accurate, strongly mass conservative, and compatible. We remark that the analysis also easily extends to the EDG-HDG method considered in [16,17]. The sequential method discussed in the article linearizes the problem by time-lagging the concentration and decoupling the Stokes/Darcy and transport subproblems. We proved well-posedness and obtained a priori estimates in the energy norm. Finally, we presented numerical results demonstrating mass conservation and robustness with respect to varying permeability and optimal convergence in the L 2 -norm for one-way coupling.