A strongly conservative hybridizable discontinuous Galerkin method for the coupled time-dependent Navier–Stokes and Darcy problem

We present a strongly conservative and pressure-robust hybridizable discontinuous Galerkin method for the coupled time-dependent Navier–Stokes and Darcy problem. We show existence and uniqueness of a solution and present an optimal a priori error analysis for the fully discrete problem when using Backward Euler time stepping. The theoretical results are veriﬁed by numerical examples.

of sources and sinks, pointwise divergence-free on the elements [28].In particular, we consider an HDG method [16] that is based on the HDG method for the Navier-Stokes equations [35] and a hybridized formulation of the mixed form of the Darcy problem [2], although nonconforming formulations based on other strongly conservative discretizations, for example, [15,21,31,41], are possible.
Previously, we proved pressure-robustness of strongly conservative HDG methods for the Stokes/Darcy [9] and stationary Navier-Stokes/Darcy [6] problems, leading to a priori error estimates for the velocity that do not depend on the best approximation of the pressure scaled by the inverse of the viscosity (see [27,32] for a review of other pressure-robust discretizations).Using Backward Euler time stepping we now show existence and uniqueness of a solution and derive an a priori error estimate to the fully-discrete time-dependent problem.Compared to previous work on the time-dependent Navier-Stokes/Darcy problem [7, 10-12, 25, 43], the novel contributions of this work is therefore the introduction and analysis of a strongly conservative HDG discretization and an a priori error estimate for the velocity that is independent of pressure.
The remainder of this paper is organized as follows.We present the time-dependent Navier-Stokes/Darcy problem in Section 2 and its HDG discretization in Section 3. Consistency and well-posedness of the discrete problem are shown in Section 4 while a priori error estimates are proven in Section 5. We end this paper with numerical examples in Section 6 and conclusions in Section 7.

The Navier-Stokes and Darcy problem
We consider the time-dependent incompressible Navier-Stokes equations coupled to the Darcy equations on a polyhedral domain Ω in R dim , dim = 2, 3, and on the time interval  = (0,  ).The domain is partitioned into two non-overlapping subdomains Ω  and Ω  such that Ω = Ω  ∪ Ω  , Ω  ∩ Ω  = ∅, and Γ  := Ω  ∩ Ω  .The boundary of the domain Ω and the interface Γ  are assumed to be Lipschitz polyhedral.We define Γ  and Γ  to be the exterior boundaries of Ω  and Ω  , respectively.We partition Γ  := Γ   ∪ Γ   , with Γ   ∩ Γ   = ∅ and |Γ   | > 0 and |Γ   | > 0, and denote the outward unit normal on Γ  to Ω  ( = , ) by .The Navier-Stokes equations are given by where   : Ω  × → R dim is the velocity in Ω  ,   : Ω  × → R dim is the pressure in Ω  , () = 1 2 (∇ +(∇)  ),  > 0 is the constant fluid viscosity, and   : Ω  ×  → R dim is a body force.In Ω  the Darcy equations are given by: where   : Ω  ×  → R dim is the fluid velocity in Ω  ,   : Ω  ×  → R is the piezometric head in Ω  , and  > 0 is the permeability constant.The Navier-Stokes equations are coupled to the Darcy equations by the following interface conditions −2((  )) where  is the unit normal vector on Γ  pointing from Ω  to Ω  , ()  :=  − ( • ) is the tangential component of a vector , and  > 0 is an experimentally determined dimensionless constant.Note that equation (3a) ensures continuity of the normal component of the velocity across the interface, equation (3b) is the Beavers-Joseph-Saffman law [4,38], and equation (3c) is a balance of forces.We assume the following initial and boundary conditions: = 0 on Γ  × , (4b) where  0 : Ω  → R dim is a solenoidal initial velocity field.We close this section by introducing  : Ω ×  → R dim and  : Ω ×  → R to be the functions such that | Ω  =   and | Ω  =   for  = , .
3. The HDG method

