## Analysis of continuous data assimilation with large (or even infinite) nudging parameters

Amanda E. Diegel , ∗ Xuejian Li , 1 Leo G. Rebholz 2

## Abstract

This paper considers continuous data assimilation (CDA) in partial differential equation (PDE) discretizations where nudging parameters can be taken arbitrarily large. We prove that solutions are long-time optimally accurate for such parameters for the heat and NavierStokes equations (using implicit time stepping methods), with error bounds that do not grow as the nudging parameter gets large. Existing theoretical results either prove optimal accuracy but with the error scaled by the nudging parameter, or suboptimal accuracy that is independent of it. The key idea to the improved analysis is to decompose the error based on a weighted inner product that incorporates the (symmetric by construction) nudging term, and prove that the projection error from this weighted inner product is optimal and independent of the nudging parameter. We apply the idea to BDF2 - finite element discretizations of the heat equation and Navier-Stokes equations to show that with CDA, they will admit optimal long-time accurate solutions independent of the nudging parameter, for nudging parameters large enough. Several numerical tests are given for the heat equation, fluid transport equation, Navier-Stokes, and Cahn-Hilliard that illustrate the theory.

Keywords: continuous data assimilation, long time accuracy, finite element method, projection properties

## 1. Introduction

Data assimilation is a broad topic that generally refers to incorporating observations into a given dynamical system to better estimate some relevant quantities of interest. One common use of data assimilation is to generate an accurate initial condition for the dynamical system so that its future state can be accurately predicted. Data assimilation was originally made popular in the weather forecasting community [27, 2, 35], and various categories of data assimilation methods exist, such as Kalman filtering and particle filtering [31, 17, 3,

∗ Department of Mathematics and Statistics, Mississippi State University, MS, 39762; partially supported by NSF grant DMS 2110768.

1 School of Mathematical and Statistical Sciences, Clemson University, Clemson, SC, 29364, xuejial@clemson.edu; partially supported by NSF grant DMS 2152623.

2 School of Mathematical and Statistical Sciences, Clemson University, Clemson, SC, 29364, rebholz@clemson.edu; partially supported by NSF grant DMS 2152623.

19], 3D/4D variational methods [15, 14, 40, 37], and nudging methods [44, 5, 4]. In the past decade, the continuous data assimilation (CDA) method developed by Azouani, Olson, and Titi [7] has rapidly become popular due to its robustness, effectiveness, and strong mathematical foundation.

CDA works by adding a nudging term (also called feedback control or penalty term) that assimilates observations into dynamical systems and attempts to drive the assimilated state to the true state. In general, CDA applied to a time dependent PDE

$$& u _ { t } + F ( u ) = f, \\ & u ( 0 ) = u _ { 0 },$$

is written as (suppressing boundary conditions for the moment)

$$\frac { \partial w } { \partial t } + F ( w ) + \mu ( I _ { H } w - I _ { H } u ) = f, \quad \sum _ { \lambda \lambda } \lambda ^ { ( 1 ) }$$

$$w ( 0 ) = 0.$$

Here, µ I ( H w -I H u ) is the nudging term, u is the partially observed true solution, I H is an interpolation or observation operator associated with the measurement location with H representing the sparsity of observations, µ is the nudging parameter emphasizing the accuracy of observations, and F is an operator specific to the PDE at hand. We note that an initial condition of 0 is chosen for w , although long-time accuracy is typically achieved in many CDA analyses considering arbitrary initial conditions. Of course, the more accurate the initial condition, the quicker the convergence, in general. While CDA is similar to traditional nudging methods, a fundamental difference is the use of the interpolation operator which is a minor modification but has a critical impact. Researchers have proved in a variety of settings and systems that by using sparse observations and appropriate I H , applying CDA leads to an assimilated solution that converges to the true solution (up to optimal discretization error) exponentially fast in time, independent of the initial condition (see e.g.[12, 13, 16, 18, 8, 1, 9, 30] and references therein).

In addition to the strong mathematical foundation, another reason CDA is attractive is its straightforward implementation. The CDA computation is carried out the same way as solving the associated evolutionary PDE, only with extra attention to the nudging term that can be easily handled in most cases. Studies have been performed using the CDA implementation based on different interpolants I H and discretization schemes that make implementation simple (e.g. only altering the diagonal of the system matrix at each time step) [41, 43, 33, 29, 20, 21], and maintain the expected exponential convergence of discrete CDA to an optimally accurate solution.

The goal of this paper is to prove optimal long-time convergence of discrete CDA with no upper restriction on the nudging parameter µ (i.e. no scaling of the error with µ ). There have been several studies that have shown excellent results with CDA with nudging parameter of 1000 and larger applied to Navier-Stokes and Cahn-Hilliard solvers where implicit time stepping is used [41, 22, 16, 25]; some of these works also test varying µ and show no deterioration of accuracy as µ gets larger. However, to date, no theoretical results exist to

explain these good results; current known results for CDA that are optimal either depend on µ or have an upper bound restriction on µ . Allowing µ to be very large can bring convenience of CDA implementation by allowing for a direct Dirichlet enforcement assimilation. Very large µ is numerically equivalent in the linear systems to a direct enforcement of w x ( i ) = u x ( i ) at all measurement nodes x i associated with the determination of I H u and I H w . Such an implementation has been mentioned in [36], but without any proof of accuracy. In this work, we provide a rigorous analysis to support this result in the context of CDA for evolutionary PDEs. We first propose projection operators tailored for use with CDA, and prove an optimal estimate for its approximation error, independent of the nudging parameter µ . With these new projections, we are able to prove that CDA with arbitrarily large parameters find optimal long-time accurate solutions with respect to the L 2 -norm with error bounds independent of µ for BDF2 - finite element discretizations of the heat equation and incompressible NavierStokes equations. We note that the tools and techniques developed and used for the heat and Navier-Stokes equations are extendable to many dissipative PDEs.

The paper is organized as follows. In Section 2, we introduce notation and provide mathematical preliminaries for a smooth analysis to follow. In Section 3, we state the new CDA Poisson and CDA Stokes projections and analyze their properties. In Section 4, we present the CDA convergence analysis for the discretized heat equation. In Section 5, we focus on a CDA convergence analysis of discretized Navier-Stokes equations. In Section 6, we provide several numerical tests that support the theory. Finally, we draw conclusions in Section 7.

## 2. Preliminaries

Throughout this paper, C will denote a generic constant, possibly changing at each instance, that is independent of mesh size h , time step ∆ , nudging parameter t µ , and problem data. We consider Ω to be a generic polygonal or polyhedral domain with smooth boundary ∂ Ω. Let H k (Ω) denote the Sobolev space W k, 2 (Ω) and L 2 (Ω) = W 0 2 , (Ω). We use ( · , · ) to denote the L 2 inner product that induces the L 2 norm ∥ · ∥ , while all other inner products or norms will be appropriately subscripted. Define the following function spaces:

$$\begin{array} { c c c } \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ & H _ { 0 } ^ { 1 } \left ( \Omega \right ) \coloneqq \{ v \in H ^ { 1 } \left ( \Omega \right ) \colon v = 0 & \text{on} \ \partial \Omega \}, \\ & X \coloneqq \left ( H _ { 0 } ^ { 1 } \left ( \Omega \right ) \right ) ^ { d }, & d = 2, 3, \\ & Q \coloneqq \{ v \in L ^ { 2 } ( \Omega ) \ \colon \int _ { \Omega } v d x = 0 \}, \\ & V \coloneqq \{ v \in X \ \colon \left ( \nabla \cdot v, q \right ) = 0 & \forall q \in Q \}. \end{array}$$

Recall that the Poincar´ inequality holds on e X : there exists a constant C P dependent only on the domain Ω satisfying

$$\| v \| \leq C _ { P } \| \nabla v \|,$$

for all v ∈ X .

## 2.1. Finite element preliminaries

We consider finite element (FE) methods for the spatial discretization of the CDA algorithm. Denote by τ h (Ω) a shape-regular and conforming triangulation of the interested domain Ω. We use Y h = H 1 0 (Ω) ∩ P k ( τ h (Ω)) as the FE discretization space for the heat equation. Let X h × Q h ⊂ X × Q be a pair of inf-sup stable FE spaces for Navier-Stokes equation, such as X h = X ∩ P k ( τ h (Ω)) and Q h = Q ∩ P k -1 ( τ h (Ω)) as Taylor-Hood or Scott-Vogelius elements.

## 2.2. BDF2 preliminaries

The time stepping method used in our analysis and most of our numerical tests is BDF2, and we give here some preliminary results for BDF2 analysis. We will utitlize the G-norm [23]: ∥ x ∥ 2 G := ( x, Gx ) R 2 , where G = [ 1 2 -1 -1 5 2 ] . It is easy to show that for x ∈ ( L 2 (Ω)) , the 2 G-norm is equivalent to the ( L 2 (Ω)) 2 norm:

$$C _ { \ell } \| x \| _ { G } \leq \| x \| \leq C _ { u } \| x \| _ { G }.$$

Here, C ℓ = 3 -2 √ 2 and C u = 3 + 2 √ 2. For v i ∈ L , i 2 = n -1 , n, n +1, we have the identity [23]:

$$\left ( \frac { 1 } { 2 } \left ( 3 v ^ { n + 1 } - 4 v ^ { n } + v ^ { n - 1 } \right ), v ^ { n + 1 } \right ) & = \frac { 1 } { 2 } \left ( \| [ v ^ { n + 1 } ; v ^ { n } ] \| _ { G } ^ { 2 } - \| [ v ^ { n } ; v ^ { n - 1 } ] \| _ { G } ^ { 2 } \right ) \\ & + \frac { 1 } { 4 } \| v ^ { n + 1 } - 2 v ^ { n } + v ^ { n - 1 } \| ^ { 2 }. \\ \wedge \wedge \, \sim \, \sim \,. \quad. \ \.$$

## 2.3. CDA preliminaries

Denote by τ H (Ω) a coarse mesh of Ω, whose nodes are locations of the measurements or observables of the true solution. Assume h &lt;&lt; H and the nodes of mesh τ H (Ω) are a subset of the vertices of the mesh τ h (Ω) and that if I h is the nodal interpolant operator on τ h , then I H h I u = I H u for any u ∈ H 1 (Ω). Also, I H must satisfy the usual CDA estimates: there exists a constant C I independent of H satisfying

$$\| I _ { H } v - v \| \leq C _ { I } H \| \nabla v \|, \quad \forall v \in H ^ { 1 } ( \Omega ),$$

$$\| I _ { H } v \| \leq C _ { I } \| v \|, \quad \forall v \in H ^ { 1 } ( \Omega ).$$

We note that if I H were the nodal interpolant, then C I in (5) could be large due to an inverse dependence on h for v ∈ X h or Y h and so this is not an ideal interpolant choice in most cases. Instead, algebraic nudging [41], the Scott-Zhang interpolant, or the L 2 projection onto the coarse space could be used [20].

