Error analysis of fully discrete mixed finite element data assimilation schemes for the Navier-Stokes equations

# Error analysis of fully discrete mixed finite element data assimilation schemes for the Navier-Stokes equations

Bosco García-Archilla Departamento de Matemática Aplicada II, Universidad de Sevilla, Sevilla, Spain. Research is supported by Spanish MINECO under grant MTM2015-65608-P (bosco@esi.us.es).    Julia Novo Departamento de Matemáticas, Universidad Autónoma de Madrid, Spain. Research is supported by Spanish MINECO under grant MTM2016-78995-P (AEI/FEDER, UE) and VA024P17 (Junta de Castilla y Leon, ES) cofinanced by FEDER funds (julia.novo@uam.es).
###### Abstract

In this paper we consider fully discrete approximations with inf-sup stable mixed finite element methods in space to approximate the Navier-Stokes equations. A continuous downscaling data assimilation algorithm is analyzed in which measurements on a coarse scale are given represented by different types of interpolation operators. For the time discretization an implicit Euler scheme, an implicit and a semi-implicit second order backward differentiation formula are considered. Uniform in time error estimates are obtained for all the methods for the error between the fully discrete approximation and the reference solution corresponding to the measurements. For the spatial discretization we consider both the Galerkin method and the Galerkin method with grad-div stabilization. For the last scheme error bounds in which the constants do not depend on inverse powers of the viscosity are obtained.

AMS subject classifications. 35Q30, 65M12, 65M15, 65M20, 65M60, 65M70,
76B75.
Keywords. data assimilation, downscaling, Navier-Stokes equations, uniform-in-time error estimates, fully discrete schemes, mixed finite elements methods.

## 1 Introduction

Data assimilation refers to a class of techniques that combine experimental data and simulation in order to obtain better predictions in a physical system. There is a vast literature on data assimilation methods (see e.g., [4], [12], [33], [35], [39], and the references therein). One of these techniques is nudging in which a penalty term is added with the aim of driving the approximate solution towards coarse mesh observations of the data. In a recent work [6], a new approach, known as continuous data assimilation, is introduced for a large class of dissipative partial differential equations, including Rayleigh-Bénard convection [17], the planetary geostrophic ocean dynamics model [18], etc. (see also references therein). Continuous data assimilation has also been used in numerical studies, for example, with the Chafee-Infante reaction-diffusion equation, the Kuramoto-Sivashinsky equation (in the context of feedback control) [36], Rayleigh-Bénard convection equations [3], [16], and the Navier-Stokes equations [24], [27]. However, there is still less numerical analysis of this technique. The present work concerns with the numerical analysis of continuous data assimilation for the Navier-Stokes equations for fully discrete schemes with inf-sup stable mixed finite element methods (MFE) in space.

We consider the Navier-Stokes equations (NSE)

 ∂tu−νΔu+(u⋅∇)u+∇p =f in  (0,T]×Ω, ∇⋅u =0 in  (0,T]×Ω, (1)

in a bounded domain , . In (1), is the velocity field, the kinematic pressure, the kinematic viscosity coefficient, and represents the accelerations due to external body forces acting on the fluid. The Navier-Stokes equations (1) must be complemented with boundary conditions. For simplicity, we only consider homogeneous Dirichlet boundary conditions on .

As in [37] we consider given coarse spatial mesh measurements, corresponding to a solution of (1), observed at a coarse spatial mesh. We assume that the measurements are continuous in time and error-free and we denote by the operator used for interpolating these measurements, where denotes the resolution of the coarse spatial mesh. Since no initial condition for is available one cannot simulate equation (1) directly. To overcome this difficulty it was suggested in [6] to consider instead a solution  of the following system

 ∂tv−νΔv+(v⋅∇)v+∇~p = f−β(IH(v)−IH(u)), in  (0,T]×Ω, ∇⋅v = 0, in  (0,T]×Ω, (2)

where is the nudging parameter.

