Approximate continuous data assimilation of the 2D Navier-Stokes equations via the Voigt-regularization with observable data

We propose a data assimilation algorithm for the 2D Navier-Stokes equations, based on the Azouani, Olson, and Titi (AOT) algorithm, but applied to the 2D Navier-Stokes-Voigt equations. Adapting the AOT algorithm to regularized versions of Navier-Stokes has been done before, but the innovation of this work is to drive the assimilation equation with observational data, rather than data from a regularized system. We first prove that this new system is globally well-posed. Moreover, we prove that for any admissible initial data, the $L^2$ and $H^1$ norms of error are bounded by a constant times a power of the Voigt-regularization parameter $\alpha>0$, plus a term which decays exponentially fast in time. In particular, the large-time error goes to zero algebraically as $\alpha$ goes to zero. Assuming more smoothness on the initial data and forcing, we also prove similar results for the $H^2$ norm.


Introduction
In real weather simulations, it is common practice to use some type of "filtered" or "regularized" model rather than the original equations to stabilize under-resolved simulations. For instance, one might use a large-eddy simulation (LES) model or an α-model, rather than the pure Navier-Stokes equations (see, e.g., [7,81,97] and the references therein). Somewhat separate from this approach is a class of techniques known as data assimilation, which aim to improve the accuracy of simulations by incorporating observational data into the dynamical system being solved. Combining these approaches is a natural idea, but a problem arises immediately because data from observations obviously does not originate from the regularized model, but from the underlying physical model 1 . This mismatch, and the error that results from it, is the subject of the present work. In particular, we show that, in the setting we consider, the error decays exponentially fast in time, but only down to a level which is determined by the regularizing parameter α. However, we also show that this lower level (which results from the mismatch between models) goes to zero as a power of α when α → 0.
We focus on the method of continuous data assimilation proposed by Azouani, Olson, and Titi in [4,5] (see also [18,54,90] for early ideas in this direction), which we refer to as the AOT algorithm. For simplicity, we consider the 2D Navier-Stokes (NS) equations to be the underlying physical model, and take the "filtered" model to be the 2D Navier-Stokes-Voigt (NS-Voigt) equations. While the 2D NS equations are already globally well-posed in the sense of possessing unique smooth solutions for all time, numerical simulations of these equations could still be under-resolved, and therefore, it is reasonable to consider filtered or smoothed versions of the equations. We therefore consider an algorithm where data arises from the NS equations, but it is incorporated into the NS-Voigt equations via an AOT-type feedbackcontrol term. We show analytically that, so long as the parameters satisfy certain bounds, the error in the solutions decays exponentially fast in time until it reaches a certain level, and that moreover, this error level goes to zero as the Voigt-regularization parameter α goes to zero.
We note that the phenomenon of exponential decay of the error down to a level which can be controlled has recently been observed in a different context. Namely, in [20], the effect of using mismatched viscosity terms is examined. One viscosity (equivalently, Reynolds number) is used for the "physical" equation, and a different viscosity (or Reynolds number) is used for the "assimilation" equation. Such a situation might arise, for example, when the true viscosity of the system is unknown, or when simulations use a smaller Reynolds number than the physical Reynolds number to stabilize the simulation. It is proven in [20] that the error decays exponentially down to a level that is controlled by the difference in the two viscosities, and that this level decreases to zeros as one viscosity approaches the other. This result parallels the results of this paper, albeit in a very different context.
The NS-Voigt equations, and their inviscid counterpart, called the Euler-Voigt equations have been studied analytically and extended in a wide variety of contexts (see, e.g., [6, 8-10, 15, 16, 19, 21, 22, 29, 37, 49, 50, 63, 64, 67, 70, 71, 76, 77, 82, 84, 88, 91-93, 95], and the references therein). Computational studies were carried out in [36,68,74,79]. The NS-Voigt equations were first proposed by Oskolkov in [92,93] as a model for Kelvin-Voigt fluids, but were later viewed as a regularization for the NS equations in [19], where also the Euler-Voigt equations were first introduced and studied. The Voigt-regularization is related to the wider class of α-models, including the NS-α (NS-α) model and the Leray-α model, which were first proposed and studied in [19, 24-27, 45, 56, 58]. The Voigt model enjoys two features that the other α-models are not known to have in the 3D case. First, it is known to be globally well-posed in the inviscid case. Second, in the viscous case, it is well-posed in the physical case of "no-slip" homogeneous Dirichlet boundary conditions, with no need to impose artificial boundary conditions. Although we only work in 2D, we focus on the NS-Voigt model due in part to these attractive features, and also due to its simplicity. The study of the AOT-algorithm applied to other models but still driven by observable data, and the errors resulting from the mismatch, will be the subject of a forthcoming work.
In [1], it was shown that the AOT data assimilation scheme for the 3D NS-α model (with any admissible initial data) has solutions which converge to solutions of the 3D NS-α (with given admissible initial data). A similar result was later proved in the context of the 3D-Leray-α model in [44]. We note that in both of these cases, the reconstruction of the NS-α solution is done via data that do not come from the NS equations themselves, but from the NS-α model. In this sense, the observations are not physical, since the they are generated from NS-α, which itself is not physical, but is used as a regularized version of Navier-Stokes; however, [1] and [44] were important steps in demonstrating that AOT-type data assimilation is not limited to 2D equations, but can be extended to 3D in certain contexts. That is to say, the barrier to extending theoretical results about the AOT method for 2D NS to 3D NS is not directly due to the dimensionality, but rather it is likely due to the nature of the equations in three dimensions, which is the same barrier standing in the way of resolving the major open problem of global existence for strong solutions of the 3D NS equations.
It is worth noting that classical data assimilation is largely focused on statistical approaches utilizing the Kalman filter [65] and 3D/4D-Var methods, and their many variations (see, e.g., [30,66,78,83], and the references therein). The AOT algorithm, which is also referred to in the literature as continuous data assimilation, differs markedly from the Kalman filter approach. Rather than starting from the numerical level and employing statistical tools, AOT data assimilation manifests at the PDE level by way of a feedback-control term which penalizes deviations from interpolations of observable data. The use of interpolation is a key difference between the AOT method and the so-called nudging or Newtonian relaxation methods introduced in [3,55]. See, e.g., [69] for an overview of nudging methods. We note that a method that resembles the AOT method in some respects was introduced in [14] in the context of stochastic differential equations. Much recent literature has built upon the AOT algorithm and its ideas; see, e.g., [1, 2, 11-13, 20, 23, 39-43, 47, 48, 51-53, 57, 59-61, 70, 72, 85-87, 94, 96]. Computational trials of the AOT algorithm and its variants were carried out in the case of the NS equations [52,75,80], the 2D Bénard convection equations [2], and the 1D Kuramoto-Sivashinsky equations [72,85].
We now describe the main algorithm proposed in this work. First, consider the Navier-Stokes equations over the 2D periodic domain Ω = T 2 = 1 For the purposes of the present work, we take this system to be the correct underlying physical system. In practice, the initial data u 0 might not be known at every point in the domain; indeed, it might not be known at all. However, we assume that we have some information about u(t), such as measurements 2 of its values at certain points in the domain, from which we may construct an approximation of u(t), namely I h (u(t)), where I h is a linear operator based on O(h −1 ) data points, satisfying (1.3) below. Thus, we will consider a new system with an AOT-type feedback-control term . Moreover, we assume that we have some reason to consider a regularized version of the original physical system; for example, it may be that our simulation is under-resolved, leading to numerical errors. For the reasons described above, we use the Voigt regularization for this purpose. Thus, we propose the following system over the domain Ω.
We take I h to be a linear operator satisfying for some constants c 1 , h > 0, and all φ ∈ H 1 . Many common interpolants satisfy this bound. For example, piecewise linear interpolation on a mesh with minimum spacing h between grid points (see, e.g., [38], Corollary 1.109), and low-mode Fourier truncation preserving a number of modes proportional to 1/h (see, e.g., [17,98]), both satisfy (1.3). We now state our global well-posedness result regarding system (1.2) and provide in Section 3 a formal proof (using a priori estimates) of global existence of solution in H and V .
We recall the dimensionless Grashof number G, defined as Here, for simplicity, we assume the measurements to be taken at every time t, and to have perfect accuracy, but see [23] for an adaption of the AOT algorithm to the case of discrete-in-time measurements, and see [11] for adaption to the case of stochastically noisy data.
where λ 1 is the first eigenvalue of the Stokes operator A defined in Section 2. The next two theorems are our main results concerning the convergence of the data assimilation algorithm.
Let u and v be the corresponding solutions of (1.1) and (1.2), respectively. Let h > 0 and µ > 0 be such that for some constant C > 0. Then, for some t 0 > 0, and a.e. t > t 0 and any α > 0 such that In particular, Theorem 1.2 says that, for large times, the regularized solution approximates the true solution to the 2D NS equations, up to some error, and this error goes to zero as a power of α as α → 0. The next theorem provides faster rate of convergence for the H 1 norm, and also convergence in the H 2 norm, provided the data is smoother, and the nudging parameter and interpolation parameters, µ and h, satisfies stricter inequalities.
for some constant C > 0 depending on ν, G and f H 2 . Then, for some t 0 > 0, for a.e. t > 0 and any α > 0 such that the following inequality holds.