We also need the following lemma for the long-time analysis that is proven in [32].

Lemma 2.1. Suppose constants r &gt; 1 and B ≥ 0 . If a real number sequence { a n } n ≥ 0 satisfies

$$\ r a _ { n + 1 } \leq a _ { n } + B,$$

then

$$a _ { n + 1 } \leq a _ { 0 } \left ( \frac { 1 } { r } \right ) ^ { n + 1 } + \frac { B } { r - 1 }.$$

## 3. CDA Poisson and CDA Stokes projections

In this section we will present and analyze Poisson and Stokes interpolants tailored to use with CDA. The estimates derived in this section are the keys that allow for the optimal long-time accuracy results that are independent of the nudging parameter to hold in the sections that follow.

## 3.1. CDA Poisson projection

↦

The CDA Poisson projection is defined as follows: Given u ∈ H 1 0 (Ω) and I H u , find u h := P u h → u h in Y h satisfying

$$\kappa ( \nabla u _ { h }, \nabla v _ { h } ) + \mu ( I _ { H } u _ { h }, I _ { H } v _ { h } ) = \kappa ( \nabla u, \nabla v _ { h } ) + \mu ( I _ { H } u, I _ { H } v _ { h } ), \quad \forall v _ { h } \in Y _ { h }, \quad \quad ( 9 )$$

where κ &gt; 0 is a positive constant.

Lemma 3.1. Given u ∈ H 1 0 (Ω) ∩ H k +1 (Ω) and µ ≥ 0 , P u h satisfies

$$\| \nabla ( u - P _ { h } u ) \| \leq C h ^ { k } \| u \| _ { H ^ { k + 1 } ( \Omega ) },$$

$$\| u - P _ { h } u \| & \leq C h ^ { k + 1 } \| u \| _ { H ^ { k + 1 } ( \Omega ) },$$

where C is independent of u , h , and µ .

Proof. We first rewrite the equation (9) as follows:

$$( \nabla ( u - u _ { h } ), \nabla v _ { h } ) + \frac { \mu } { \kappa } ( I _ { H } ( u - u _ { h } ), I _ { H } v _ { h } ) = 0, \quad \forall v _ { h } \in Y _ { h }.$$

Set e = u -u h = ( u -I h u ) + ( I h u -u h ) =: η + ϕ , where I h u ∈ Y h is the nodal interpolant of u . Properties of I h and I H imply that

$$I _ { H } \eta = I _ { H } ( u - I _ { h } u ) = I _ { H } u - I _ { H } I _ { h } u = I _ { H } u - I _ { H } u = 0.$$

Using this, after choosing v h = ϕ in (12), then yields

$$\| \nabla \phi \| ^ { 2 } + \frac { \mu } { \kappa } \| I _ { H } \phi \| ^ { 2 } = ( \nabla \eta, \nabla \phi ) \leq \frac { 1 } { 2 } \| \nabla \phi \| ^ { 2 } + \frac { 1 } { 2 } \| \nabla \eta \| ^ { 2 },$$

thanks to Young's inequality. Reducing (13) gives

$$\| \nabla \phi \| ^ { 2 } + \frac { 2 \mu } { \kappa } \| I _ { H } \phi \| ^ { 2 } \leq \| \nabla \eta \| ^ { 2 }.$$

With the triangle inequality and standard approximation theory, we get

$$\| \nabla e \| ^ { 2 } + \frac { \mu } { \kappa } \| I _ { H } e \| ^ { 2 } \leq 4 \| \nabla ( u - I _ { h } u ) \| ^ { 2 } \leq C h ^ { 2 k } \| u \| _ { H ^ { k + 1 } ( \Omega ) } ^ { 2 },$$

where C is independent of h and µ . This proves the first result of the lemma.

For the L 2 estimate, we begin by defining ˆ ξ ∈ C 2 (Ω) ∩ C 0 (Ω) satisfying

$$- \Delta \xi = e, \quad \text{on } \Omega.$$

Let ψ be any harmonic function (e.g. a harmonic polynomial [6]) that satisfies ψ x ( i ) = ˆ ( ξ x i ) for all the measurement points x i . Then we have for ξ := ˆ ξ -ψ that

$$- \Delta \xi = e, \ I _ { H } \xi = 0.$$

Testing the first ξ equation with e yields

$$\| e \| ^ { 2 } = ( \nabla \xi, \nabla e ).$$

From (12) with v h = I h ξ , we have that

$$( \nabla e, \nabla I _ { h } \xi ) + \mu ( I _ { H } e, I _ { H } I _ { h } \xi ) = 0,$$

and so combining this with (16) gives

$$\| e \| ^ { 2 } = ( \nabla e, \nabla ( \xi - I _ { h } \xi ) ) - \mu ( I _ { H } e, I _ { H } I _ { h } \xi ).$$

Since the measurement nodes { x i } M i =1 are nodes on τ h and τ H ⊂ τ h , I H h I ξ = I H ξ = 0. Hence we have the estimate

$$\| e \| & = ( \nabla e, \nabla ( \xi - I _ { h } \xi ) ) \\ & \leq \| \nabla e \| \| \nabla ( \xi - I _ { h } \xi ) \| \\ & \leq C h \| \nabla e \| \| \Delta \xi \| \\ & = C h \| \nabla e \| \| e \|,$$

where C is an interpolation constant and independent of µ . This reduces

$$\| e \| \leq C h \| \nabla e \|,$$

and using the bound (15) for ∥∇ ∥ e proves the second result of the lemma.

## 3.2. CDA Stokes projection

↦

We now consider the CDA Stokes projection, which is defined as follows: Given u ∈ X and I H u , find u h := P S h u → u h in V h = { v ∈ X h : ( ∇· v, q ) = 0 ∀ q ∈ Q h } satisfying

$$\nu ( \nabla u _ { h }, \nabla v _ { h } ) + \mu ( I _ { H } u _ { h }, I _ { H } v _ { h } ) = \nu ( \nabla u, \nabla v ) + \mu ( I _ { H } u, I _ { H } v _ { h } ), \quad \forall v _ { h } \in V _ { h }. \quad \quad ( 1 7 )$$

At first, this projection may appear to be very similar to the CDA Poisson projection with Y h replaced by V h . However, there is a subtle difference regarding the fact that I H does not necessarily map into the (discretely divergence free space) V h . Hence, there is a more delicate analysis needed to prove that this projection is optimal in L 2 and H 1 .

For simplicity, we will consider Scott-Vogelius elements with X h = P k ( τ h ) ∩ X and Q h = P k -1 ∩ L 2 0 (Ω) with k = 2, but the same idea can be applied for k &gt; 2 in 2D and for k ≥ 3 in 3D (for lower order elements than these, the following construction would not work). The key idea is to alter the nodal interpolation I h similar to what was done in [26] as follows: I S h will be a nodal interpolant at all the vertices of τ h (and therefore also at all the measurement nodes so that I H h I S = I H ), but the value at the midedges is defined so that the total flux through any element boundary is 0. Hence on each edge e , where { x i } 3 i =1 are the nodes on e (i=1 and 3 correspond to vertices, i=2 to the midedge), and ⃗ n is the unit outward normal, ψ i ( x ) is the finite element nodal basis that belongs to support point x i , ⃗ e j is a d dimension unit vector with j th component 1 and others 0, ψ j i ( x ) is a d dimension vector with j th component as ψ i ( x ) and others 0, we make the following correction to the midedge value:

$$I _ { h } ^ { S } u ( x _ { 2 } ) = u ( x _ { 2 } ) + \frac { \int _ { e } u ( s ) \cdot \vec { n } d s - \sum _ { i, j } \int _ { e } ( ( \vec { e } _ { j } \cdot u ( x _ { i } ) ) \psi _ { i } ^ { j } ( s ) ) \cdot \vec { n } d s } { \sum _ { j = 1 } ^ { d } \int _ { e } \psi _ { 2 } ^ { j } ( s ) \cdot \vec { n } d s } \vec { n }.$$

With this definition of I S h , one can calculate that the total flux through each element E is zero, and thus by the divergence theorem that ∫ E ∇ · I S h u dx = 0. With this, since ∇· X h ⊂ Q h , we get that ∥∇ · I S h u ∥ L 2 ( E ) = 0 and thus ∥∇ · I S h u ∥ = 0. Therefore, it holds that I S h : H 1 → V h .

In addition to being divergence-free, it is proven in [26] that I S h has the following property.

Lemma 3.2. For u ∈ L ∞ (Ω) , the interpolant I S h satisfies

$$\max _ { x \in e } | u ( x ) - I _ { h } ^ { S } u ( x ) | \leq C h ^ { k + 1 }.$$

Since all nodes of τ h are on edges of the mesh, we immediately get that the max error at any node is Ch k +1 . This can be used with standard approximation theory to get that

$$\| u ( x ) - I _ { h } ^ { S } ( u ) \| \leq C h ^ { k + 1 } \| u \| _ { H ^ { k + 1 } ( \Omega ) },$$

$$\| \nabla \left ( u - I _ { h } ^ { S } u \right ) \| \leq C h ^ { k } \| u \| _ { H ^ { k + 1 } ( \Omega ) }.$$

With these approximation properties and with I S h mapping into V h , we have the following result for the CDA-Stokes projection:

Lemma 3.3. Given u ∈ X ∩ H k +1 (Ω) and µ ≥ 0 , the Stokes projection P S h u satisfies

$$\| \nabla ( u - P _ { h } ^ { S } u ) \| \leq C h ^ { k } \| u \| _ { H ^ { k + 1 } ( \Omega ) },$$

$$\| u - P _ { h } ^ { S } u \| \leq C h ^ { k + 1 } \| u \| _ { H ^ { k + 1 } ( \Omega ) },$$

where C is independent of u , h , and µ .

Proof. Follow the proof of Lemma 3.1 with I h replaced by I S h and Y h replaced by V h .

## 4. Analysis of large parameter CDA for BDF2-FEM discretization of heat equation

We now consider an analysis of CDA with large nudging parameters applied to the BDF2FEM discretization of the heat equation with form

$$u _ { t } - \kappa \Delta u & = f, \\ u ( \cdot, 0 ) & = u _ { 0 }, \\ u | _ { \partial \Omega } & = 0,$$

where κ &gt; 0 is the thermal diffusivity constant.

Introducing the bilinear form a ( · , · ) : H 1 0 × H 1 0 → R defined by a ϕ, v ( ) = ( κ ∇ ϕ, ∇ v ), the CDA algorithm we consider herein for the heat equation using a BDF2-FEM discretization is as follows:

Algorithm 4.1. Let ∆ t be an uniform time step and Y h be the FE space. Given any initials condition w 0 h ∈ L 2 (Ω) , w 1 h ∈ Y h , source f ∈ L ∞ (0 , ∞ ; L 2 (Ω)) , I H satisfying (5)(6), observations { I H u n +1 } n ≥ 1 ∈ L 2 (Ω) , and nudging parameter µ &gt; 0 , find { w n h } ∈ Y h for n = 1 2 3 , , , · · · , satisfying

$$\left ( \frac { 3 w _ { h } ^ { n + 1 } - 4 w _ { h } ^ { n } + w _ { h } ^ { n - 1 } } { 2 \Delta t }, v _ { h } \right ) + a ( w _ { h } ^ { n + 1 }, v _ { h } ) + \mu \left ( I _ { H } ( w _ { h } ^ { n + 1 } - u ^ { n + 1 } ), I _ { H } v _ { h } ) = ( f, v _ { h } ), \quad ( 2 5 )$$

for all v h ∈ Y h .

We now prove that this method has temporally and spatially optimal long-time accuracy with respect to the L 2 -norm for any µ &gt; 0.

Theorem 4.1. Let u ∈ L ∞ (0 , ∞ ; H k +1 (Ω)) be the solution of the heat equation (24) where u , u t tt , u ttt ∈ L ∞ (0 , ∞ ; L 2 (Ω)) , u 0 = u 0 ∈ L 2 (Ω) , u 1 ∈ H 1 0 (Ω) , and f ∈ L ∞ (0 , ∞ ; L 2 (Ω)) . Then, for any given H &gt; 0 and µ ≥ 1 2 C H 2 I 2 , the difference between the CDA solution of (25) and the true solution satisfies:

$$& \| w _ { h } ^ { n + 1 } - u ^ { n + 1 } \| ^ { 2 } \leq C \left ( \frac { 1 } { 1 + \lambda \Delta t } \right ) ^ { n + 1 } \left ( \| [ w _ { h } ^ { 1 } - u ^ { 1 } ; w _ { h } ^ { 0 } - u ^ { 0 } ] \| _ { G } ^ { 2 } + \frac { \kappa \Delta t } { 4 } \| \nabla ( w _ { h } ^ { 1 } - u ^ { 1 } ) \| ^ { 2 } \right ) \\ & \quad + \frac { C } { \lambda } \left ( \Delta t ^ { 4 } + h ^ { 2 k + 2 } \right ) \\. & \quad. & \quad. & \quad. & \quad \kappa c ^ { 2 } \quad \kappa c ^ { 2 } \cdot \quad. & \quad. & \quad.$$

for any n ≥ 1 . Here, λ = min ∆ { t -1 , κC 2 ℓ 2 C H 2 I 2 , κC 2 ℓ 4 C 2 P } and the constant C is independent of µ , h , and ∆ t .

Remark 4.2. In (4.1), the nudging parameter µ does not appear on any right side terms, which implies that there will be no adverse effect by considering arbitrarily large µ in the CDA algorithm.

Proof. Using Taylor's theorem, the solution of the heat equation (24) satisfies

$$\left ( \frac { 3 u ^ { n + 1 } - 4 u ^ { n } + u ^ { n - 1 } } { 2 \Delta t }, v _ { h } \right ) + a ( u ^ { n + 1 }, v _ { h } ) = ( f, v _ { h } ) + \frac { \Delta t ^ { 2 } } { 3 } ( u _ { t t t } ( t ^ { * } ), v _ { h } ), \ \forall v _ { h } \in Y _ { h }, \quad ( 2 7 )$$

where t ∗ ∈ [ t n -1 , t n +1 ]. Subtracting equation (25) from (27) gives the error equation

$$\left ( \frac { 3 e ^ { n + 1 } - 4 e ^ { n } + e ^ { n - 1 } } { 2 \Delta t }, v _ { h } \right ) & + a ( e ^ { n + 1 }, v _ { h } ) + \mu \left ( I _ { H } e ^ { n + 1 }, I _ { H } v _ { h } \right ) \\ & = \frac { \Delta t ^ { 2 } } { 3 } ( u _ { t t t } ( t ^ { * } ), v _ { h } ), \ \forall v _ { h } \in Y _ { h },$$

where e n = u n -w n h . We decompose e n as: e n = ( u n -P u h n ) + ( P u h n -w n h ) = η n + ϕ n . Then with the projection equation (9), we have that a η, v ( h ) + µ I ( H η, I H h v ) = 0, so we can rewrite (28) as

$$& \left ( \frac { 3 \phi ^ { n + 1 } - 4 \phi ^ { n } + \phi ^ { n - 1 } } { 2 \Delta t }, v _ { h } \right ) + a ( \phi ^ { n + 1 }, v _ { h } ) + \mu \left ( I _ { H } \phi ^ { n + 1 }, I _ { H } v _ { h } \right ) \\ & = - \left ( \frac { 3 \eta ^ { n + 1 } - 4 \eta ^ { n } + \eta ^ { n - 1 } } { 2 \Delta t }, v _ { h } \right ) + \frac { \Delta t ^ { 2 } } { 3 } ( u _ { t t t } ( t ^ { * } ), v _ { h } ).$$

Taking v h = ϕ n +1 in equation (29), using the BDF2 identity (4), and multiplying by 2∆ , t we have

$$& \| [ \phi ^ { n + 1 } ; \phi ^ { n } ] \| _ { G } ^ { 2 } - \| [ \phi ^ { n } ; \phi ^ { n - 1 } ] \| _ { G } ^ { 2 } + \frac { 1 } { 2 } \| \phi ^ { n + 1 } - 2 \phi ^ { n } + \phi ^ { n - 1 } \| ^ { 2 } + 2 \kappa \Delta t \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ & + 2 \mu \Delta t \| I _ { H } \phi ^ { n + 1 } \| ^ { 2 } = - \left ( 3 \eta ^ { n + 1 } - 4 \eta ^ { n } + \eta ^ { n - 1 }, \phi ^ { n + 1 } \right ) + \frac { 2 \Delta t ^ { 3 } } { 3 } ( u _ { t t t } ( t ^ { * } ), \phi ^ { n + 1 } ).$$

Following the analysis from [41] for the time derivative term and using the Cauchy-Schwarz and Young's inequalities on the right-hand side terms of (30), we obtain

$$\lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad\lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \ \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \ \alpha \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \alpha \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \quad \lambda \ \lambda
 \quad \quad \quad$$

Using the interpolation property (5), we provide the lower bound of the left-hand side terms

of equation (30) by