In [23] the continuous in time data assimilation algorithm is analyzed and two different methods are considered: the Galerkin method and the Galerkin method with grad-div stabilization. In this paper we extend the results in [23] to the fully discrete case. For the time discretization of equation (1), we consider the fully implicit backward Euler method and the second order backward differentiation formula (BDF2), both in the fully implicit and semi-implictit cases. For the spatial discretization we consider inf-sup stable mixed finite elements. As in [23] we consider both the Galerkin method and the Galerkin method with grad-div stabilization. Although grad-div stabilization was originally proposed in [19] to improve the conservation of mass in finite element methods, it was observed in the simulation of turbulent flows in [32], [40] that grad-div stabilization has the effect of producing stable (non-oscillating) simulations.

For the three time discretization methods that we consider and the two different spatial discretizations (Galerkin method with or without grad-div stabilization) we prove uniform-in-time error estimates for the approximation of the unknown reference solution, , corresponding to the coarse-mesh measurement . As in [23], for the Galerkin method without stabilization, the spatial error bounds we prove are optimal, in the sense that the rate of convergence is that of the best interpolant. In the case where grad-div stabilization is added, as in [13], [14], we get error bounds where the error constants do not depend on inverse powers of the viscosity parameter . This fact is of importance in many applications where viscosity is orders of magnitude smaller than the velocity (or with large Reynolds number).

We now comment on the literature on numerical methods for (1). In [37], a semidiscrete postprocessed Galerkin spectral method for the two-dimensional Navier-Stokes equations is studied. Under suitable conditions on the nudging parameter , the coarse mesh size , and the degrees of freedom in the spectral method, uniform-in-time error estimates are obtained for the error between the numerical approximation to  and . The use of a postprocessing technique introduced in [21] [22], gives a method with a higher convergence rate than the standard spectral Galerkin method. A fully-discrete method for the spatial discretization in [37] is analyzed in [30], where the (implicit and semi-implicit) backward Euler method is used for time discretization and uniform-in-time error estimates are obtained with the same convergence rate in space as in [37].

Other related works are  [34] and [38]. In [38] they only analyze linear problems and, for the proof of the results on the Navier-Stokes equations they present, they refer to [34] with some differences that they point out. They also present a wide collection of numerical experiments. In [34], the authors consider fully discrete approximations to equation (1) where the spatial discretization is performed with a MFE Galerkin method with grad-div stabilization. A semi-implicit BDF2 scheme in time is analyzed in [34], and, as in [30], [37], [23] and the present paper uniform-in-time error bounds are obtained. In the present paper, apart from the semi-implicit BDF2 scheme of [34] we also analyze the implicit Euler and the implicit BDF2 schemes. Respect to the spatial errors we obtain the same results as in [23]. More precisely, comparing with [34], we remove the dependence on inverse powers on on the error constants of the Galerkin method with grad-div stabilization. Also, for the standard Galerkin method, although with constants depending on inverse powers on , we get a rate of convergence for the method in space one unit larger than the method in [34]. These results are sharp in space, as it can be checked in the numerical experiments of [23]. This means that the Galerkin method with grad-div stabilization has a rate of convergence in the norm of the velocity using polynomials of degree and that error constants are independent on  for the grad-div stabilized method, and dependent in the case of the standard method. The analysis in [34] is restricted to being an interpolant for non smooth functions (Clément, Scott-Zhang, etc), since explicit use is made of bounds (24) and (25), which are not valid for nodal (Lagrange) interpolation. In the present paper, as in [23], we prove error bounds both for the case in which is an interpolant for non smooth functions but also for the case in which is a standard Lagrange interpolant. To our knowledge reference [23] and the present paper are the only references in the literature where such kind of bounds are proved.

It is important to mention that, compared with [34][38], [30] and [37], and as in [23], we do not need to assume an upper bound on the nudging parameter . The authors of [34] had observed (see [34, Remark 3.8]) that the upper bound on  they required in the analysis does not hold in the numerical experiments. This fact is corroborated by the numerical experiments in [23] that show numerical evidence of the advantage of increasing the value of the nudging parameter well above the upper bound in required in previous works in the literature. The analysis in the present paper, as that in [23], does not demand any upper bound on the nudging parameter .

