On the inadequacy of nudging data assimilation algorithms for non-dissipative systems: A Computational Examination of the Korteweg de-Vries and Lorenz equations

Accurate computational prediction of nonlinear evolution systems encounters two major challenges. That of initialization and of long-time accuracy, namely, how do we initialize the system and how do we ensure that our computed solution remains accurate for long intervals of time. Data assimilation is a class of techniques that attempts to resolve both of these issues by incorporating observational data into the underlying evolutionary physical model. The purpose of data assimilation is to alleviate the reliance on knowing the complete state of the initial data and achieve synchronization of our simulated solution with the true dynamics by utilizing adequately collected sparse spatio-temporal observational data (measurements).

This work focuses on characterizing the class of systems for which one should not expect the implementation of the Azouani-Olson-Titi (AOT) continuous data assimilation algorithm, introduced in [4, 5]. Since its development, this algorithm has been the subject of much study, both analytically [1, 7, 8, 9, 10, 11, 16, 18, 21, 20, 27, 29, 31, 32, 33, 34, 35, 38, 39, 42, 41, 43, 51, 53, 54, 55, 60, 66, 67, 69, 71, 73, 80] and computationally [2, 19, 24, 26, 40, 44, 61, 62, 65, 59]. We note that the AOT algorithm has been referenced a few times in the literature under the names interpolated nudging or simply nudging, among others. In the context of classical methods of data assimilation, nudging refers to a method developed by Anthes and Hokes [3, 50] in the 1970’s for numerical weather prediction. The AOT algorithm is superficially similar to nudging, but is technically and conceptually different. In the classical algorithm of Anthes and Hokes one nudges the system only at the discrete nodal points of the observational data. On the other hand, the AOT algorithm works by constructing a spatial approximation interpolation operator based on the observational data and then nudges at every point in the vicinity of the observation nodal point using the interpolation operator. Nevertheless, for simplicity, we will refer to the AOT algorithm as “nudging” for the remainder of this paper. Notably, the AOT data assimilation algorithm was designed for dissipative evolutionary nonlinear systems motivated by, and capitalizing on, the fact that the long-term dynamics of such systems are determined by the dynamics of large spatial scales (see, e.g., [25] and the references therein). Such systems include for example the 2D Navier-Stokes equations, the 2D Rayleigh-Bernard convection problem, as well as the 1D Kuramoto-Sivashinsky equation, and some nonlinear reaction-diffusion systems, to name some (cf. [75]).

In this work, we are concerned with the application of the nudging algorithm to systems that are not known to be dissipative. For instance, the shallow water equations have long been a simplified model for the ocean to which data assimilation algorithms have been applied. Contrary to the case of dissipative dynamical systems [25], it is unknown whether these equations have finitely many determining parameters, such as finitely many Fourier modes or nodal values which determine the solution uniquely. Yet researchers still apply data assimilation to these systems, see, e.g., [14, 58, 76]. The property, for example, of finitely many determining Fourier modes [37, 56, 25] is an essential first litmus test for data assimilation to work, as it ensures the existence of a minimum length scale at which measurements and observations determine the asymptotic behavior of the system we wish to predict. Without this property, there is no guarantee that any amount of observations will be sufficient to recover the dynamics of the true solution one wishes to predict.

Another classical example of a system that is not dissipative is the 2D damped-forced incompressible Euler equations, which have long been thought to lack the property of having finitely many determining modes. Unfortunately, we cannot even test data assimilation methods computationally on the 2D Euler equations, even though they are globally well-posed. It was shown in [30], that if we consider a given solution, $u$, to the 2D incompressible Euler equations (equipped with periodic boundary conditions) that is supported over finitely many Fourier modes for all times $t\in(0,\infty)$, then $u$ is necessarily constant in time. This indicates that solutions that are evolving in time develop arbitrarily small length scales as time advances, and thus one cannot accurately simulate the long time dynamics of the equations with a fixed spatial resolution. Moreover, utilizing adaptive mesh schemes is not enough to resolve this issue, as the computer running the simulations will eventually run out of memory and be unable to adequately refine the mesh. We encounter similar issues when even attempting to simulate the damped and driven 2D incompressible Euler system, even though it is known to have a global weak attractor, which is not necessarily finite-dimensional [47].These equations are of interest as the damping term dissipates energy and enstrophy, yet it does not introduce any smoothing effects unlike the addition of viscosity. For an in-depth look at the analytical properties of this system, see e.g. [6, 22, 23, 47, 52] and the references within.