$$& \text{on} \varproj \ \mathfrak { v } \\ & 2 \kappa \Delta t \| \nabla \phi ^ { n + 1 } \| ^ { 2 } + 2 \mu \Delta t \| I _ { H } \phi ^ { n + 1 } \| ^ { 2 } = \kappa \Delta t \| \nabla \phi ^ { n + 1 } \| ^ { 2 } + 2 \mu \Delta t \| I _ { H } \phi ^ { n + 1 } \| ^ { 2 } + \kappa \Delta t \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ & \quad \geq \frac { \kappa \Delta t } { C _ { I } ^ { 2 } H ^ { 2 } } \| \phi ^ { n + 1 } - I _ { H } \phi ^ { n + 1 } \| ^ { 2 } + 2 \mu \Delta t \| I _ { H } \phi ^ { n + 1 } \| ^ { 2 } + \kappa \Delta t \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ & \quad \geq \beta \left ( \| \phi ^ { n + 1 } - I _ { H } \phi ^ { n + 1 } \| ^ { 2 } + \| I _ { H } \phi ^ { n + 1 } \| ^ { 2 } \right ) + \kappa \Delta t \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ & \quad \geq \frac { \beta } { 2 } \| \phi ^ { n + 1 } \| ^ { 2 } + \kappa \Delta t \| \nabla \phi ^ { n + 1 } \| ^ { 2 }, \\ & \quad \text{where} \ \mathfrak { B } = \min \{ 2 \mu \Delta t \ \stackrel { \kappa \Delta t } { \longleftarrow } \ \mathfrak { v } = \stackrel { \kappa \Delta t } { \longleftarrow } \ \mathfrak { h } \text{allowinv} \ n \text{ to the sufficient large} \ i \ e \ \ n > \frac { \kappa } { \kappa }$$

where β = min 2 { µ ∆ t, κ ∆ t C H 2 I 2 } = κ ∆ t C H 2 I 2 by allowing µ to be sufficient large, i.e., µ ≥ κ 2 C H 2 I 2 . Combining (30), (31), and (32) results in

$$\text{comming} \gtrsim \emptyset _ { J }, \, \infty _ { J }, \, \text{and} \, \infty _ { J } \text{couss} \, m \\ \| [ \phi ^ { n + 1 } ; \phi ^ { n } ] \| _ { G } ^ { 2 } & + \frac { \kappa \Delta t } { 4 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } + \frac { \kappa \Delta t } { 4 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } + \frac { \kappa \Delta t } { 2 C _ { I } ^ { 2 } H ^ { 2 } } \| \phi ^ { n + 1 } \| ^ { 2 } + \frac { \kappa \Delta t } { 4 } \| \nabla \phi ^ { n } \| ^ { 2 } \\ & \leq \| [ \phi ^ { n } ; \phi ^ { n - 1 } ] \| _ { G } ^ { 2 } + \frac { \kappa \Delta t } { 4 } \| \nabla \phi ^ { n } \| ^ { 2 } + C \kappa ^ { - 1 } \Delta t ^ { 5 } \| u _ { t t t } ( t ^ { * } ) \| ^ { 2 } \\ & + C \kappa ^ { - 1 } \Delta t \left ( \| \eta _ { t } \| _ { L ^ { \infty } ( 0, \infty ; L ^ { 2 } ) } ^ { 2 } + \int _ { t ^ { n - 1 } } ^ { t ^ { n + 1 } } \| \eta _ { t t } \| ^ { 2 } d t \right ), \\. \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots$$

where we added term κ ∆ t 4 ∥∇ ϕ n ∥ 2 on both sides of (33).

After using the Poincar´ inequality and norm equivalence (3), we obtain e

$$\text{ave using one to ensure equivalency and no import equivalence } \varphi, \, \text{we would} \\ \| [ \phi ^ { n + 1 } ; \phi ^ { n } ] \| _ { G } ^ { 2 } + \frac { \kappa \Delta t } { 4 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } + \lambda \Delta t \left ( \| [ \phi ^ { n + 1 } ; \phi ^ { n } ] \| _ { G } ^ { 2 } + \frac { \kappa \Delta t } { 4 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \right ) \leq \\ \| [ \phi ^ { n } ; \phi ^ { n - 1 } ] \| _ { G } ^ { 2 } + \frac { \kappa \Delta t } { 4 } \| \nabla \phi ^ { n } \| ^ { 2 } + C \kappa ^ { - 1 } \Delta t ^ { 5 } \| u _ { t t t } ( t ^ { * } ) \| ^ { 2 } \\ + C \kappa ^ { - 1 } \Delta t \left ( \| \eta _ { t } \| _ { L ^ { \infty } ( 0, \infty ; L ^ { 2 } ) } ^ { 2 } + \int _ { t ^ { n - 1 } } ^ { t ^ { n + 1 } } \| \eta _ { t t } \| ^ { 2 } d t \right ). \\ \sim _ { \sim 2 } \quad \sim _ { \sim 2 }$$

Here, λ = min ∆ { t -1 , κC 2 ℓ 2 C H 2 I 2 , κC 2 ℓ 4 C 2 P } . Using Lemma 2.1 and the estimate on η , we have

$$\| [ \phi ^ { n + 1 } ; \phi ^ { n } ] \| _ { G } ^ { 2 } + \frac { \kappa \Delta t } { 4 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } & \leq \left ( \frac { 1 } { 1 + \lambda \Delta t } \right ) ^ { n + 1 } \left ( \| [ \phi ^ { 1 } ; \phi ^ { 0 } ] \| _ { G } ^ { 2 } + \frac { \kappa \Delta t } { 4 } \| \nabla \phi ^ { 1 } \| ^ { 2 } \right ) \\ & \quad + \frac { C } { \lambda } \left ( \Delta t ^ { 4 } + h ^ { 2 k + 2 } \right ). \\ \vdots \quad \Pi \quad. \quad e \quad. \quad \dots \quad. \quad. \quad. \quad. \quad. \quad. \quad.. \quad........................................................................................................................................................................................................$$

Finally, the proof of Theorem 4.1 is completed by applying the triangle inequality and another application of the norm equivalence property (3).

## 5. Discrete CDA analysis for Navier-Stokes equations

We consider a BDF2-FEM algorithm with CDA for the NSE, and assume divergence-free elements ( X ,Q h h ) that also satisfy the discrete inf-sup condition. We assume u solves the

Navier-Stokes equations:

$$u _ { t } + u \cdot \nabla u + \nabla p - \nu \Delta u & = f, \\ \nabla \cdot u & = 0, \\ u ( \cdot, 0 ) & = u _ { 0 }, \\ u | _ { \partial \Omega } & = 0.$$

We define the following bilinear and trilinear forms, respectively:

$$a ^ { N } ( \phi, v ) = ( \nu \nabla \phi, \nabla v ), \quad \forall \phi, v \in X,$$

$$b ( \phi, w, v ) = ( \phi \cdot \nabla w, v ), \quad \forall \phi, w, v \in X.$$

By construction, b ϕ, v, v ( ) = 0 for all ϕ, v ∈ X with ∥∇ · ϕ ∥ = 0, and additionally have the following upper bounds [34]: for any ϕ, w, v ∈ X it holds that for some C &gt; 0 dependent only on the size of Ω,

$$b ( \phi, w, v ) \leq C \| \nabla \phi \| \| \nabla w \| \| \nabla v \|,$$

$$b ( \phi, w, v ) \leq C \left ( \| w \| _ { L ^ { \infty } } + \| \nabla w \| _ { L ^ { 3 } } \right ) \| \phi \| \| \nabla v \|.$$

Algorithm 5.1. Let ∆ t &gt; 0 be the time step size. Given initials condition w 0 h ∈ L 2 (Ω) , w 1 h ∈ X h , source f ∈ L ∞ (0 , ∞ ; L 2 (Ω)) , I H satisfying (5)-(6), observations { I H u n +1 } n ≥ 1 ∈ L 2 (Ω) , and nudging parameter µ &gt; 0 , find { w n h } ∈ X h for n = 1 2 3 , , , · · · , satisfying

$$\left ( \frac { 3 w _ { h } ^ { n + 1 } - 4 w _ { h } ^ { n } + w _ { h } ^ { n - 1 } } { 2 \Delta t }, v _ { h } \right ) + a ^ { N } ( w _ { h } ^ { n + 1 }, v _ { h } ) + b ( 2 w _ { h } ^ { n } - w _ { h } ^ { n - 1 }, w _ { h } ^ { n + 1 }, v _ { h } ) \\ + \mu \left ( I _ { H } ( w _ { h } ^ { n + 1 } - u ^ { n + 1 } ), I _ { H } v _ { h } \right ) = ( f ^ { n + 1 }, v _ { h } ), \quad \forall v _ { h } \in X _ { h }.$$

We now prove a long time optimal L 2 accuracy result for Algorithm 5.1. This improves on the analysis in [41, 32, 20] by providing optimal accuracy in space and time that does not scale with µ when µ gets large.

Theorem 5.1. Let u ∈ L ∞ (0 , T ; H k +1 (Ω)) solve the Navier-Stokes equations (35), u tt ∈ L ∞ (0 , T ; H 1 (Ω)) , u ttt ∈ L ∞ (0 , T ; L 2 (Ω)) , u 0 = u 0 ∈ L 2 (Ω) , u 1 ∈ X , and f ∈ L ∞ (0 , T ; L 2 (Ω)) . Assume the time step ∆ t satisfies

$$\frac { 1 } { 2 } - C \nu ^ { - 1 } \Delta t \left ( h ^ { 2 k + 2 } + h ^ { 2 k } + \| u ^ { n + 1 } \| _ { L ^ { \infty } } ^ { 2 } + \| \nabla u ^ { n + 1 } \| _ { L ^ { 3 } } ^ { 2 } \right ) > 0 \quad a n d \quad \mu \geq \frac { \nu } { 4 C _ { I } ^ { 2 } H ^ { 2 } }.$$

If

$$\frac { \nu } { 4 C _ { I } ^ { 2 } H ^ { 2 } } - C \nu ^ { - 1 } \Delta t \left ( h ^ { 2 k + 2 } + h ^ { 2 k } + \| u ^ { n + 1 } \| _ { L ^ { \infty } } ^ { 2 } + \| \nabla u ^ { n + 1 } \| _ { L ^ { 3 } } ^ { 2 } \right ) > 0,$$

then the difference between the CDA solution of (40) and the true Navier-Stokes equations' solution satisfies: for any n ≥ 1

$$\| w _ { h } ^ { n + 1 } - u ^ { n + 1 } \| ^ { 2 } & \leq C \left ( \frac { 1 } { 1 + \lambda \Delta t } \right ) ^ { n + 1 } \left ( \| [ w _ { h } ^ { 1 } - u ^ { 1 } ; w _ { h } ^ { 0 } - u ^ { 0 } ] \| _ { G } ^ { 2 } + \frac { \kappa \Delta t } { 4 } \| \nabla ( w _ { h } ^ { 1 } - u ^ { 1 } ) \| ^ { 2 } \right ) \\ & \quad + \frac { C \nu ^ { - 1 } \left ( h ^ { 2 k + 2 } + h ^ { 4 k } + \Delta t ^ { 4 } \right ) } { \lambda }. \\ & \quad \cdot \quad \cdot \quad \cdot$$

Here, λ = min { ∆ t -1 , νC 2 l 4 C 2 P , ( ν 4 C H 2 I 2 -Cν -1 ( h 2 k +2 + h 2 k + ∥ u n +1 ∥ 2 L ∞ + ∥∇ u n +1 ∥ 2 L 3 ) ) C 2 l } , and note the constants are independent of µ .

Remark 5.2. Note that for ∆ t and H sufficiently small, λ = νC 2 l 4 C 2 P . Then for n large enough, the L 2 accuracy is O ( h k +1 +∆ ) t 2 .

Proof. Using Taylor's theorem, it can be shown that the true solution of the Navier-Stokes equations (35) satisfies

$$& \left ( \frac { 3 u ^ { n + 1 } - 4 u ^ { n } + u ^ { n - 1 } } { 2 \Delta t }, v \right ) + a ^ { N } ( u ^ { n + 1 }, v ) + b ( 2 u ^ { n } - u ^ { n - 1 }, u ^ { n + 1 }, v ) \\ & = ( f ^ { n + 1 }, v ) + \frac { \Delta t ^ { 2 } } { 3 } ( u _ { t t t } ( t ^ { * } ), v ) + \Delta t ^ { 2 } b ( u _ { t t } ( t ^ { * * } ), u ^ { n + 1 }, v ), \ \forall v \in X, \\ & \, \cdot \quad \,. \quad \,. \quad \,. \quad \,.$$

where t ∗ , t ∗∗ ∈ [ t n -1 , t n +1 ]. Subtracting equation (40) from (43) gives an error equation:

$$\Lambda _ { 1 } ^ { ( 1 + 1 ) } & \Lambda _ { 2 } ^ { ( 1 + 1 ) } \Lambda _ { 3 } ^ { ( 1 + 1 ) } \Lambda _ { 4 } ^ { ( 1 + 1 ) } \Lambda _ { 5 } ^ { ( 1 + 1 ) } \Lambda _ { 6 } ^ { ( 1 + 1 ) } \Lambda _ { 7 } ^ { ( 1 + 1 ) } \Lambda _ { 8 } ^ { ( 1 + 1 ) } \Lambda _ { 9 } ^ { ( 1 + 1 ) } \\ & \quad \left ( \frac { 3 e ^ { n + 1 } - 4 e ^ { n } + e ^ { n - 1 } } { 2 \Delta t }, v \right ) + a ^ { N } ( e ^ { n + 1 }, v _ { h } ) - b ( 2 e ^ { n } - e ^ { n - 1 }, e ^ { n + 1 }, v _ { h } ) \\ & \quad + b ( 2 u ^ { n } - u ^ { n - 1 }, e ^ { n + 1 }, v _ { h } ) + b ( 2 e ^ { n } - e ^ { n - 1 }, u ^ { n + 1 }, v _ { h } ) + \mu \left ( I _ { H } e ^ { n + 1 }, I _ { H } v _ { h } \right ) \\ & = \frac { \Delta t ^ { 2 } } { 3 } ( u _ { t t t } ( t ^ { * } ), v _ { h } ) + \Delta t ^ { 2 } b ( u _ { t t } ( t ^ { * * } ), u ^ { n + 1 }, v _ { h } ), \quad \forall v _ { h } \in X _ { h }, \\ \dots \Lambda _ { 2 } \Lambda _ { 3 } ^ { n } - \Lambda _ { 4 } ^ { n } \quad \dots ^ { n } \quad \dots ^ { n } \quad \dots ^ { n } \quad \dots \Lambda _ { 2 } \Lambda _ { 3 } ^ { n } \end{pmatrix}$$

where e n = u n -w n h . We then decompose e n as

$$e ^ { n } = \left ( u ^ { n } - P _ { h } ^ { S } u ^ { n } \right ) + \left ( P _ { h } ^ { S } u ^ { n } - w _ { h } ^ { n } \right ) = \eta ^ { n } + \phi ^ { n }$$

and, with the help of the projection equation (17), we rewrite (44) as:

$$& \text{and, with the neip or the projection equation } ( 1 / \text{), we rewrite } ( 4 4 ) \text{ as:} \\ & \left ( \frac { 3 \phi ^ { n + 1 } - 4 \phi ^ { n } + \phi ^ { n - 1 } } { 2 \Delta t }, v \right ) + \left ( \frac { 3 \eta ^ { n + 1 } - 4 \eta ^ { n } + \eta ^ { n - 1 } } { 2 \Delta t }, v \right ) + a ^ { N } ( \phi ^ { n + 1 }, v _ { h } ) \\ & \quad - b ( 2 \phi ^ { n } - \phi ^ { n - 1 }, \phi ^ { n + 1 }, v _ { h } ) - b ( 2 \phi ^ { n } - \phi ^ { n - 1 }, \eta ^ { n + 1 }, v _ { h } ) - b ( 2 \eta ^ { n } - \eta ^ { n - 1 }, \phi ^ { n + 1 }, v _ { h } ) \\ & \quad - b ( 2 \eta ^ { n } - \eta ^ { n - 1 }, \eta ^ { n + 1 }, v _ { h } ) + b ( 2 u ^ { n } - u ^ { n - 1 }, \phi ^ { n + 1 }, v _ { h } ) + b ( 2 u ^ { n } - u ^ { n - 1 }, \eta ^ { n + 1 }, v _ { h } ) \\ & \quad + b ( 2 \phi ^ { n } - \phi ^ { n - 1 }, u ^ { n + 1 }, v _ { h } ) + b ( 2 \eta ^ { n } - \eta ^ { n - 1 }, u ^ { n + 1 }, v _ { h } ) + \mu \left ( I _ { H } \phi ^ { n + 1 }, I _ { H } v _ { h } \right ) \\ & \quad = \frac { \Delta t ^ { 2 } } { 3 } ( u _ { t t t } ( t ^ { * } ), v _ { h } ) + \Delta t ^ { 2 } b ( u _ { t t } ( t ^ { * * } ), u ^ { n + 1 }, v _ { h } ), \ \forall v _ { h } \in X _ { h }. \\ & \quad \cdot \cdot \cdot \cdot$$

Taking v h = ϕ n +1 and using the identity b ( · , v, v ) = 0 , ∀ v ∈ X , we have

$$& \text{Taking } v _ { h } = \phi ^ { n + 1 } \text{ and using the identity } b ( \cdot, v, v ) = 0, \ \forall v \in X, \, \text{we have} \\ & \frac { 1 } { 2 \Delta t } ( \| [ \phi ^ { n + 1 } ; \, \phi ^ { n } ] \| _ { G } ^ { 2 } - \| [ \phi ^ { n } ; \, \phi ^ { n - 1 } ] \| _ { G } ^ { 2 } + \frac { 1 } { 2 } \| \phi ^ { n + 1 } - 2 \phi ^ { n } + \phi ^ { n - 1 } \| ^ { 2 } ) + \nu \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ & \quad + \mu \| I _ { H } \phi ^ { n + 1 } \| ^ { 2 } = - \left ( \frac { 3 \eta ^ { n + 1 } - 4 \eta ^ { n } + \eta ^ { n - 1 } } { 2 \Delta t }, \phi ^ { n + 1 } \right ) + \Delta t ^ { 2 } b ( u _ { t t } ( t ^ { * * } ), \, u ^ { n + 1 }, \, \phi ^ { n + 1 } ) \\ & \quad + b ( 2 \phi ^ { n } - \phi ^ { n - 1 }, \, \eta ^ { n + 1 }, \, \phi ^ { n + 1 } ) + b ( 2 \eta ^ { n } - \eta ^ { n - 1 }, \, \eta ^ { n + 1 }, \, \phi ^ { n + 1 } ) - b ( 2 u ^ { n } - u ^ { n - 1 }, \, \eta ^ { n + 1 }, \, \phi ^ { n + 1 } ) \\ & \quad - b ( 2 \phi ^ { n } - \phi ^ { n - 1 }, \, u ^ { n + 1 }, \, \phi ^ { n + 1 } ) - b ( 2 \eta ^ { n } - \eta ^ { n - 1 }, \, u ^ { n + 1 }, \, \phi ^ { n + 1 } ) + \frac { \Delta t ^ { 2 } } { 3 } ( u _ { t t t } ( t ^ { * } ), \, \phi ^ { n + 1 } ). \\ & \text{Applying standard inequalities such as Young's and the Cauchy-Schwarz inequalities, prop-}. \\ & \cdot \, \dots \, \dots \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdots$$

Applying standard inequalities such as Young's and the Cauchy-Schwarz inequalities, properties (38)-(39), we bound the right hand side term of (46) as follows:

$$& - b ( 2 \phi ^ { n - } - \phi ^ { n - 1 }, u ^ { n + 1 }, \phi ^ { n + 1 } ) + \frac { \Delta t ^ { 2 } } { 3 ( u _ { t u } ( t ^ { * } ), \phi ^ { n + 1 } ) }. \\ & \text{Applying standard inequalities such as Young's and the Cauchy-Schwarz inequalities, prop-} \\ & \text{erties (38)-(39), we bound the right hand side term of (46) as follows: } \begin{matrix} \Delta t ^ { 2 } _ { 1 } \\ \frac { \Delta t ^ { 2 } } { 3 } ( u _ { t u } ( t ^ { * } ), \phi ^ { n + 1 } ) \leq C \nu ^ { - 1 } \Delta t ^ { 4 } \| u _ { t u } ( t ^ { * } ) \| ^ { 2 } + \frac { \nu } { 2 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ \Delta t ^ { 2 } b ( u _ { t u } ( t ^ { * } ), u ^ { n + 1 }, \phi ^ { n + 1 } ) \leq C \nu ^ { - 1 } \Delta t ^ { 4 } \| \nabla u ^ { n + 1 } \| ^ { 2 } \| \nabla u _ { t u } ( t ^ { * } ) \| ^ { 2 } + \frac { \nu } { 2 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ - \left ( \frac { 3 \eta ^ { n + 1 } - 4 \eta ^ { n } + \eta ^ { n - 1 } } { 2 \Delta t } \right ) \leq C \nu ^ { - 1 } \left ( \| \eta \| _ { L ^ { 2 } _ { \infty } ( \infty, L ^ { 2 } ) + 1 } \int _ { n - 1 } ^ { n + 1 } \| \eta \| ^ { 2 } d t \right ) + \frac { \nu } { 2 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ b ( 2 \phi ^ { n - } - \phi ^ { n - 1 }, \eta ^ { n + 1 }, \phi ^ { n + 1 } ) = - b ( \phi ^ { n + 1 } - 2 \phi ^ { n } + \phi ^ { n - 1 }, \eta ^ { n + 1 }, \phi ^ { n + 1 } ) + b ( \phi ^ { n + 1 }, \eta ^ { n + 1 }, \phi ^ { n + 1 } ) \\ & \leq C \nu ^ { - 1 } \left ( \| \eta ^ { n + 1 } \| _ { L ^ { 2 } _ { \infty } + \| \nabla \eta ^ { n + 1 } \| _ { L ^ { 2 } _ { \infty } } \right ) \| \phi ^ { n + 1 } - 2 \phi ^ { n } + \phi ^ { n - 1 } \| ^ { 2 } \\ & \quad + C \nu ^ { - 1 } \left ( \| \eta ^ { n + 1 } \| _ { L ^ { 2 } _ { \infty } + \| \nabla \eta ^ { n + 1 } \| ^ { 2 } \right ) \| \phi ^ { n + 1 } \| ^ { 2 } + \frac { \nu } { 1 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ b ( 2 \eta ^ { n - } - \eta ^ { n - 1 }, \eta ^ { n + 1 }, \phi ^ { n + 1 } ) \leq C \nu ^ { - 1 } \| \nabla ( 2 \eta ^ { n - } - \eta ^ { n - 1 } ) \| ^ { 2 } \| \nabla \eta ^ { n + 1 } \| ^ { 2 } + \frac { \nu } { 2 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ - b ( 2 u ^ { n - } - u ^ { n - 1 }, \eta ^ { n + 1 }, \phi ^ { n + 1 } ) = b ( 2 u ^ { n - } - u ^ { n - 1 }, \phi ^ { n + 1 }, \eta ^ { n + 1 } ) \\ & \leq C \nu ^ { - 1 } \| 2 u ^ { n - } - u ^ { n - 1 } \| ^ { 2 } \| \eta ^ { n + 1 } \| _ { L ^ { 2 } _ { \infty } + \frac { \nu } { 2 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ - b ( 2 \phi ^ { n - } - \phi ^ { n - 1 }, u ^ { n + 1 }, \phi ^ { n + 1 } ) = b ( \phi ^ { n + 1 } - 2 \phi ^ { n } + \phi ^ { n - 1 }, u ^ { n + 1 }, \phi ^ { n + 1 } ) - b ( \phi ^ { n + 1 }, u ^ { n + 1 }, \phi ^ { n + 1 } ) \\ & \leq C \nu ^ { - 1 } \left ( \| u ^ { n + 1 } \| _ { L ^ { 2 } _ { \infty } + \| \nabla u ^ { n + 1 } \| _ { L ^ { 2 } _ { \infty } } \right ) \| \phi ^ { n + 1 } \| ^ { 2 } + \frac { \nu } { 1 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ - b ( 2 \eta ^ { n - } - \eta ^ { n - 1 }, u ^ { n + 1 }, \phi ^ { n + 1 } ) = - b ( 2 \eta ^ { n - } - \eta ^ { n - 1 }, u ^ { n + 1 }, \phi ^ { n + 1 } ) \\ & \leq C \nu ^ { - 1 } \left ( \| u ^ { n + 1 } \| _ { L ^ { 2 } _ { \infty } + \| \nabla u ^ { n + 1 } \| _ { L ^ { 2 } _ { \infty } } \right ) \| \eta ^ { n - } - \eta ^ { n - 1 } \| ^ { 2 } \\ & \quad + \frac { \nu } { 2 } \| \nabla \phi ^ { n + 1 } \| ^ { 2 }. \\ & \quad \text{We use the interpolant inequality (5) to lower bound $\nu \| \nabla \phi ^ { n + 1 } \| ^ { 2 } + \mu \| I _ { H } \phi ^ { n + 1 } \| ^ { 2 }, following
 & \text{the analysis of [20], as follows: } \begin{matrix} \nu \| \nabla \phi ^ { n + 1 } \| ^ { 2 } + \mu \| I _ { H } \phi ^ { n + 1 } \| ^ { 2 } \geq \frac { \nu } { 2 \nu } H ^ { 2 } \| \phi ^ { n + 1 } - I _ { H } \phi ^ { n + 1 } \| ^ { 2 } + \mu \| I _ { H } \phi ^ { n + 1 } \| ^ { 2 } + \frac { 3 \nu } { 4 \nu } \| \nabla \phi ^ { n + 1 } \| ^ { 2 } \\ \end{matrix}.$$

We use the interpolant inequality (5) to lower bound ν ∥∇ ϕ n +1 ∥ 2 + µ I ∥ H ϕ n +1 ∥ 2 , following the analysis of [20], as follows:

$$& \nu \| \nabla \phi ^ { n + 1 } \| ^ { 2$$

where γ = min { µ, ν 4 C H 2 I 2 } = ν 4 C H 2 I 2 by allowing µ to be sufficient large, i.e., µ ≥ ν 4 C H 2 I 2 .

Combining (46), (47), and (48), multiplying by 2∆ , and taking into account the assumpt tions on the regularity of the solution u while utilizing standard approximation properties, leads to

$$\text{ leads to} \\ & \| [ \phi ^ { n + 1 } ; \phi$$

Choose H and ∆ t so that ν 4 C H 2 I 2 -Cν -1 ∆ t ( h 2 k +2 + h 2 k + ∥ u n +1 ∥ 2 L ∞ + ∥∇ u n +1 ∥ 2 L 3 ) &gt; 0 and 1 2 -Cν -1 ∆ t ( h 2 k +2 + h 2 k + ∥ u n +1 ∥ 2 L ∞ + ∥∇ u n +1 ∥ 2 L 3 ) &gt; 0 are satisfied. With these choices, using the Poincar´ inequality, G-norm equivalence (3), CDA Stokes projection error estimate, e and proceeding similarly as (33)-(34), leads to

$$\begin{array} {$$

where λ = min { ∆ t -1 , νC 2 l 4 C 2 P , ( ν 4 C H 2 I 2 -Cν -1 ( h 2 k +2 + h 2 k + ∥ u n +1 ∥ 2 L ∞ + ∥∇ u n +1 ∥ 2 L 3 ) ) C 2 l } . Then, the proof is completed by applying the triangle inequality.

## 6. Numerical Experiments

We now illustrate the effectiveness of the proposed method for several test problems. For all of our tests, we take the L 2 projection onto piecewise constants as I H , implemented with algebraic nudging [41].

## 6.1. The heat equation

For our first test, we pick a problem with domain Ω = (0 1) , , 2 κ = 1, and analytical solution

$$u = \sin ( t + 2 \pi x + \pi y ).$$

The solution then determines f = u t -∆ , and the boundary conditions are taken to be u Dirichlet on all of ∂ Ω and equal to the solution.

We compute with the BDF2-FEM CDA algorithm above using H = 1 9 and µ = ∞ (direct implementation) with P 2 elements on barycenter refinements of uniform triangular meshes. We note that we also tested with µ = 10 5 , 10 8 and saw no difference in results between these solutions and those found using direct implementation.

The theory above from Theorem 4.1 shows that for µ sufficiently large and for n large enough,

$$\| w _ { h } ^ { n } - u ( t ^ { n } ) \|$$

We illustrate this theory by first taking small T = 0 3 and ∆ . t = 0 001, and varying . h . Errors and rates at the final time are shown in Table 1, and we observe third order convergence in space, which is consistent with the theory since we use P 2 elements. To illustrate temporal convergence rates, we take h = 1 256 , T = 1, and vary ∆ . t Errors and rates at the final time are shown in Table 2, and we observe rates consistent with second order temporal convergence (the last error and rate are reduced due to temporal error no longer being completely dominant over spatial error with the smallest choice of ∆ ). t Plots of the error measured in the L 2 -norm for all of these tests are shown in Figure 1, and the results are consistent with the theory - exponential convergence in time up to the discretization error, then long-time optimal accuracy, i.e. O h ( 3 + ∆ t 2 ) since P 2 elements are used for this test, with respect to the L 2 -norm.

Figure 1: Shown above are L 2 error versus time for the tests using the analytical solution for the heat equation and BDF2-CDA algorithm, for varying h (left) and varying ∆ t (right).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00396v1-figure-f04032ea-f896-44db-8fd1-fa1f553d186f.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00396v1-figure-6deb46c0-dda1-4bc2-8482-62872148a7d3.jpg)

Table 1: Shown above are errors at the final time and spatial convergence rates for T = 0 3 and ∆ . t = 0 001 . ( M = T/ ∆ ). t

| h     |   ∥ w M h - u (0 . 3) ∥ | rate   |
|-------|-------------------------|--------|
| 1/32  |               0.000469  | -      |
| 1/64  |               4.947e-05 | 3.25   |
| 1/128 |               5.865e-06 | 3.08   |
| 1/256 |               7.236e-07 | 3.02   |

|   ∆ t |   ∥ w M h - u (1) ∥ | rate   |
|-------|---------------------|--------|
| 0.2   |           6.985e-05 | -      |
| 0.1   |           1.87e-05  | 1.90   |
| 0.05  |           4.712e-06 | 1.99   |
| 0.025 |           1.372e-06 | 1.78   |

Table 2: Shown above are errors and temporal convergence rates for h = 1 256 , T = 1 0 and varying ∆ . t ( M = T/ ∆ ). t

## 6.2. A contaminant transport model

In this numerical experiment, a CDA-FEM of the heat equation with added linear transport that models contaminant transport in a river is investigated. The linear transport equation takes the form

$$u _ { t } + U \cdot \nabla u - \kappa \Delta u \ = \ f,$$

$$u ( 0 ) \ = \ u _ { 0 },$$

where U is a given divergence free velocity field, κ is the (material dependent) diffusion coefficient, and here u represents the concentration of a contaminant. This system is similar to heat equation with an additional linear term U · ∇ u , and thus CDA can be applied in the same way as above for the heat equation with (essentially) the same convergence results.

We follow a test problem from [28], the domain for which is the area enclosed by the curves y = sin( x ), y = 1 + sin( x ), x = 0 and x = 4 π as left and right boundaries. Using the mesh τ h shown in Figure 2, we use ( P , P 2 1 ) Taylor-Hood velocity-pressure elements (which gives 13,928 degrees of freedom (dof)) to compute the Stokes equations with viscosity 0 01 . and f = 0, using no-slip boundary conditions on the top and bottom boundaries, a plug inflow of u in = 3, and a zero-traction outflow enforced with the do-nothing condition. We take our transport velocity U to be the velocity solution of this discrete Stokes problem.

Figure 2: Shown above are the fine mesh used for the finite element computations (left) and the coarse mesh nodes used for CDA (right).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00396v1-figure-2217dadb-b0ce-4088-a75d-986e3a1c08f0.jpg)

We use BDF2-FEM with CDA to solve (51)-(52) equipped with a boundary condition of 0 at the inflow ( c in = 0), the transport velocity U coming from the discrete Stokes problem, transport term discretization

$$b ^ { * } ( U, w _ { h } ^ { n + 1 }, v ) \coloneqq b ( U, w _ { h } ^ { n + 1 }, v _ { h } ) + \frac { 1 } { 2 } ( ( \nabla \cdot U ) w _ { h } ^ { n + 1 }, v _ { h } )$$

Figure 3: Shown above are plots of convergence of the BDF2-CDA solutions to the true solution, for varying µ .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00396v1-figure-5d35b5fa-3dcc-4afc-bced-df485cb8d45b.jpg)

where b U, w ( n +1 h , v h ) := ( U · ∇ w n +1 h , v h ), κ = 0 01, and zero Neumann conditions . ∇ · c n = 0 at all other boundaries. For the initial condition, we take zero contaminant except for two 'blobs', which are represented with c = 3 inside the circles with radius 0 1 centered at (1 . , 1 5) . and (5 , -0 5). . We first compute a direct numerical simulation (DNS) for the concentration c using BDF2-FEM, P 2 on the same mesh as is used for the Stokes problem, and time step ∆ = 0 02. t .

Finally, we compute the CDA solution, using BDF2-CDA with the same parameters as the DNS except zero initial conditions, taking the DNS solution u as the true solution, and testing the algorithm using several choices of µ including µ = ∞ which we implement with direct enforcement. The 37 measurement nodes are shown in Figure 2 and we note that these nodes are also vertices on the finite element mesh. Convergence of the CDA solution to the true (DNS) solution is shown in Figure 3, for varying µ . We observe very similar convergence for each choice of µ =0.1, 1, 10000, and ∞ , with the plots for µ = 10000 and ∞ (direct enforcement) laying on top of each other. For the case µ = ∞ , we show contour plots of the DA and true (DNS) solution at t=0, 1, 2, and 5 in Figure 4, and observe that the agreement between the solutions increases as time-stepping moves forward and that the solutions are visually indistinguishable at the final stopping time.

## 6.3. The Navier-Stokes equations

To illustrate the theory for the Navier-Stokes equations, we choose a test problem from [42, 39] for simulating 2D Kelvin-Helmholtz instability. The domain is taken to be the unit square and enforce periodic boundary conditions on x = 0 1 (representing an infinite , extension in the horizontal direction). At y = 0 1, we enforce a no penetration and free slip , condition. The initial condition for the benchmark problem is set as

$$u _ { 0 } ( x, y ) = \left ( \begin{array} { c } u _ { \infty } \tanh \left ( \frac { 2 y - 1 } { \delta _ { 0 } } \right ) \\ 0 \end{array} \right ) + c _ { n } \left ( \begin{array} { c } \partial _ { y } \psi ( x, y ) \\ - \partial _ { x } \psi ( x, y ) \end{array} \right ),$$

Figure 4: Contour plots of DA and DNS velocity magnitudes at times 0, 1, 2 and 5.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00396v1-figure-02a4aea1-cb46-4fa2-a419-c41e36e35a28.jpg)

Figure 5: Shown above are convergence plots of the CDA solution to the true solution for the KelvinHelmholtz test problem.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00396v1-figure-54ef54ef-518f-4801-8ffa-7a849c833b32.jpg)

where δ 0 = 1 28 is the initial vorticity thickness, u ∞ = 1 is set as the reference velocity, c n is a noise/scaling factor that is taken to be 10 -3 , and

$$\psi ( x, y ) = u _ { \infty } \exp \left ( - \frac { ( y - 0. 5 ) ^ { 2 } } { \delta _ { 0 } ^ { 2 } } \right ) \left ( \cos ( 8 \pi x ) + \cos ( 2 0 \pi x ) \right ).$$

The Reynolds number for this test is defined by Re = δ 0 u ∞ ν = 1 28 ν , and ν is defined by selecting Re . We use Re =100 for this test.

We compute solutions for the BDF2-CDA for the NSE similar to the previous experiments: first compute a reference solution (using the initial condition above), and then use this reference solution as the true solution to nudge towards. We use ( P , P 2 disc 1 ) Scott-Vogelius elements on a barycenter refinement of an h = 1 96 uniform triangular mesh, and a time / step size of ∆ = 0 01. t . For the CDA, we take the initial condition to be w 0 h = w 1 h = 0, and compute using Algorithm 5.1 with µ = ∞ (implemented with direct enforcement) and varying H (=1/25, 1/20, and 1/15) using N = 1 H 2 measurement points. Convergence of the CDA solutions to the true solutions are shown in Figure 5, and we observe that the method is converging in L 2 for all three choices of H with faster convergence as H decreases. Contour plots of the solution's absolute vorticity magnitudes are shown in Figure 6, along with those of the true solution. At t = 2, none of the CDA solutions appear accurate. By t = 5, the N=20 and 25 solutions are already approaching the true solution, and by t = 8 their solutions are visually indistinguishable from the true solution. At t = 9, there are still visible errors in the N = 15 solution, demonstrating that the inclusion of more information yields faster convergence.

## 6.4. The Cahn-Hilliard equation

In this subsection, we briefly demonstrate via a numerical experiment that similar CDA projections could be established for different types of evolutionary PDEs. Specifically, we consider the Cahn-Hilliard (CH) equation [10, 11, 38]

$$\partial _ { t } \varphi - \Delta \left ( \varphi ^ { 3 } - \varphi \right ) + \varepsilon ^ { 2 } \Delta ^ { 2 } \varphi & = 0, \\ \varphi ( \cdot, 0 ) & = \varphi _ { 0 },$$

where φ represents the order parameter which takes on values between -1 and 1 and is often interpreted as a concentration of one component in a two component system. The states φ = ± 1 indicate phases of pure concentration and ε &gt; 0 is understood to represent an interfacial width between the two phases. The following CDA-FEM algorithm for the CH equation (53) was analyzed in [16]: Given ϕ n -1 h ∈ Z h and true solution φ ∈ L ∞ (0 , T ; H 2 N (Ω)), find ϕ n h ∈ Z h such that

$$\left ( \frac { \phi _ { h } ^ { n } - \phi _ { h } ^ { n - 1 } } { \Delta t }, \psi \right ) + \left ( \nabla \left ( ( \phi _ { h } ^ { n } ) ^ { 3 } - \phi _ { h } ^ { n - 1 } \right ), \nabla \psi \right ) + \varepsilon ^ { 2 } a _ { h } ^ { I P } \left ( \phi _ { h } ^ { n }, \psi \right ) + \mu \left ( I _ { H } \left ( \phi _ { h } ^ { n } - \varphi ^ { n } \right ), \psi \right ) = 0,$$

Figure 6: Contour plots of DA and DNS absolute vorticity magnitudes at times 2, 5, 8, and 9, reading from top to bottom.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00396v1-figure-614a1dce-b5f1-48c3-9276-957105b51d39.jpg)

for all ψ ∈ Z h with Z h := H 1 (Ω) ∩ P k ( τ h (Ω)). Here, the bilinear form a IP h ( · , · ) has the following definition:

$$\lambda _ { h } ^ { I P } \left ( w, v \right ) \coloneqq & \sum _ { K \in \tau _ { h } } \int _ { K } \left ( \nabla ^ { 2 } w \, \colon \nabla ^ { 2 } v \right ) \, d x + \sum _ { e \in \mathcal { E } _ { h } } \int _ { e } \left \{ \left \{ \frac { \partial ^ { 2 } w } { \partial n _ { e } ^ { 2 } } \right \} \right \} \left [ \frac { \partial v } { \partial n _ { e } } \right ] \right ] d S \\ & + \sum _ { e \in \mathcal { E } _ { h } } \int _ { e } \left \{ \left \{ \frac { \partial ^ { 2 } v } { \partial n _ { e } ^ { 2 } } \right \} \right \} \left [ \frac { \partial w } { \partial n _ { e } } \right ] d S + \sigma \sum _ { e \in \mathcal { E } _ { h } } \frac { 1 } { | e | } \int _ { e } \left [ \left \{ \frac { \partial w } { \partial n _ { e } } \right ] \right ] \left [ \frac { \partial v } { \partial n _ { e } } \right ] d S,$$