Preliminaries
All through this paper we denote the usual Lebesgue and Sobolev spaces by L p for 1 ≤ p ≤ ∞ and H s ≡ W s,2 for s > 0, respectively. Let V be the set of all divergence-free, mean-free, smooth vector fields on T 2 . We follow the standard convention of denoting by H and V the closures of V in L 2 and H 1 , respectively, with inner products respectively, associated with the norms u H = (u, u) 1/2 and u V = (∇u, ∇u) 1/2 . For the sake of convenience, we use u L 2 and u H 1 to denote the above norms in H and V , respectively.
The following materials are standard in the study of fluid dynamics, in particular for the Navier-Stokes equations and related PDEs, and we refer to reader to [28,99] for more details. We define the Stokes operator A −P σ ∆ with domain D(A) H 2 ∩ V , where P σ is the orthogonal Leray-Helmholtz projector from L 2 to H. Notice that under periodic boundary conditions, we have A = −∆P σ . Moreover, the Stokes operator can be extended as a linear operator from V to V ′ as It is well-known that A −1 : H ֒→ D(A) is a positive-definite, self-adjoint, and compact operator from H into itself, thus, A −1 possesses an orthonormal basis of positive eigenfunctions {w k } ∞ k=1 in H, corresponding to a sequence of non-increasing sequence of eigenvalues. Therefore, A has non-decreasing eigenvalues λ k , i.e., 0 ≤ λ 1 ≤ λ 2 , . . .. Then, we have the following Poincaré inqualities: for u ∈ D(A). Thus, ∇u L 2 is equivalent to u H 1 .
Also, we state the following Ladyzhenskaya inequality u 2 L 4 ≤ c u L 2 ∇u L 2 for all u ∈ V , which is a variation of the following interpolation result that is frequently used in this paper (see, e.g., [89] for a detailed proof). Assume 1 ≤ q, r ≤ ∞, and 0 < γ < 1. For v ∈ L q (T n ), such that ∂ α v ∈ L r (T n ), for |α| = m, then (2.1) Next, for any u, v ∈ V, we use the standard notation for the bilinear term which can be extended to a continuous map B : for smooth functions u, v, w ∈ V. We will use the following important properties of the map B, also known as Ladyzhenskaya-Sobolev inequalities which can be derived from the aforementioned Nirenberg inequality. Detailed proofs can be found in, e.g., [28,46]. Lemma 2.1. For the operator B, the following identities hold. See, e.g., [28,46,99] for details.
Regarding the pressure term, we recall a corollary of a deep result of deRham; namely, that for any distribution F, the equality F = ∇p holds for some distribution p if and only if F, w = 0 for all w ∈ V. See [100] for an elementary proof of the corollary.
In order to prove Theorem 1.1, we summarize a series of well-known results about 2D NS system in the next theorem and refer the readers to [28,[31][32][33][34][35]99] for more details.
Theorem 2.2. With initial data u 0 ∈ V and external forcing f ∈ H, the 2D NS system (1.1) possess a unique global solution u ∈ L ∞ (0, ∞; V ) ∩ L 2 (0, ∞; D(A)). In particular, for some t 0 and t > t 0 , we have Moreover, if f ∈ L 2 (0, T ; H), then the H 2 norm of u is uniformly bounded after certain time t 0 > 0, i.e., we have If we further assume that u 0 ∈ H 2 ∩ V and f ∈ H 2 , then the H 3 norm of u is uniformly bounded; namely, Also, we use the following Brezis-Gallouet inequality in our paper.
Remark 2.5. In fact we will use the above lemma with a slight variation, i.e., which was proved in [20] by a straight-forward minor modification of the proof in [62].