Notation
Let  = , .We denote by   ℎ = {} a conforming triangulation of Ω  of shape-regular simplices .Our goal is to present a strongly conservative and pressure robust HDG discretization.A necessary requirement for our discretization to achieve this is that the velocity field is globally (div; Ω)-conforming [35,36].For this reason we assume that  ℎ =   ℎ ∪   ℎ is a matching simplicial mesh, i.e.,   ℎ and   ℎ match at the interface.Without this assumption, technical complications may arise to obtain a globally (div; Ω)-conforming velocity field.We denote by ℎ  the diameter of  ∈  ℎ and define the meshsize as ℎ := max ∈ ℎ ℎ  .A face  is an interior face if for two elements  + and  − in  ℎ ,  =  + ∩  − , and a boundary face if  ∈  lies on the boundary Ω.The set of all facets in Ω and Ω  are denoted by, respectively, ℱ ℎ and ℱ  ℎ , while the set of all facets on the interface Γ  is denoted by ℱ  ℎ .The set of all facets on Γ  are denoted by ℱ , ℎ while the set of all facets interior to Ω  are denoted by ℱ , ℎ .The sets of facets on Γ   and Γ   are denoted by, respectively, ℱ , ℎ and ℱ , ℎ .By Γ 0 and Γ  0 we denote the union of facets in Ω and Ω  .The outward unit normal vector on  for any element  ∈   ℎ is denoted by   .On the interface Γ  ,  =   = −  .We will drop the superscript  if the definition of the outward unit normal vector is clear.
Denoting by   () the space of polynomials of total degree  on a domain , we define the following finite element spaces for the velocity approximation: For notational purposes, we write  ℎ = ( ℎ , vℎ ) ∈  ℎ :=  ℎ × Xℎ and   ℎ = (  ℎ , vℎ ) ∈   ℎ :=   ℎ × Xℎ .Furthermore, for the pressure approximation we define the finite element spaces We write For scalar functions  and , we define Similar notation is used for vector-and matrix-valued functions.

The semi-discrete problem
An HDG method for the stationary Navier-Stokes and Darcy problem was proposed in [6].Its extension to the time-dependent problem is given by: Let  ,0 ℎ ∈   ℎ ∩ (div; Ω  ) be the initial condition for the velocity in The different forms are defined as: where  > 0 is a penalty parameter and where   ℎ is the linear part of  ℎ .For the velocity-pressure coupling we have, for  = , , the forms:

The fully-discrete problem
Using backward Euler time-stepping, and lagging the convective velocity in the nonlinear term, we obtain the following linear implicit discretization: Let  ,0 ℎ ∈   ℎ ∩ (div; Ω  ) be the initial condition for the velocity in Ω  such that ∇ •  ,0 ℎ = 0 pointwise on each Remark 3.1.As observed previously in [6] for the stationary Navier-Stokes and Darcy problem, the velocity solution to equation (7) satisfies the following properties: (i) it is exactly divergence-free on elements in Ω On Ω  and Ω  , we define the function spaces On Ω, we then define The trace space of   on facets in Γ  0 is denoted by X.If  ∈   , we denote its trace by ū :=   () where   :   → X is the trace operator restricting functions in   to Γ  0 .Similarly, the trace space of   on facets Γ  0 is denoted by Q ,    :   → Q is the trace operator, and if  ∈   , then q :=    () ∈ Q .
Using the notation  :=  × X and  :=  × Q × Q , we define As in [6], we define the following norms on the extended function spaces: where Here  •  is the usual jump operator and ‖‖ [8].Finally, we will also require the following two norms on the pressure in Ω  : The following inequalities will be used in the remainder of this paper (see [42], Eq. (5.5), [22], Thm.4.4 and Prop.4.5, and [17], Lem.1.46): where   ,   ,  , , and   are positive constants independent of ℎ and ∆.
For  ℎ we have: Due to the use of different function spaces, the inf-sup condition equation (9a) is different from that proven in [6].We therefore prove equation (9a) in Appendix A. Equation (9b) is proven in Lemma 3 of [6].For   ℎ ,   , and   , we have that for all ,  ∈ (ℎ), where    > 0 is a constant independent of ℎ and ∆.For  ℎ ∈  ℎ we have where the first inequality holds for  large enough and where    > 0 is a constant independent of ℎ and ∆.A direct consequence of equations (10) and (11) is that where    := max(   ,  −1 ,  −1/2 ) > 0 and    := min(   ,  −1 ,  −1/2 ) > 0 are constants independent of ℎ and ∆, and where equation (12b) holds for  large enough.