For the error analysis of the fully discrete method with the implicit Euler scheme we do not need to assume any restriction on the size of the time step (Theorem 3.6 below). For the case of the implicit BDF2 and semi-implicit BDF2 (Theorems 3.9 and 3.12 below) we only need to assume is smaller than a constant (that depends on norms of the theoretical exact solution).

The rest of the paper is as follows. In Section 2 we state some preliminaries and notation. In Section 3 we analyze the fully discrete schemes. First of all, some general results are stated and proved and then the error analysis of the implicit Euler method, the implicit BDF2 method and the semi-implicit BDF2 schemes is carried out. In Section 4 some numerical experiments are shown.

## 2 Preliminaries and Notation

We denote by  the standard Sobolev space of functions defined on the domain with distributional derivatives of order up to in . By  we denote the standard seminorm, and, following [11], for  we will define the norm  by

 ∥f∥ps,p,D=s∑j=0|D|p(j−s)d|f|pj,p,D,

where denotes the Lebesgue measure of . Observe that is scale invariant. If is not a positive integer, is defined by interpolation [1]. In the case , we have . As it is customary, will be endowed with the product norm and, since no confusion will arise, it will also be denoted by . When , we will use to denote the space . By we denote the closure in of the set of infinitely differentiable functions with compact support in . The inner product of or will be denoted by and the corresponding norm by . The norm of the dual space of is denoted by .

We will use the following Sobolev’s inequality [1]: For , let and  be such that . Then, there exists a positive scale invariant constant such that

 ∥v∥Lq′(Ω)≤cs|Ω|sd−1p+1q′∥v∥Ws,p(Ω),1q′≥1q,∀v∈Ws,p(Ω). (3)

If the above relation is valid for . A similar inequality holds for vector-valued functions.

We will denote by and the Hilbert spaces , , endowed with the inner product of and respectively.

The following interpolation inequalities will also be used (see, e.g., [11, formula (6.7)] and [20, Exercise II.2.9])

 (4)

(where, for simplicity, by enlarging the constants if necessary, we may take the constant  in (4) equal to  in (3) for ) and Agmon’s inequality

 (5)

The case is a direct consequence of [2, Theorem 3.9]. For , a proof for domains of class  can be found in [11, Lemma 4.10].

We will make use of Poincaré’s inequality,

 ∥v∥0≤cP|Ω|1/d∥∇v∥0,∀v∈H10(Ω), (6)

where the constant  can be taken . If we denote by , then from (6) it follows that

 ∥v∥1≤(^cP)1/2∥∇v∥0,∀v∈H10(Ω). (7)

The constants , , and  in the inequalities above are all scale-invariant. This will also be the case of all constants in the present paper unless explicitly stated otherwise. We will also use the well-known property, see [31, Lemma 3.179]

 ∥∇⋅v∥0≤∥∇v∥0,v∈H10(Ω)d. (8)

Let , be a family of partitions of suitable domains , where denotes the maximum diameter of the elements , and are the mappings from the reference simplex onto . We shall assume that the partitions are shape-regular and quasi-uniform. Let , we consider the finite-element spaces

 Sh,r = {χh∈C(¯¯¯¯Ωh)∣∣χh|τhj∘ϕhj∈Pr−1(τ0)}⊂H1(Ωh), S0h,r = Sh,r∩H10(Ωh),

where denotes the space of polynomials of degree at most on . For , stands for the space of piecewise constants.

When has polygonal or polyhedral boundary and mappings  from the reference simplex are affine. When has a smooth boundary, and for the purpose of analysis, we will assume that exactly matches , as it is done for example in [10], [41]. At a price of a more involved analysis, though, the discrepancies between and can also be included in the analysis (see, e.g., [5], [42]).