Global well-posedness of system (1.2)
In this section, we provide the proof of Theorem 1.1. Here, we work formally, but we note that following formal derivation including the a priori estimates can be made rigorous by using a Galerkin approximation then passing to the limit in an appropriate sense under certain compactness conditions by using Aubin-Lions Lemma (c.f. [28,99]). Since such approach is standard for both 2D and 3D NS system as well as other turbulence models such NS-Voigt system and NS-α model, we omit the details here but refer the readers to [19,28,43,73,76,99] and the reference therein for more detailed discussions.
Proof of Theorem 1.1.
Taking a formal inner-product of (1.2a) with v, integrating by parts, and employing (2.2b), we obtain The first term on the right side of (3.1) is bounded by Using (1.3), we estimate the second term on the right side of (3.1) by where we used the Cauchy-Schwarz, Poincaré, and Young inequalities.
Regarding the last term on the right side of (3.1), we use (1.3) to obtain Then, by combining all the above estimates, we obtain where the constant C depends on λ 1 , ν, u H 1 , f L 2 , c 1 , h, and µ. Therefore, Grönwall inequality implies for all T > 0. In order to obtain the H 2 bounds on v, we multiply (1.2a) by Av and integrate by parts to obtain 1 2 where we used (2.2b). Then, we estimate the three terms on the right side of (3.3). Using Cauchy-Schwarz inequality, the first term is bounded by By (1.3), the second term is bounded by Av 2 L 2 , while the last term is estimated similarly as Combining the above estimates, we obtain where the constant C depends on λ 1 , ν, u H 1 , f L 2 , c 1 , h, and µ. Hence, Grönwall inequality implies v ∈ L ∞ (0, T ; H 2 ∩ V ). Next, we show that the time derivative dv/dt is uniformly bounded in V . In fact, for any test function φ ∈ V , we have where we integrated by parts and used ∇ · φ = 0. Since v is bounded in V , the first term on the right side of (3.4) is bounded by