with σ ≥ 1 known as a penalty parameter. The jumps and averages that appear in (55) are defined as follows. For an interior edge e shared by two triangles K ± where n e points from K -to K + , we define on the edge e

$$\left [ \frac { \partial v } { \partial n _ { e } } \right ] = n _ { e } \cdot ( \nabla v _ { + } - \nabla v _ { - } ) \,, \quad \left \{ \left \{ \frac { \partial ^ { 2 } v } { \partial n _ { e } ^ { 2 } } \right \} \right \} = \frac { 1 } { 2 } \left ( \frac { \partial ^ { 2 } v _ { - } } { \partial n _ { e } ^ { 2 } } + \frac { \partial ^ { 2 } v _ { + } } { \partial n _ { e } ^ { 2 } } \right ),$$

where ∂ v 2 ∂n 2 e = n e · ( ∇ 2 v ) n e and v ± = v | K ± . For a boundary edge e , we take n e to be the unit normal pointing towards the outside of Ω and define on the edge e

$$\left [ \frac { \partial v } { \partial n _ { e } } \right ] = - n _ { e } \cdot \nabla v _ { K }, \quad \left \{ \left \{ \frac { \partial ^ { 2 } v } { \partial n _ { e } ^ { 2 } } \right \} \right \} = n _ { e } \cdot \left ( \nabla ^ { 2 } v \right ) n _ { e }.$$