Existence and uniqueness
We start this section with some auxiliary results.Lemma 4.2.For  , ℎ and  , ℎ that satisfy equation (7), there exists a   > 0, independent of ℎ and ∆, such that Proof.We will prove equation ( 16) in three dimensions only noting that the proof in two dimensions is similar.
To ease notation we will drop the "time" superscript .The proof follows the proof of Lemma 2.1 from [34] with modifications made to take into account Brezzi-Douglas-Marini (BDM) elements and HDG facet functions.The local degrees of freedom for the BDM element are ( [5], Prop.2.3.2): where 3 ⊂  −2 () and since   ℎ − p ℎ ∈   (), we obtain from equation (18) that Setting now   ℎ = 0 and  ℎ = 0 in equation (7), and after integration by parts, we find for all  ℎ ∈   ℎ that: Choose  ℎ =  ℎ , with  ℎ defined in equation (18).By equations ( 19), (20), and the definition of To find out more about ‖ ℎ ‖ Ω  , let us define the norm ‖•‖ 0,ℎ for functions in   ℎ ∩ (div; Ω  ): Consider now a single element  and denote by ℱ  the set of faces of .In an approach similar to that used in the proof of Lemma 4.4 from [33], we have: where the first line on the right hand side is by using the degrees of freedom equation ( 17), the second by definition of  ℎ given by equation (18), and the last is by the Cauchy-Schwarz inequality.Therefore, after summing equation ( 23) over all  in   ℎ : The result follows after combining this with equation ( 21).
An immediate consequence of equation (8c) and Lemma 4.2 is that if  , ℎ and  , ℎ satisfy equation ( 7), then for 1 where The following result, which was shown in Theorem 5.2 from [12], will be used to prove the next lemma: there exists a constant  > 0, independent of ℎ and ∆, such that where ̃︀ and p, ℎ be (part of ) the solution to equation (7).There exists a constant   > 0, independent of ℎ and ∆, such that for all  ≥ 1 Proof.For ease of notation we will drop the "time" superscript .Then, note that Since   ℎ is a solution to equation (7), by Remark 3.1 we know that ∇•  ℎ = 0 and Therefore, using equation ( 25), Next, using equation (8e) and Lemma 4.2, we note that The result follows by combining equations ( 24) and ( 27)- (29).
For the remainder of this section we define Let us remark that by equation (8d ) −1 for all   ℎ ∈   ℎ .This result, in combination with equation (15), is used in the following lemma to prove a well-posedness result.Proof.Consider equation (7) for the solution at time level  +1 which we write here as: Given  , ℎ ∈   ℎ we remark that, by equation (8d) with  = 2 and equation ( 15), Furthermore, by equations (14b) and (8b), we obtain the following boundedness result: where )︁ .
We are now ready to prove a bound on  +1 ℎ .
We end this section by proving existence and uniqueness for all time levels under a suitable data assumption.
A bound for ||| +1 ℎ |||  is given by Lemma 4.6 and the data assumption equation ( 54).Together with equation ( 38) we obtain Note that the data assumption equation (54) implies that The result follows from equations ( 56) and (57).
By a triangle inequality and properties of the interpolant Π  and projection Π , we obtain the following velocity error estimate that is independent of the pressure.

Numerical examples
We implement the fully discrete HDG method equation (7) in Netgen/NGSolve [39,40].For all examples we choose the penalty parameter as  = 8 2 (see [1,37]), where  is the polynomial degree in the approximation spaces.