We shall denote by the MFE pair known as Hood–Taylor elements [8, 44] when , where

 Xh,r=(S0h,r)d,Qh,r−1=Sh,r−1∩L2(Ωh)/R,r≥3,

and, when , the MFE pair known as the mini-element [9] where , and . Here, is spanned by the bubble functions , , defined by , if  and 0 elsewhere, where denote the barycentric coordinates of . All these elements satisfy a uniform inf-sup condition (see [8]), that is, for a constant independent of the mesh size the following inequality holds

 infqh∈Qh,r−1supvh∈Xh,r(qh,∇⋅vh)∥vh∥1∥qh∥L2/R≥βis. (9)

To approximate the velocity we consider the discrete divergence-free space

 Vh,r=Xh,r∩{χh∈H10(Ωh)d∣(qh,∇⋅χh)=0∀qh∈Qh,r−1}.

For each fixed time , notice that the solution of (1) is also the solution of a Stokes problem with right-hand side . We will denote by its MFE approximation, solution of

 ν(∇sh,∇φh)−(qh,∇⋅φh) = ν(∇u,∇φh)−(p,∇⋅φh) (10) = (f−ut−(u⋅∇u),φh)∀φh∈Xh,r, (∇⋅sh,ψh) = 0∀ψh∈Qh,r−1.

Notice that is the discrete Stokes projection of the solution of (1) (see [28]), and it satisfies satisfies that for all ,

 ν(∇Sh(u),∇φh)=ν(∇u,∇φh)−(p,∇⋅φh)=(f−ut−(u⋅∇)u,φh).

The following bound holds (see e.g., [29]):

 ∥u−sh∥0+h∥u−sh∥1≤CNj(u,p)hj,1≤j≤r, (11)

where here and in the sequel, for and  we denote

 (12)

and, when and  depend on

 ¯¯¯¯¯Nj(v,q)=supτ≥0Nr(v(τ),q(τ)), Mj(v)=supτ≥0∥v(τ)∥j, Mj(p)=supτ≥0∥p(τ)∥Hj/R. (13)

Assuming that is of class , with , and using standard duality arguments and (11), one obtains

 ∥u−sh∥−s≤CNr(u,p)hr+s,0≤s≤min(r−2,1). (14)

We also have the following bounds [23, Lemma 3.7]

 ∥sh∥∞ ≤D0((∥u∥d−2∥u∥2)1/2+(N1(u,p)Nd−1(u,p))1/2) (15) ∥∇sh∥L2d/(d−1) ≤D0(N1(u,p)N2(u,p))1/2, (16)

where the constant  does not depend on . And, arguing as in [23, Lemma 3.7] one can also prove

 ∥∇sh∥∞ ≤D0((∥u∥d−1∥u∥3)1/2+(N2(u,p)Nd(u,p))1/2). (17)

We also consider a modified Stokes projection that was introduced in [13] and that we denote by satisfying

 ν(∇smh,∇φh)=(f−ut−(u⋅∇)u−∇p,φh),∀φh∈Vh,r, (18)

and the following error bound, see [13]:

 ∥u−smh∥0+h∥u−smh∥1≤C∥u∥jhj,1≤j≤r. (19)

From [10], we also have

 ∥∇(u−smh)∥∞≤C∥∇u∥∞, (20)

where does not depend on . Also we have the following bounds [23, Lemma 3.8]

 ∥smh∥∞ ≤D1(∥u∥d−2∥u∥2)1/2, (21) ∥∇smh∥L2d/(d−1) ≤D1(∥u∥1∥u∥2)1/2, (22)

where the constant  is independent of .

We will denote by the projection of the pressure onto . It is well-known that

 ∥p−πhp∥0≤Chj−1∥p∥Hj−1/R,1≤j≤r. (23)

We will assume that the interpolation operator is stable in , that is,

 ∥IHu∥0≤c0∥u∥0,∀u∈L2(Ω)d, (24)

and that it satisfies the following approximation property,

 ∥u−IHu∥0≤cIH∥∇u∥0,∀u∈H10(Ω)d. (25)