As the damped and driven Euler equations appear to be ill-suited for computational study with regard to the efficacy of data assimilation, for the remainder of this paper we will instead consider the 1D Korteweg de-Vries (KdV) equation as well as a partially dissipative variant of the Lorenz 1963 system. We consider the KdV equation in particular, as it was proven in [45, 46] that the system has a finite-dimensional global attractor in the presence of a damping and forcing terms. Moreover, it was shown in [55] that the nudging algorithm recovers the true solution given sufficiently many observed Fourier modes. Without the presence of forcing and damping, the KdV equation is a dispersive equation that has infinitely many conserved quantities, the $L^{2}$ norm being among them. The lack of dissipation and these conservation properties make this equation ideal to consider as a relatively simple 1D model to explore why data assimilation fails to recover solutions to systems without finitely many determining modes. The 1963 Lorenz system is a chaotic system that has been used to test methods of data assimilation. We note that the AOT algorithm has applied to the fully dissipative 1963 Lorenz system in [12, 72, 17].

The remainder of this manuscript is organized as follows: in Section 2 we discuss various mathematical preliminaries regarding the formulation of the KdV equation and we show that, unlike the damped and driven case, the standard version of KdV does not have the finitely many determining modes property. In Section 3 we establish the numerical framework we use to investigate the efficacy of data assimilation computationally. Section 4 is split into two subsections. In Section 4.1 we show the numerical construction of pathological examples where nudging fails to recover the true solution for the KdV equations. In Section 4.2 we verify the effectiveness of nudging on the counterexamples from Section 4.1 when the system is dissipative. In Section 5 we consider a partially dissipative variant of the Lorenz system and we show that nudging fails to recover the reference solution in this setting. We end this work in Section 6 with a discussion of the conclusions we reach from our numerical results.

The nudging algorithm is given as follows. Consider the following evolution system:

where here $F$ is a given differential operator that is potentially nonlocal and nonlinear. We will specifically be considering solutions of Equation 2.1 that lie in, e.g., $L^{2}$ subject to periodic boundary conditions.

In the remainder of this section we will utilize $H^{m}_{\text{per}}$, with $m$ being a non-negative integer, to denote the closure of the set of $L$-periodic trigonometric polynomials on $\mathbf{R}$, with respect to the standard Hilbert norm:

We note that here $H^{0}_{\text{per}}=L^{2}_{\text{per}}$, with $L^{2}_{\text{per}}$ denoting the space of $L^{2}$ functions subject to $L$-periodic boundary conditions.

We assume that Equation 2.1 is globally well-posed and for such systems we define the related nudged system, introduced in [4, 5], as follows:

Here $I_{h}$ is a linear spatial interpolant with associated length scale $h$ and $\mu$ is a constant feedback-control nudging parameter. For our purposes we will consider only $I_{h}\mathrel{\mathop{\mathchar 58\relax}}=P_{M}$, where $P_{M}$ is given to be a projection operator onto the lowest $M$ Fourier modes, i.e. $h\sim\frac{L}{M}$, where $L$ is a typical large spatial scale for the underlying evolutionary equation Equation 2.1. As we require $P_{M}u$ to be well defined we will only consider solutions to Equation 2.1 that are in $L^{2}_{\text{per}}$ for all positive time.

For the purpose of this study we will let $F$ be given such that we obtain the KdV equation:

Here $u\mathrel{\mathop{\mathchar 58\relax}}\mathbf{R}\times[0,\infty)\to\mathbf{R}$ represents the amplitude of a wave and $\delta$ is the dispersive coefficient. We consider the KdV equation equipped with the periodic boundary condition $u(x,t)=u(x+2,t)$, for all $x\in\mathbf{R}$ and all time $t\geq 0$ with spatial mean zero. Hence, the typical large length scale in this case is $L=2$. We note that choosing the length of the spatial domain to be 2 is arbitrary and was chosen as this choice of domain appears in the literature, see, e.g., [79].

The corresponding nudged system is given as follows:

Here $\mu>0$ is a feedback-control nudging parameter and $P_{M}$ is given to be the projection operator onto the lowest $M$ Fourier modes. Here $u$ is the unknown reference solution of Equation 2.6, for which the first $M$ Fourier modes are observed (given).

We note that below and through the remainder of the paper, we will use the notation $\hat{f}$ to reference the Fourier modes of an arbitrary function $f$ in $L^{2}_{\text{per}}$. That is,

For simplicity of presentation, we omit the time index and, unless otherwise stated, we write $\hat{f}_{k}\mathrel{\mathop{\mathchar 58\relax}}=\hat{f}(k,t)$.

We will now state and prove some propositions that will be useful in our numerical results.

Now, using these propositions, we show that Equation 2.6 does not have the finitely many determining modes property.

In a future work, we will extend this argument to the damped and driven Euler equations, which introduce additional complexities due to the external forcing and damping terms.

In this section we will detail the numerical methods we utilize to simulate Equations 2.6 and 2.7. All simulations depicted in this study were performed Matlab (R2023a) using our own code. For this computational study we implemented the KdV equation on $\mathbf{R}$, subject to periodic B.C. with period 2 with pseudo-spectral methods. Specifically, we implement the equations using an integrating factor method to compute the linear term exactly with a fourth-order Runge-Kutta scheme for the timestepping with the nonlinear term computed explicitly. For an overview of integrating factor schemes see e.g. [57, 77] and the references contained within. Our particular choice of numerical scheme, which is given Algorithm 1, was chosen as it allows for larger timesteps than traditional schemes for KdV, such as the Zabusky-Kruskal explicit leapfrog scheme seen in the literature [79].