In [16], it was once again determined that an upper bound to the nudging parameter µ was required in the analysis of the FEM CDA scheme. However, the following numerical experiment suggests otherwise, in which we consider a square domain Ω = (0 1) , 2 and take τ h (Ω) to be a regular triangulation of Ω consisting of right isosceles triangles that are a quasi-uniform family. (We use a family of meshes τ h (Ω) such that no triangle in the mesh has more than one edge on the boundary.)

As the Firedrake Project [24] was utilized this experiment, we now describe how the data assimilation term µ I ( H ( ϕ n h -φ n ) , I H ψ ) was implemented. Specifically, a true solution φ was obtained at all times by selecting a cross shaped region as initial conditions, setting the interfacial width parameter ε = 0 002 and the nudging parameter . µ = 0, and solving the CH equation using the C 0 interior penalty FEM (54). A data assimilation grid size H was chosen and grid points were identified and located on the finite element mesh. Finally, the initial conditions for the numerical solution ϕ 0 was set to random initial conditions.

Figure 7 demonstrates the effectiveness of the CDA-FEM for various grid sizes H with a large nudging parameter of µ = 10 . 9 Specifically, we choose five different grid sizes H = 0.011049, 0.015625, 0.03125, 0.0625, and 0.125, which correspond respectively to 8,100, 4,096, 1,024, 256 and 64 grid points, while the fine mesh uses piecewise quadratics and has 33,025 grid points. Theorem 4.3 from [16] provides a sufficient condition that the grid size should be chosen as H = O ( ε 2 ) but our experiments suggest that a grid size much coarser than that will produce good results if a large enough nudging parameter is chosen. This also suggests that the new theory developed herein should allow for a relaxation of the restriction on the size of the data assimilation grid.