The Bernardi–Girault [7], Girault–Lions [25], or the Scott–Zhang [43] interpolation operators satisfy (25) and (24). Notice that the interpolation can be on piecewise constants.

Finally, we will denote by the Lagrange interpolant of a continuous function . In Subsection 3.4 we consider the case in which for which bounds (24) and (25) do not hold.

## 3 Fully discrete schemes

We consider the following method to approximate (1). For we define satisfying for all

 (dtu(n)h,φh)+ν(∇u(n)h, ∇φh)+bh(u(n)h,u(n)h,φh)+μ(∇⋅u(n)h,∇⋅φh) +(∇p(n)h,φh) =(f(n),φh)−β(IH(u(n)h)−IH(u(n)),IHφh), (∇⋅u(n)h,ψh) =0. (26)

In (3) is a finite difference approximation to the time derivative at time , where is the time step. Also, is a stabilization parameter that can be zero in case we do not stabilize the divergence or different from zero in case we add grad-div stabilization and is defined in the following way

 bh(uh,vh,φh)=((uh⋅∇)vh,φh)+12(∇⋅(uh)vh,φh),∀uh,vh,φh∈Xh,r.

It is straightforward to verify that enjoys the skew-symmetry property

 bh(u,v,w)=−bh(u,w,v)∀u,v,w∈H10(Ω)d. (27)

Let us observe that taking from (3) we get

 (dtu(n)h,φh)+ν(∇u(n)h,∇ φh)+bh(u(n)h,u(n)h,φh)+μ(∇⋅u(n)h,∇⋅φh)= (f(n),φh)−β(IH(u(n)h)−IH(u(n)),IHφh). (28)