We note that in the simulations of Equation 2.6, depicted below, we utilized different spatial resolutions depending on the particular choice of initial profile and dispersion coefficient, $\delta^{2}$. We ensured in our simulations that $\Delta t\lesssim{\Delta x}^{2}$, which was consistent with the implementation of this timestepping scheme in chapter 10 of [77] and appeared to be numerically stable for all spatial resolutions we utilized in this work. Here $\Delta t=t_{n+1}-t_{n}$ and $\Delta x=\frac{2}{N}$ represent the temporal and spatial resolution scales with $N$ indicating the spatial resolution, given in powers of $2$. While we do not explicitly state the spatial resolution for each simulation, one can see in the spectral plots, i.e. Figures 1, 2, 3, 4, 5 and 6, that the x-axis covers $[1,N/2]$. In this study, $N$ ranges between $2^{7}$ and $2^{11}$, depending on the amount of Fourier modes required to accurately simulate a given solution. $N$ is chosen specifically to ensure that the solution is well resolved, which here means that Fourier coefficients of the reference solution decay to machine precision before the dealiasing cutoff, i.e. the vertical red line featured in Figures 1, 2, 3, 4, 5 and 6. That is, the momentum decays exponentially in wave number, we take $N$ large enough such that the momentum decays to machine precision (given by $2.216e-16$) before the wave number $\frac{N}{3}$. We utilize 2/3 truncation for dealiasing to accurately calculate the nonlinear term, i.e., we set the last 1/3 of the Fourier modes to zero to prevent aliasing effects. Additionally, we note that while the red vertical line on spectral plots (see Figures 6(b), 6(d), 4(b) and 5(b) denotes the dealiasing cutoff, the vertical blue line denotes the observational window cutoff, that is, the Fourier modes to the left of and including the blue vertical line are accurately observed and used in the computation of the nudging term.

The nudging term is implemented numerically in Algorithm 1 as an explicit term. This implementation of the feedback-control term that we used induces a CFL constraint required for the numerical scheme to remain stable, $\mu\lesssim\frac{2}{\Delta t}$. This CFL constraint can be removed by utilizing an implicit formulation of the nudging term, which we do not consider in this work. It is worth noting that the nudging term cannot be implemented numerically using multi-stage methods, such as RK4, due to the need to interpolate in time (see e.g. [68]). The simulations featured in this work primarily use a value of $\mu=100$, which was found to be sufficient for this study. For an in-depth examination of the role the value of the nudging parameter plays in convergence, see [15].

To test the effectiveness of nudging on the KdV equation Equation 2.6 we utilized an “identical twin” experimental design. This experimental design is standard in the validation of data assimilation methods, where we first initialize a reference solution. This solution will be treated as the true solution which we will attempt to recover using observations. The observations will be given by the lowest $M$ Fourier modes. We note that the observations are considered to be perfect with no distortions (i.e., no error measurements) or stochastic noise.

We note that typically in the data assimilation literature for dissipative evolutionary equations, the equations are evolved forward in time long enough such that the solution remains in an absorbing set in a relevant phase space. That is, we run solutions long enough forward in time such that they are uniformly bounded with respect to time in $L^{2}$, $H^{1}$, or another relevant norm for all future times, see e.g., [74, 75] for a full definition of absorbing sets as well as their properties. In this study we do not evolve the equations forward in time to arbitrarily large times, as solutions are not dissipative and so there is not necessarily an absorbing set that attracts all bounded solutions after a long enough time. In fact, as the $L^{2}$ norm is a conserved quantity for Equation 2.6, such an absorbing set cannot exist in $L^{2}$. We opted to simply evolve forward the initial profile a given amount of time, typically $t=1$, in order to ensure that the reference solution has developed interesting dynamics outside the observation domain. While solutions can develop small length scales not captured by the observed modes, they typically have a finite number of active Fourier modes, with the total number of modes active at a given time dependent on the initial condition, the current time, the total momentum in the system, and the dispersion parameter, $\delta^{2}$. An active Fourier mode in this sense means that $\absolutevalue{\hat{u}_{k}}>\epsilon$, where $k\in\mathbf{Z}$ is the wavemode index and $\epsilon$ indicates machine precision.

As the focus of this work is examining the applicability of the AOT algorithm to dynamical systems that are not dissipative, we now turn to examine a nondissipative variant of the Lorenz 1963 system. Lorenz introduced in 1963 a simplified/reduced model of atmospheric convection [64], and so it is of particular interest for studying methods of data assimilation, see e.g. [49, 12, 29, 48, 70, 63] and the references within. The 1963 Lorenz system is given by:

Here $\sigma>0$ is the Prandtl number, $r>0$ is the Rayleigh number, and $b>0$ is a geometric factor. It is well-known that, for certain parameter choices (e.g., $\sigma=10$, $r=28$, $b=8/3$), that the Lorenz system has a chaotic global attractor [78].

It is also known that applying the AOT algorithm, or a synchronization scheme, as in [49, 48], by observing only the $X$ component is sufficient to synchronize the nudged solution with the reference solution (see e.g. [49, 12, 48]). In particular, we can consider the following three synchronization schemes from [48] (which was motivated by [68]) to understand when data assimilation works for this system.

The above equations, Equations 5.2, 5.3 and 5.4, correspond to the regimes when the observed values of $X$, $Y$, or $Z$, respectively, are observed and inserted directly into the equations for the remaining, unobserved components. That is, e.g. in Equation 5.2 we observe the true value of $X$ continuously in time and insert that into the evolution equations of $y$ and $z$. In these equations we note that numerically we should not be simulating the evolution of $X$, $Y$, or $Z$ as they appear in Equations 5.2, 5.3 and 5.4, respectively, but rather these values are observed (given) from the true solution in Equation 5.1 and are inserted directly into the equations for the remaining variables. Analysis on Equations 5.2, 5.3 and 5.4 was conducted in [48] where it was shown that the continuous coupling on $X$ or $Y$ alone was sufficient to obtain synchronization between the assimilated equations and the reference solution, however the coupling seen in Equation 5.4 was not sufficient to obtain convergence in either the $x$ or $y$ components to the corresponding components $X$ or $Y$ of the true solution. We do note that in [48] the equations for Equation 5.1 were written in a slightly different form, however this can be easily reconciled by using the transformation $\tilde{Z}=Z+(r+\sigma)$.

For the sake of completeness, we will discuss the convergence results for the above systems, Equations 5.2, 5.3 and 5.4 proven in [48] and discuss how the results break down in the absence of dissipation. Let

where $(X,Y,Z)$ is the given solution of Equation 5.1 and $(x,y,z)$ is the solution to one of the assimilated equation, Equations 5.2, 5.3 and 5.4. Here we define $\delta$ with the understanding that the variable in $(x,y,z)$ missing from the assimilated equation is taken to be equal to the corresponding variable in Equation 5.1, e.g. when considering Equation 5.2 we define $x(t)\mathrel{\mathop{\mathchar 58\relax}}=X(t)$ for all time, $t$ and have $\delta=(0,Y-y,Z-z)$. The convergence results for the synchronization algorithms, Equations 5.2, 5.3 and 5.4 are given as follows:

The proofs of 5.0.1 are derived by taking the difference of the assimilated equation (Equation 5.2 or Equation 5.3) with Equation 5.1, summing the squares of the derivatives, and applying Grönwall’s inequality to the resulting quantity. The $J$ estimate in Equation 5.7 arises from the use of the global bound

which was derived in [28] and was used to bound the contribution from the nonlinear term in the $z$ equation.

We note that the analysis done in the proof of 5.0.1 is dependent on the particular parameters and requires $\sigma>0$, $b>0$ and $r>0$. To illustrate how these sorts of estimates break down in the absence of dissipation, we consider the case when both $X$ and $Y$ are observed, i.e., are given from the true solution of Equation 5.1. This leads to the following system of equations

Now, if we take $w=z-Z$, where $Z$ solves Equation 5.1, we see from Equations 5.1 and 5.10 that $w$ satisfies the following equation

From here, we can see that if $b>0$, we can multiply both sides by $w$ and apply Grönwall’s inequality to obtain

which yields convergence of $z$ to $Z$ exponentially fast in time. If we instead consider $b=0$, removing the dissipation on the $Z$ component of Equation 5.1, then we now have that $\dot{w}=0$, and so $w(t)=z_{0}-Z_{0}$ holds for all time. Therefore, we have that the absence of dissipation on the $Z$ component prevents synchronization between the assimilated solution and the reference solution, even when both the $X$ and $Y$ components are observed continuously in time.

We note that the above discussion has been on the synchronization data assimilation scheme, given by the direct insertion, as it was first introduced in [68], of the observed quantities into the model equations. We now consider the nudged Lorenz system for the observed $X$ and $Y$ components of Equation 5.1, given as follows:

Here $\mu>0$ is the nudging coefficient, which controls the strength of the feedback-control terms. For simplicity, we utilize the same nudging parameter, $\mu$, in the equations for $x$ and $y$, however this is not strictly necessary.

At any fixed time, $t$, the solutions, $(x(t;\mu),y(t;\mu),z(t;\mu))$, to the nudged equations, Equation 5.13, are expected to converge to the solution $(X(t),Y(t),z(t))$ of the synchronization algorithm, given by Equation 5.10 as $\mu\to\infty$. This can be justified rigorously by applying the argument present in [15] to the Lorenz 1963 system, provided both methods are initialized with the same initial data. An analogous statement was made in [72] for the Lorenz system in particular (with $b>0$), where it was shown that the solution, $x(t;\mu)$, of Equation 5.13, converges to $X(t)$ at any fixed time, $t$, as we send $\mu\to\infty$ when the nudging term is present only on the equation for $x$ in Equation 5.13. In the case where $b=0$, if we assume a priori that the solutions of Equation 5.13 converge to those of Equation 5.10 as $\mu\to\infty$, one can deduce using the reverse triangle inequality, that $\norm{\delta(t)}^{2}_{\ell^{2}}\geq|z_{0}-Z_{0}|^{2}$, where $\delta$ is given as defined previously in Equation 5.5. Therefore, we expect the nudged Lorenz system, Equation 5.13, to fail to recover the true solution in the same way that Equation 5.10 fails when dissipation is removed from the equations for $z$.

We now demonstrate numerically that this dissipation mechanism (when $\sigma,b>0$) of the Lorenz system is necessary for data assimilation methods to work. In particular, we will be considering a partially dissipative variant of the Lorenz system by taking $b=0$, i.e., eliminating the dissipation mechanism only on the $Z$ component. The analytical results regarding the nudging algorithm applied to the Lorenz equations have been done when $b>0$ and are qualitatively similar to those we reference above in 5.0.1 from [48], see e.g. [12, 72, 17]. That is, the solutions of Equation 5.13 synchronize with those of Equation 5.1 exponentially fast in time for large enough nudging parameter values and with the exponential rate depending on the particular parameters of Equation 5.1. While nudging is only required on either the $X$ or $Y$ components in order to drive convergence, we will utilize nudging on both the $X$ and $Y$ components to emphasize that, even with additional nudging, the method is not effective when the dissipation on the $z$-component is eliminated (i.e. when $b=0$).

We have numerically implemented the Lorenz system as in [12, 72] using an explicit Euler scheme for the time discretization. As in [72], we will utilize $\Delta t=1e-4$ for our timestep. We note that in many works the reference solution is initialized by running the solution forward in time until it is approximately on the global attractor (see e.g. [12, 49]), however as we are attempting to test nudging for non-dissipative or partially dissipative systems, which need not possess global attractors, we will instead be utilizing the initial data $(X,Y,Z)=(30,40,50)$ as in [72].

First, we verify the validity of nudging on the Lorenz system with full dissipation. Using the initial data and parameter choices discussed previously we apply nudging to the Lorenz system with $b=8/3$. The results can be seen in Figure 7, where we have applied nudging to only the first two components of Equation 5.13. Here we can see that all variables converge to an error level of approximately $1e-12$ with spurious dips below this error level in each component. Note that here we see exponential convergence to the true solution, as we expect in this case.

Next, we examine what happens when we apply nudging for the partially dissipative Lorenz system, i.e., with $b=0$. The results can be seen in Figure 8. Note that here we still see exponential decay to the true solution in both $x$ and $y$. While nudging appears to affect $z$ initially, this effect quickly fades after a short interval of time and the error appears to stabilize to be approximately constant $1e+1$. The error level that is attained appears to be directly linked to the difference in the initial conditions for the $z$ component. To illustrate this, we plotted the $\ell^{2}$ (Euclidean norm) error over time when both the reference and nudged equations are initialized with the same initial data, except for the $z$ component. Here $Z(t_{0})=50$ and $z(t_{0})=50+\delta$, where $\delta$ ranges in values from $1e-1$ to $1e-12$. The results of these tests can be seen in Figure 9 and they illustrate that the $\ell^{2}$ error is highly sensitive to small perturbations in the initial condition for $z$. We note that the behavior seen in Figure 9 appears to align with what one should expect from our earlier discussion of Equation 5.11 for the synchronization algorithm. We note that we obtain an initial exponential decay in the error in the $z$ component which is driven by the convergence in both $x$ and $y$. This initial convergence levels off and becomes approximately constant as $x$ and $y$ converge to $X$ and $Y$, respectively.

Based on the results of the simulation presented in Section 4, we conclude that the AOT data assimilation algorithm is unsuitable for non-dissipative or partially dissipative dynamical systems, such as the KdV equation, Equation 2.6, and the partially dissipative Lorenz system, Equation 5.1. In fact, we conclude that the pathologic counterexamples used in Section 4.1 present an insurmountable challenge for any method of data assimilation to recover. In particular, we note that counterexample we utilize in Figure 1 will fail for any method of data assimilation, as the true solution has no discernible impact on the observed modes, and so it is impossible to reliably tell the true solution apart from the zero solution (or one of an infinite number of solutions one can readily construct) based solely on the observations. Moreover, the counter examples considered in Section 4.1 are not specific to the KdV equation. The particular solutions considered in 2.0.5 were constructed due to the periodic invariance induced by the convolutional Fourier mode structure of the nonlinear term (see 2.0.4), and thus one can construct similar counter examples for other non-dissipative systems with this same nonlinearity, such as the inviscid Burgers equations.

In addition to showing that data assimilation simply doesn’t work for certain pathological counterexamples, we note that the extent to which data assimilation does work for the classical KdV equation may be highly dependent on the choice of initial data for the nudged system, $v^{in}$. Indeed, while it was proved in [4] that for 2D NSE the choice of initial data $v^{in}$ was arbitrary, we have here that data assimilation will simply fail if not initialized properly, as it does in Figure 1. The problem of initialization is one of the core strengths of the nudging algorithm, and so it is worth noting that this property fails to hold for a system without finitely many determining modes.

We note that data assimilation could be made to work if one were to assume more of the true solution. For instance, if one were to assume that the solution was supported on an annulus in Fourier space, one could eliminate all the examples that we study here. We caution anyone against doing this, however, as imposing these sorts of restrictions on the solution can trivialize the dynamics of equations in a way that may not be physical.