After integration by parts, the second term is bounded by
while Cauchy-Schwarz inequality implies that the third term is bounded by Regarding the last two terms, we used (1.3) and proceed as and similar argument gives Thus, we conclude that (I − α 2 ∆)dv/dt is uniformly bounded in V ′ . By inverting the Helmholtz operator (I − α 2 ∆), we obtain the uniform boundedness of dv/dt in V . By using a difference test functionφ ∈ D(A), similar estimates provide that (I −α 2 ∆)dv/dt is uniformly bounded in D(A) ′ , which in turn indicates that dv/dt is uniformly bounded in H. Therefore, together with the H 2 bounds on v, one can apply standard compactness arguments and pass to the limit and get the existence of a solution in V .
The uniqueness is proven in a similar way to standard arguments for Navier-Stokes system, see, e.g., [28,99]. Hence, for the sake of simplicity, we omit the details here. The proof of Theorem 1.1 is thus complete.

Proof of L 2 convergence rate
In this section, we provide the proof of Theorem 1.2.
Proof of Theorem 1.2.
By denoting the difference w = u − v, Π = p − q and taking the difference of (1.1a) and (1.2a), we get (4.1a) Multiply (4.1a) by w, integrate by parts over Ω, and we obtain where we used (2.2b). For the first term on the right side of (4.2), we integrate by parts and use (1.1a)in order to get Then, we estimate the four integrals on the right side of (4.3).
Regarding the first integral, we integrate by parts and proceed as After integration by parts, we bound the second integral by where we used Lemma 2.1 and Lemma 2.3. The third integral is zero due to ∇ · w = 0. As for the last integral, we integrate by parts and apply Cauchy-Schwarz inequality and obtain Next, using Hölder inequality, we bound the second term on the right side of (4.2) by where we used Lemma 2.1 and Young's inequality. Finally, the last term of (4.2) is estimated as Combining all the above estimates and denoting and Y (t) := ∇w(t) 2 L 2 , we obtain which is equivalent to where the constant C u,f depends on ν, u H 3 , u H 2 and f H 1 . Then, by our choice of h, we have By the fact that w 2 Thus, in view of Theorem 2.2, for t > t 0 , due to our choice of α and µ, from which we also infer that ξ(t) := µ 2 − C ν u 2 H 1 satisfies the condition in Lemma 2.4. Therefore, we get for all t ≥ t 0 . Thus, proof of Theorem 1.1 is complete.