Figure 8 demonstrates the effectiveness of the CDA-FEM for various values of the nudging parameter µ . For this experiment, we set the data assimilation grid to be H = 0 03125 . , and choose four different values for the nudging parameter µ = 2500 5000 10000 , , , and 10 . 9 Theorem 4.3 from [16] admits a sufficient condition that the appropriate value for the nudging parameter µ is at least 1 / ε 2 = 2500, but if µ is too large then H needs to be very small. However, our experiments show that good results can also be obtained for much larger values of µ without requiring H to be very small.

Figure 7: The L 2 error between the true solution and the solution to the data assimilation finite element method for various data assimilation grid sizes. The mesh size is h = √ 2 64 / and the time step size is ∆ = t 0 002 . / 32 . All other parameters are defined in the text.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00396v1-figure-64ff5cc2-0af7-4ff3-972c-ae743375ac73.jpg)

Figure 8: The L 2 error between the true solution and the solution to the CDA-FEM for various values of the nudging parameter. The mesh size is h = √ 2 64 / and the time step size is ∆ t = 0 002 . / 32 . All other parameters are defined in the text.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00396v1-figure-7db43ea6-9703-4f65-935c-0f77571ef479.jpg)

## 7. Conclusion

In this paper, we proved optimal long time convergence in L 2 for finite element methods with continuous data assimilation without requiring any upper restrictions on the nudging parameter µ, and without any negative influence of µ on the error estimate. The goal was achieved by proposing a CDA Poisson projection along with an optimal projection error estimate, in which a harmonic approximation and nodal interpolant played key roles in the the analysis. We then applied the CDA Poisson projection to CDA convergence analysis for two typical evolution PDEs: the heat equation and the Navier-Stokes equations. The key contribution here is the fact that with these new projection, we were able to eliminate

the scaling of the error with µ as µ gets large. This in turn allows for choosing large or very large µ , with the latter being numerically equivalent to direct Dirichlet enforcement of the measurement data into the discrete solution.Multiple numerical tests illustrate the theory.

Future work will consider the case of how to choose µ for noisy data. In many applications of CDA, the data measurements will have error present, and it is an open question of how µ should be chosen if certain properties of the noise are known, such as signal to noise ratio or distribution of the noise.

## 8. Acknowledgements

Authors Xuejian Li and Leo G. Rebholz were partially supported by NSF grant No. DMS2152623. Author Amanda E. Diegel was partially supported by NSF grant No. DMS-2110768. Also, the authors thank Dr. Quyuan Lin for helpful discussions regarding to this work.

## References

- [1] D. Albanez, H. Nussenzveig Lopes, and E. Titi. Continuous data assimilation for the three-dimensional Navier-Stokesα model. Asymptotic Anal. , 97(1-2):139-164, 2016.
- [2] R.A. Anthes. Data assimilation and initialization of hurricane prediction models. J. Atmos. Sci. , 31(3):702-719, 1974.
- [3] M.S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp. A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE T. Signal Proces. , 50(2):174-188, 2002.
- [4] D. Auroux and M. Nodet. The back and forth nudging algorithm for data assimilation problems: theoretical results on transport equations. ESAIM Control Optim. Calc. Var. , 18(2):318-342, 2012.
- [5] D. Aurouxa and J. Blum. Back and forth nudging algorithm for data assimilation problems. C. R. Math. Acad. Sci. Paris , 340(12):873-878, 2005.
- [6] S. Axler, P. Bourdon, and W. Ramey. Harmonic Function Theory . Springer, New York, 2001.
- [7] A. Azouani, E. Olson, and E.S. Titi. Continuous data assimilation using general interpolant observables. J. Nonlinear Sci. , 24:277-304, 2014.
- [8] H. Bessaih, E. Olson, and E.S. Titi. Continuous data assimilation with stochastically noisy data. Nonlinearity , 28(3):729-753, 2015.
- [9] A. Biswas, J. Hudson, A. Larios, and Y. Pei. Continuous data assimilation for the 2D magnetohydrodynamic equations using one component of the velocity and magnetic fields. Asymptot. Anal. , 108(1-2):1-43, 2018.