Rates of convergence
In this section we verify the rates of convergence by the method of manufactured solutions.For this we consider the domains Ω  = (0, 1) × (0, 0.5) and Ω  = (0, 1) × (−0.5, 0).The interface is given by Γ To construct a manufactured solution, we consider the following inhomogeneous boundary conditions and modified interface conditions: where   ,   ,   ,   ,   ,   , and   , and the functions   and   in equations (1a) and (2b) are chosen such that the exact solution is given by: The initial condition for the velocity is set by first solving the stationary Stokes-Darcy problem with the above boundary/interface conditions and functions   and   .In our simulations we choose  = 10 −4 and  = 1.We consider polynomial degrees  = 1 (corresponding to approximating the cell pressure by piecewise constants and the other unknowns by piecewise linear polynomials) and  = 2 (in which the cell pressure is approximated by piecewise linears and the other unknowns by piecewise quadratic polynomials).We compare results obtained by choosing  = 10 −1 ,  = 10 −3 , and  = 10 −5 .Let us define   :=  −  ℎ and, similar to [22], From Corollary 5.2 we expect that, for smooth enough solutions, ‖  ‖  = (ℎ  + ∆).The spatial rate of convergence is indeed observed in Table 1 (to obtain these results we chose our time step as ∆ = 0.8ℎ +1 and set  = (0, 0.1)).Table 1 also lists the  2 -norm of   and   :=  −  ℎ .For the velocity we observe that ‖  ‖ Ω = (ℎ +1 ) for  = 10 −1 and ‖  ‖ Ω ≈ (ℎ +1/2 ) for  = 10 −5 .For  = 10 −3 we have that ‖  ‖ Ω lies between (ℎ +1/2 ) and (ℎ +1 ), depending on whether  = 1 or  = 2.The slower convergence in the  2 -norm for  = 10 −5 is not surprising; the flow problem is advection dominated and analysis of HDG methods for the scalar advection equation reveals a priori error estimates for the solution to be (ℎ +1/2 ), see Lemma 4.8 from [42].We furthermore observe optimal rates of convergence for the pressure: We next consider the temporal rates of convergence.For this we consider a fine mesh with 9508 cells and set  = 2 and  = (0, 1).In Table 2 we vary the time step and present the errors and rates of convergence.All errors are (∆).
Finally, let us remark that despite our analysis holding only under the small data assumption (see Eq. ( 54)), we are nevertheless able to compute the solution for very small values of viscosity.From Tables 1 and 2 we even observe that the variation in ‖  ‖  for the different values of  is small, despite the upper bound in Corollary 5.2 depending on  and  −1 .

Surface/subsurface flow with nonuniform permeability field
In this example we consider surface/subsurface flow.For this example we divide the domain Ω = (0, 1) × (−0.5, 0.5) into two subdomains Ω  and Ω  .We consider a case where the interface Γ  = Ω  ∩Ω  is not horizontal (see Fig. 1a).Furthermore, let Γ   = {︀  ∈ Γ  :  2 = −0.5 }︀ , and Γ   = Γ  ∖Γ   .We then impose the following boundary conditions: and set   = 0 and   = 0. We consider both  = 10 −1 and  = 10 −3 together with  = 0.5, and choose the permeability to be piecewise constant such that  −1  = 10 − with  ∈ [2, 6] a random number that is chosen differently in each element of the mesh in Ω  .(The analysis presented in this paper assumes a constant permeability, but noting that 0 <  min ≤ () ≤  max the analysis is easily extended to this situation.)A plot of the permeability is given in Figure 1b.To set the initial condition for the velocity in Ω  we solve the stationary Stokes-Darcy problem.We compute the solution on a mesh consisting of 91 720 elements, using  = 2, a time step of ∆ = 0.01, and on the time interval  = (0, 10).Plots of the velocity and pressure fields at different time levels are shown Table 1.Errors and spatial rates of convergence for a manufactured solution (see Sect. 6.1).Results are for  = 1 and  = 2 with parameters  = 10 −4 ,  = 1, and  ∈ {10 −1 , 10 −3 , 10 −5 }.
Here   =  −  ℎ and   =  −  ℎ .The rate of convergence is denoted by . in Figures 2 and 3, both for  = 10 −1 and  = 10 −3 .The velocity fields at  = 0 and  = 10 for both values of viscosity are similar: flow in Ω  away from the interface is more or less horizontal while in Ω  flow finds its way through the permeability maze in the direction of negative pressure gradient.At  = 5.2 (when the inflow magnitude of the velocity is close to its minimum), the behavior of the velocity fields when  = 10 −1 and  = 10 −3 are significantly different: when  = 10 −1 the velocity field is similar to that at  = 0 and  = 10, but when  = 10 −3 we obtain a large area of circulation.The pressure fields are similar for the two values of viscosity and follow a more or less linear profile in Ω  .Pressure variations in Ω  are small.

Conclusions
We presented a strongly conservative HDG method for the coupled time-dependent Navier-Stokes and Darcy problem.Existence and uniqueness of a solution to the fully discrete problem were proven assuming a small data assumption.We furthermore determined a pressure-independent a priori error estimate for the discrete velocity.This estimate is optimal in space in the combined discrete  1 -norm on Ω  and (div)-norm on Ω  , and optimal in time.Our analysis is supported by numerical examples.Table 2. Errors and temporal rates of convergence for a manufactured solution (see Sect.

Appendix B. Useful inequalities
Let  be a sufficiently smooth function.Using Taylor's theorem in integral form, it is shown in Lemma 7.67 from [26] that A minor modification of the proof of equation (B.1) leads to: We also have, by the fundamental theorem of Calculus and the Cauchy-Schwarz inequality, that