For the analysis below, we introduce the values  and , defined as follows

 ¯¯¯μ={0,\rm if~{}μ=0,1,\rm otherwise,¯¯¯k={0,\rm if~{}μ=0,1/μ,\rm otherwise. (29)
###### Lemma 3.1

The following bound holds for , and ,

 |bh(^vh,vh,eh)−bh(^wh,wh,eh)|≤ 1δ^L(wh,ε)∥eh∥20 +δ(∥∇w∥∞2∥^eh∥20+ε4∥∇⋅^eh∥20),

where , and

 ^L(wh,ε)=(∥∇wh∥∞2+∥wh∥2∞4ε), (30)
###### Proof:

Adding and applying the skew-symmetry property (27) we have

 |bh(^vh,vh,eh)−bh(^wh,wh ,eh)|=|bh(^eh,wh,eh)|≤∥∇wh∥∞∥^eh∥0∥eh∥0 +12∥∇⋅^eh∥0∥wh∥∞∥eh∥0 ≤∥∇wh∥∞2δ∥eh∥20+δ∥∇wh∥∞2∥^eh∥20 +∥wh∥2∞4δϵ∥eh∥20+δϵ4∥∇⋅^eh∥20.

The proof of the following lemma is similar to the proof of [23, Lemma 3.1].

###### Lemma 3.2

Let be the finite element approximation defined in (3) and let , , in be sequences satisfying

 (dtw(n)h,φh)+ν(∇w(n)h, ∇φh)+bh(w(n)h,w(n)h,φh)+μ(∇⋅w(n)h,∇⋅φh)= (f(n),φh)+(τ(n)h,φh)+¯¯¯μ(θ(n)h,∇⋅φh), (31)

Assume that the quantity  defined in (39), below, is bounded. Then, if and satisfies condition (44), below, the following bounds hold for ,

 (dte(n)h,e(n)h)+γ2∥e(n)h∥20+μ2∥∇⋅e(n)h∥20≤¯¯¯k∥θ(n)h∥20+β2c20∥u(n)−w(n)h∥20 +((1−¯¯¯μ)^cPν+¯μ2L)∥τ(n)h∥2−1+¯¯¯μ, (32)

where, and  are defined in (29), and is defined in (47) below.

###### Proof:

Subtracting (3.2) from (3) we get the error equation

 (dte(n)h,φh)+ν(∇e(n)h,∇φh)+β(IHe(n)h,IHφh) (33) +bh(u(n)h,u(n)h,φh)−bh(w(n)h,w(n)h,φh)+μ(∇⋅e(n)h,∇⋅φh)= β(IHu(n)−IHw(n)h,IHφh)+(τ(n)h,φh)+¯¯¯μ(θ(n)h,∇⋅φh),

for all Taking in (33) we get

 (dte(n)h,e(n)h)+ν∥∇e(n)h∥20+β∥IHe(n)h∥20+ μ∥∇⋅e(n)h∥20 (34) ≤|bh(u(n)h,u(n)h,e(n)h−bh(w(n)h,w(n)h,eh (n))|+β|(IHu(n)−IHw(n)h,IHe(n)h)| +|(τ(n)h,e(n)h)|+|¯¯¯μ(θ(n)h,∇⋅e(n)h)|.

We will bound the terms on the right-hand side of (34). For the nonlinear term and the truncation errors we argue differently depending on whether or .

If , applying Lemma 3.1 with , , and , we have

 |bh(u(n)h,u(n)h,e(n)h)−bh(w(n)h,w(n)h,e(n)h)| ≤(^L(w(n)h,ν)+12∥∇w(n)h∥∞)∥e(n)h∥20+ν4∥∇⋅e(n)h∥20, ≤(^L(w(n)h,ν)+12∥∇w(n)h∥∞)∥e(n)h∥20+ν4∥∇e(n)h∥20 (35)

where in the last inequality we have applied (8). For the term when , using (7) we get

 |(τ(n)h,e(n)h)| ≤∥τ(n)h∥−1∥e(n)h∥1≤(^cP)1/2∥τ(n)h∥−1∥∇e(n)h∥0 ≤^cPν∥τ(n)h∥2−1+ν4∥∇e(n)h∥20. (36)

When , we bound the nonlinear term by applying Lemma 3.1 with , , and , that is

 |bh(u(n)h,u(n)h,e(n)h)−bh(w(n)h,w(n)h,e(n)h)| ≤(^L(w(n)h,μ)+12∥∇w(n)h∥∞)∥e(n)h∥20+μ4∥∇⋅e(n)h∥20, (37)

In the sequel we denote

 L = maxn≥0(^L(w(n)h,ν)+12∥∇w(n)h∥∞)ifμ=0, (38) L = (39)

Observe that in the case , the left-hand side of (3) can be bounded by , and, in the case  the left-hand side of (3) is bounded by .

Next, we bound the last two terms of the right-hand side of (33) when . We have

 |(τ(n)h,e(n)h)|+|¯¯¯μ(θ(n)h,∇⋅e(n)h)|≤ 12L∥τ(n)h∥20+L2∥e(n)h∥20 +¯¯¯k∥θ(n)h∥20+μ4∥∇⋅e(n)h∥20, (40)

where is defined in (29).

For the second term on the right-hand side of (34) applying (24) we get

 β|(IHu(n)−IHw(n)h,IHeh)| ≤ βc0∥u(n)−w(n)h∥0∥IHe(n)h∥0 (41) ≤ β2c20∥u(n)−w(n)h∥20+β2∥IHe(n)h∥20.

Inserting (3), (3), (3), (3) and (41) into (34) we get

 (dte(n)h,e(n)h)+(1+¯¯¯μ)ν2∥∇e(n)h∥20+β2∥IHe(n)h∥20+μ2∥∇⋅e(n)h∥20 (42) ≤L∥eh∥20+¯¯¯k∥θ(n)h∥20+β2c20∥u(n)−w(n)h∥20+((1−¯¯¯μ)^cPν+¯¯¯μ2L)∥τ(n)h∥2−1+¯¯¯μ.

Now we bound

 L∥e(n)h∥20≤2L∥IHe(n)h∥20+2L∥(I−IH)e(n)h∥20.

Assuming that and taking into account that we get

 (dte(