Proof of H 1 convergence rate
In this section, we provide the proof of Theorem 1.3.
Proof of Theorem 1.3.
Multiply equation (4.1a) by Aw, integrate by parts over Ω, and we get 1 2 where we used (2.4) and (2.5). In order to bound the first term on the right side of (5.1), we integrate by parts and use the equation of u, and obtain Next, we estimate the four integrals on the right side of (5.2). After integration by parts, the first one is bounded by For the second integral, using integration by parts and Lemma 2.1, we have where we also used Lemma 2.3. The third integral vanished due to ∇·w = 0. As for the last integral, we integrate by parts and use Hölder inequality in order to get The second term on the right side of (5.1) is estimated as − Ω (w · ∇)w · ∆u dx ≤ C w L 4 ∇w L 4 ∆u L 2 ≤ Cλ −1/4 1 where we used Poincaré's inequality. Regarding the last term in (5.1), we have By summing up all the above estimates and denoting X(t) = ∇w(t) 2 L 2 + α 2 ∆w 2 L 2 andỸ (t) = ∆w(t) 2 L 2 , we obtain d dtX (t) +Ỹ (t) ≤C u,f α 4 + C (νλ 1 ) 1/3 u 4/3 where the time integral ofC u,f is bounded since u ∈ L 2 (0, T ; H 4 ) and f H 2 . Hence, by our choice of h, µ, and α, we have 2c 2 1 h 2 µ 2 ν < µ 2 and for t > t 0 , Thus, similar arguments together with Lemma 2.4 implỹ for all t ≥ t 0 . The proof of Theorem 1.3 is thus complete.
Remark 5.1. There seems to be a compromise between the choice of µ and α. According to Theorem 1.2, µα 2 should not be too large. However, choosing larger α provides a stronger regularizing effect to the solution of (1.2), and choosing µ larger increases the convergence rate of the data assimilation. In a future computational work, we will examine the interplay between these parameters, and their effect on the accuracy and efficiency of approximating the true solution.