| [10]   | J.W. Cahn. On spinodal decomposition. Acta Metall Mater , 9(9):795-801, 1961.                                                                                                                                                                                                |
|--------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [11]   | J.W. Cahn and J.E. Hilliard. Free energy of a nonuniform system. I. interfacial free energy. J. Chem. Phys. , 28(2):258-267, 1958.                                                                                                                                           |
| [12]   | Y. Cao, A. Giorgini, M. Jolly, and A. Pakzad. Continuous data assimilation for the 3D Ladyzhenskaya model: analysis and computations. Nonlinear Anal. Real World Appl. 68(103659):1-29, 2022.                                                                                |
| [13]   | E. Carlson, J. Hudson, and A. Larios. Parameter recovery for the 2 dimensional Navier- Stokes equations via continuous data assimilation. SIAM J.Sci. Comput. , 42(1):A250- A270, 2020.                                                                                      |
| [14]   | P. Courtier, E. Andersson, W. Heckley, J. Pailleux, D. Vasiljevic, M. Hamrud, A. Hollingsworth, F. Rabier, and M. Fisher. The ECMWF implementation of three- dimensional variational assimilation (3D-Var). I:formulation. Q. J. R. Meteorol. Soc. 124(550):1783-1807, 1998. |
| [15]   | P. Courtier, J.N. Thepaut, and A Hollingsworth. A strategy for operational im- plementation of 4D-var, using an incremental approach. Q. J. Roy. Meteor. Soc. 120(519):1367-1387, 1994.                                                                                      |
| [16]   | A.E. Diegel and L.G. Rebholz. Continuous data assimilation and long-time accuracy in a C 0 interior penalty method for the Cahn-Hilliard equation. Appl. Math. Comput. 424:127042, 2022.                                                                                     |
| [17]   | G. Evensen. The ensemble Kalman filter for combined state and parameter estimation: Monte Carlo techniques for data assimilation in large systems. IEEE Control Syst. Mag. , 29(3):83-104, 2009.                                                                             |
| [18]   | A. Farhat, M.S. Jolly, and E.S. Titi. Continuous data assimilation for the 2d B´nard e convection through velocity measurements alone. Physica D. , 303:59-66, 2015.                                                                                                         |
| [19]   | P. Fearnhead and H.R. K¨nsch. u Particle filters and data assimilation. Annu. Rev. Stat. Appl. , 5(1):421-449, 2018.                                                                                                                                                         |
| [20]   | B. Garcia-Archilla and J. Novo. Error analysis of fully discrete mixed finite element data assimilation schemes for the Navier-Stokes equations. Adv. Comput. Math. , 46(61): 1-33, 2020.                                                                                    |
| [21]   | B. Garcia-Archilla, J. Novo, and E.S. Titi. Uniform in time error estimates for a finite element method applied to a downscaling data assimilation algorithm. SIAM J. Numer. Anal. , 58:410-429, 2020.                                                                       |
| [22]   | M. Gardner, A. Larios, L.G. Rebholz, D. Vargun, and C. Zerfas. Continuous data as- similation applied to a velocity-vorticity formulation of the 2D Navier-Stokes equations. Electronic Research Archive , 29(3):2223-2247, 2021.                                            |

| [23]   | E. Hairer and G. Wanner. Solving ordinary differential equations. II , volume 14 of Springer Series in Computational Mathematics . Springer-Verlag, Berlin, second edition, 1996. Stiff and differential-algebraic problems.                                                                                                                                                                                                                                                                                                      |
|--------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [24]   | D.A. Ham, P.H.J. Kelly, L. Mitchell, C.J. Cotter, R.C. Kirby, K. Sagiyama, N. Bouziani, S. Vorderwuelbecke, T.J. Gregory, J. Betteridge, D.R. Shapero, R.W. Nixon-Hill, C.J. Ward, P.E. Farrell, P.D. Brubeck, I. Marsden, T.H. Gibson, M. Homolya, T. Sun, A.T.T. McRae, F. Luporini, A. Gregory, M. Lange, S.W. Funke, F. Rathgeber, G-T. Bercea, and G.R. Markall. Firedrake User Manual . Imperial College London and University of Oxford and Baylor University and University of Washington, first edition edition, 5 2023. |
| [25]   | E. Hawkins, L. Rebholz, and D. Vargun. Removing splitting/modeling error in projec- tion/penalty methods for Navier-Stokes simulations with continuous data assimilation. Communications in Mathematical Research , 40(1):1-29, 2024.                                                                                                                                                                                                                                                                                             |
| [26]   | T. Heister, L. G. Rebholz, and M.Y. Xiao. Flux-preserving enforcement of inhomoge- neous Dirichlet boundary conditions for strongly divergence-free mixed finite element methods for flow problems. J. Math. Anal. Appl. , 438(1):507-513, 2016.                                                                                                                                                                                                                                                                                  |
| [27]   | J.E. Hoke and R.A. Anthes. The initialization of numerical models by a dynamic- initialization technique. Monthly Weather Review , 104(12):1551-1556, 1976.                                                                                                                                                                                                                                                                                                                                                                       |
| [28]   | J. Honeycutt, H. Johnson, and S. Kelly. Using data assimilation to better predict contaminant transport in fluids. SIAM Undergraduate Research Online , 11:1-16, 2018.                                                                                                                                                                                                                                                                                                                                                            |
| [29]   | A.H. Ibdah, C.F. Mondaini, and E.S. Titi. Fully discrete numerical schemes of a data assimilation algorithm: uniform-in-time error estimates. IMA J. Numer. Anal. 40(4):2584-2625, 11 2019.                                                                                                                                                                                                                                                                                                                                       |
| [30]   | M.S. Jolly, V.R. Martinez, E.J. Olson, and E.S. Titi. Continuous data assimilation with blurred-in-time measurements of the surface quasi-geostrophic equation. Chin. Ann. Math. Ser. B , 40(5):721-764, 2019.                                                                                                                                                                                                                                                                                                                    |
| [31]   | R.E. Kalman. A new approach to linear filtering and prediction problems. J. Basic Eng. , 82(1):35-45, 1960.                                                                                                                                                                                                                                                                                                                                                                                                                       |
| [32]   | A. Larios, L.G. Rebholz, and C. Zerfas. Global in time stability and accuracy of IMEX- FEM data assimilation schemes for Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg. , 345:1077-1093, 2019.                                                                                                                                                                                                                                                                                                                       |
| [33]   | A. Larios and C. Victor. Continuous data assimilation with a moving cluster of data points for a reaction diffusion equation: Acomputational study. Commun. Comp. Phys. 29(4):1273-1298, 2021.                                                                                                                                                                                                                                                                                                                                    |
| [34]   | W. Layton. An Introduction to the Numerical Analysis of Viscous Incompressible Flows SIAM, Philadelphia, 2008.                                                                                                                                                                                                                                                                                                                                                                                                                    |

| [35]   | J. Lewis and S. Lakshmivarahan. Sasaki´ s pivotal contribution: calculus of variations applied to weather map analysis. Mon. Weather Rev. , 136(9):3553-3567, 2008.                                                                      |
|--------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [36]   | X. Li, E.V. Hawkins, L.G. Rebholz, and D. Vargun. Accelerating and enabling conver- gence of nonlinear solvers for Navier-Stokes equations by continuous data assimilation. Comput. Methods Appl. Mech. Engrg. , 416:1-17, 2023.         |
| [37]   | X. Li, G. Wei, X.M. He, and T. Lin. Variational data assimilation and its decoupled iterative numerical algorithms for Stokes-Darcy model. SIAM J. Sci. Comput. , 2023.                                                                  |
| [38]   | A. Miranville. The Cahn-Hilliard Equation: Recent Advances and Applications . SIAM, 2019.                                                                                                                                                |
| [39]   | M. Olshanskii and L. Rebholz. Local conservation laws of continuous Galerkin method for the incompressible Navier-Stokes equations in EMAC form. Computer Methods in Applied Mechanics and Engineering , 418B(116583):1-14, 2024.        |
| [40]   | F. Rawlins, S.P. Ballard, K.J. Bovis, A.M. Clayton, D. Li, G.W. Inverarity, A.C. Lorenc, and T.J. Payne. The met office global four-dimensional variational data assimilation scheme. Q. J. R. Meteorol. Soc. , 133(623):347-362, 2007.  |
| [41]   | L.G. Rebholz and C. Zerfas. Simple and efficient continuous data assimilation of evolu- tion equations via algebraic nudging. Numer. Methods Partial Differ. Equ. , 37(3):2588- 2612, 2021.                                              |
| [42]   | P. Schroeder, V. John, P. Lederer, C. Lehrenfeld, G. Lube, and J. Schoberl. On reference solutions and the sensitivity of the 2D Kelvin-Helmholtz instability problem. Computers and Mathematics with Applications , 77:1010-1028, 2019. |
| [43]   | C. Zerfas, L. Rebholz, M. Schneier, and T. Iliescu. Continuous data assimilation reduced order models of fluid flow. Comput. Methods Appl. Mech. Eng. , 357(112596):1-21, 2019.                                                          |
| [44]   | X. Zou, I.M. Navon, and F.X. Ledimet. An optimal nudging data assimilation scheme using parameter estimation. Q. J. R. Meteorol. Soc , 118(508):1163-1186, 1992.                                                                         |
Amanda Diegel Xuejian Li Leo G. Rebholz
Analysis of continuous data assimilation with large (or even infinite) nudging parameters
This paper considers continuous data assimilation (CDA) in partial
differential equation (PDE) discretizations where nudging parameters can be
taken arbitrarily large. We prove that long-time optimally accurate solutions
are obtained for such parameters for the heat and Navier-Stokes equations
(using implicit time stepping methods), with error bounds that do not grow as
the nudging parameter gets large. Existing theoretical results either prove
optimal accuracy but with the error scaled by the nudging parameter, or
suboptimal accuracy that is independent of it. The key idea to the improved
analysis is to decompose the error based on a weighted inner product that
incorporates the (symmetric by construction) nudging term, and prove that the
projection error from this weighted inner product is optimal and independent of
the nudging parameter. We apply the idea to BDF2 - finite element
discretizations of the heat equation and Navier-Stokes equations to show that
with CDA, they will admit optimal long-time accurate solutions independent of
the nudging parameter, for nudging parameters large enough. Several numerical
tests are given for the heat equation, fluid transport equation, Navier-Stokes,
and Cahn-Hilliard that illustrate the theory.