PyDA: A Hands-On Introduction to Dynamical Data Assimilation with Python

Abstract

Dynamic data assimilation offers a suite of algorithms that merge measurement data with numerical simulations to predict accurate state trajectories. Meteorological centers rely heavily on data assimilation to achieve trustworthy weather forecast. With the advance in measurement systems, as well as the reduction in sensor prices, data assimilation (DA) techniques are applicable to various fields, other than meteorology. However, beginners usually face hardships digesting the core ideas from the available sophisticated resources requiring a steep learning curve. In this tutorial, we lay out the mathematical principles behind DA with easy-to-follow Python module implementations so that this group of newcomers can quickly feel the essence of DA algorithms. We explore a series of common variational, and sequential techniques, and highlight major differences and potential extensions. We demonstrate the presented approaches using an array of fluid flow applications with varying levels of complexity.

1. Introduction

Data assimilation (DA) refers to a class of techniques that lie at the interface between computational sciences and real measurements, and aim at fusing information from both sides to provide better estimates of the system’s state. One of the very mature applications that significantly utilize DA is weather forecast, that we rely on in our daily life. In order to predict the weather (or the state of any system) in the future, a model has to be solved, most often by numerical simulations. However, a few problems rise at this point and we refer to only few of them here. First, for accurate predictions, these simulations need to be initiated from the true initial condition, which is never known exactly. For large scale systems, it is almost impossible to experimentally measure the full state of the system at a given time. For example, imagine simulating the atmospheric or oceanic flow, then you need to measure the velocity, temperature, density, etc. at every location corresponding to your numerical grid! Even in the hypothetical case when this is possible, measurements are always contaminated by noise, reducing the fidelity of your estimation. Second, the mathematical model that completely describes all the underlying processes and dynamics of the system is either unknown or hard to deal with. Then, approximate and simplified models are adopted instead. Third, the computational resources always constrain the level of accuracy in the employed schemes and enforce numerical approximations. Luckily, DA appears at the intersection of all these efforts and introduces a variety of approaches to mitigate these problem, or at least reduce their effects. In particular, DA techniques combines possibly incomplete dynamical models, prior information about initial system’s state and parameterization, and sparse and corrupted measurement data to yield optimized trajectory in order to describe the system’s dynamics and evolution.

As indicated above, dynamical data assimilation techniques have a long history in computational meteorology and geophysical fluid dynamics sciences [1,2,3]. Then comes the question: why do we write such an introductory tutorial about a historical topic? Before answering this question, we highlight a few points. DA borrows ideas from numerical modeling and analysis, linear algebra, optimization, and control. Although these topics are taught separately in almost every engineering discipline, their combination is rarely presented. We believe that incorporating DA course in engineering curricula is important nowadays as it provides a variety of global tools and ideas that can potentially be applied in many areas, not just meteorology. This was proven while administering a graduate class on “Data Assimilation in Science and Engineering” at Oklahoma State University, as students from different disciplines and backgrounds were astonished by the feasibility and utility of DA techniques to solve numerous inverse problems they are working on, not related to weather forecast. Nonetheless, the availability of beginner-friendly resources has been the major shortage that students suffered from. The majority of textbooks either derives DA algorithms from their very deep roots or surveys their historical developments, without focus on actual implementations. On the other hand, available packages are presented in a sophisticated way that optimizes data storage and handling, computational cost, and convergence. However, this level of sophistication takes a steep learning curve to understand the computational pipeline as well as the algorithmic steps, and a lot of learners fall hopeless during this journey.

Therefore, the main objective of this tutorial paper is to familiarize beginning researchers and practitioners with basic DA ideas along with easy-to-follow pieces of codes to feel the essence of DA and trigger the priceless “aha” moments. With this in mind, we choose Python as the coding language, being a popular, interpreted language, and easy to understand even with minimum programming background, although not the most computationally favored language in high performance computing (HPC) environments. Moreover, whenever possible, we utilize the built-in functions and libraries to minimize the user coding efforts. In other words, the provided codes are presented for demonstrative purposes only using an array of academic test problems, and significant modifications should be incorporated before dealing with complex applications. Meanwhile, this tutorial will give the reader a jump-start that hopefully shortens the learning curve of more advanced packages. A Python-based DA testing suite has been also designed to compare different methodologies [4].

Dynamical data assimilation techniques can be generally classified into variational DA and sequential DA. Variational data assimilation which works by setting an optimization problem defined by a cost functional along with constraints that collectively incorporate our knowledge about the system. The minimizer of the cost functional represents the DA estimate of the unknown system’s variables and/or parameters. On the other hand, in sequential methods (also known as statistical methods), the system state is evolved in time using background information until observations become available. At this instant, an update (correction) to the system’s variables and/or parameters is estimated and the solver is re-initialized with this new updated information until new measurements are collected, and so on. We give an overview of both approaches as well as basic implementation. In particular, we briefly discuss the three dimensional variational data assimilation (3DVAR) [5,6], the four dimensional variational data assimilation (4DVAR) [7,8,9,10,11,12], and forward sensitivity method (FSM) [13,14] as examples of variational approaches. Kalman filtering and its variants [15,16,17,18,19,20,21,22] are the most popular applications of sequential methods. We introduce the main ideas behind standard Kalman filter and its extensions for nonlinear and high-dimensional problems. The famous Lorenz 63 is utilized to illustrate the merit of all presented algorithms, being a simple low-order dynamical system that exhibit interesting dynamics. This is to help readers to digest the different pieces of codes and follow the computational pipeline. Then, the paper is concluded with a section that provides the deployment of selected DA approaches for dynamical systems with increasing levels of dimensionality and complexity. We highlight here that the primary purpose of this paper is to provide an introductory tutorial on the data assimilation for educational purposes. All Python implementations of the presented algorithms as well as the test cases are made publicly accessible at our GitHub repository https://github.com/Shady-Ahmed/PyDA.

2. Preliminaries

2.1. Notation

Before we dive into the technical details of dynamical data assimilation approaches, we briefly present and describe our notations and assumptions. In general, we assume that all vector-valued functions or variables are written as a column vector. Unless stated otherwise, boldfaced lowercase letters are used to denote vectors and boldface uppercase letters are reserved to matrices. We suppose that the system state at any time t is denoted as $\mathbf{u}\left(t\right)={[u_1\left(t\right),u_2\left(t\right),\dots ,u_n\left(t\right)]}^T\in {\mathbb{R}}^n$, where n is the state-space dimension. The dynamics of the system are governed by the following differential equation

(1)$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}}=\mathbf{f}(\mathbf{u};\theta ),\end{array}$

where $\mathbf{f}:{\mathbb{R}}^n\times {\mathbb{R}}^p\to {\mathbb{R}}^n$ encapsulates the model’s dynamics, with $\theta \in {\mathbb{R}}^p$ being the vector of model’s parameters and p being the number of these parameters. With a time-integration scheme applied, the discrete-time model can be written as follows,

(2)$\begin{array}{c}\hfill \mathbf{u}\left(t_{k+1}\right)=M(\mathbf{u}\left(t_k\right);\theta ),\end{array}$

where M is the one-time step transition map that evolves the state at time $t_k$ to time $t_{k+1}=t_k+\mathsf{\Delta }t$, with $\mathsf{\Delta }t$ being the time step length.

We denote the true value of the state variable as ${\mathbf{u}}_t$, which is assumed to be unknown and a good approximation of it is sought. Our prior information about the state $\mathbf{u}$ is called the background, with a subscript of b as ${\mathbf{u}}_b$. This represents our beginning knowledge, which might come from historical data, numerical simulations, or just an intelligent guess. The discrepancy between this background information and true state is denoted as ${\xi }_b={\mathbf{u}}_t-{\mathbf{u}}_b$, resulting from imperfect model, inaccurate model’s initialization, incorrect parameterization, numerical approximations, etc. From probabilistic point of view, we suppose that the background error has a zero mean and a covariance matrix of $\mathbf{B}$. This can be represented as $E\left[{\xi }_b\right]=0$ and $E\left[{\xi }_b{\xi }_b^T\right]=\mathbf{B}$, where $\mathbf{B}\in {\mathbb{R}}^{n\times n}$ is a symmetric and positive-definite matrix and the superscript T refers to the transpose operation. Moreover, we assume that unknown true state has a multivariate Gaussian distribution with a mean ${\mathbf{u}}_b$ and a covariance matrix $\mathbf{B}$ (i.e., ${\mathbf{u}}_t=\mathcal{N}({\mathbf{u}}_b,\mathbf{B})$).

We define the set of the collected measurements at a specific time $t_k$ as $\mathbf{w}\left(t_k\right)\in {\mathbb{R}}^m$, where m is the dimension of observation-space. We highlight that the observed quantity need not be the same as the state variable. For instance, if the state variable that we are trying to resolve is the temperature of sea surface, we may have access only to radiance measurements by satellites. However, those case be related to each other through Planch–Stefan’s law, for instance. Formally, we can relate the observables and the state variables as

(3)$\begin{array}{c}\hfill \mathbf{w}\left(t_k\right)=h\left({\mathbf{u}}_t\left(t_k\right)\right)+{\xi }_m,\end{array}$

where $h:{\mathbb{R}}^n\to {\mathbb{R}}^m$ defines the mapping from state-space to measurement-space and ${\xi }_m\in {\mathbf{R}}^m$ denotes the measurement noise. The mapping h can refer to the sampling (and probably interpolation) of state variables at the measurements locations, relating different quantities of interest (e.g., relating see surface temperature to emitted radiance), or both! The model’s map M and observation operator h can be linear, nonlinear, or a combination of them. Similar to the background error, the observation noise ${\xi }_m$ is assumed to possess a multivariate normal distribution, with a zero mean and a covariance matrix $\mathbf{R}\in {\mathbb{R}}^{m\times m}$, i.e., ${\xi }_m=\mathcal{N}(0,\mathbf{R})$. An extra grounding assumption is that the measurement noise and the state variables (either true or background) are uncorrelated. Furthermore, all noises are assumed to be temporally uncorrelated (i.e., white noise). Even though we consider only Gaussian distribution for the background error, and observation noise, we emphasize that there has been a lot of studies dealing with non-Gaussian data assimilation [23,24,25,26].

The objective of data assimilation is to provide an algorithm that fuses our prior information ${\mathbf{u}}_b$ and measurement data $\mathbf{w}$ to yield a better approximation of the unknown true state. This better approximation is called the analysis, and denoted as ${\mathbf{u}}_a$. The difference between this better approximation and the true state is denoted as ${\xi }_a={\mathbf{u}}_t-{\mathbf{u}}_a$.

2.2. Twin Experiment Framework

In a realistic situation, the true state values are unknown and noisy measurements are collected by sensing devices. However, for testing ideas, the ground truth need to be known beforehand such that the convergence and accuracy of the developed algorithm can be evaluated. In this sense, the concept of twin experiment has been popular in data assimilation (and inverse problems, in general) studies. First, a prototypical test case (all called toy problems!) is selected based on the similarities between its dynamics and real situations. Similar to your first “Hello World!” program, the Lorenz 63 and Lorenz 96 are often used in numerical weather forecast investigations, the one-dimensional Burgers equation is explored in computational fluid dynamics developments, the two-dimensional Kraichnan turbulence and three-dimensional Taylor–Green vortex are analyzed in turbulence studies, and so on. A reference true trajectory is computed by fixing all parameters and running the forward solver until some final time is reached. Synthetic measurements are then collected by sampling the true trajectory at some points in space and time. A mapping can be applied on the true state variables and arbitrary random noise is artificially added (e.g., a white Gaussian noise). Finally, the data assimilation technique of interest is implemented starting from false values of the state variables or the model’s parameters along with the synthetic measurement data. The output trajectory of the algorithm is thus compared against the reference solution, and the performance can be evaluated. It is always recommended that researchers get familiar with twin experiment frameworks as they provide well-structured and controlled environments for testing ideas. For instance, the influence of different measurement sparsity and/or level of noise can be cheaply assessed, without the need to locate or modify sensors.

3. Three Dimensional Variational Data Assimilation

The three dimensional variational data assimilation (3DVAR) framework can be derived from either an optimal control or Bayesian analysis points of view. The interested readers can be referred to other resources for mathematical foundations (e.g., [27]). In order to compute a good approximation of the system state, the following cost functional can be defined,

(4)$\begin{array}{c}\hfill J\left(\mathbf{u}\right)={\displaystyle \frac{1}{2}}{(\mathbf{w}-h\left(\mathbf{u}\right))}^T{\mathbf{R}}^{-1}(\mathbf{w}-h\left(\mathbf{u}\right))+{\displaystyle \frac{1}{2}}{(\mathbf{u}-{\mathbf{u}}_b)}^T{\mathbf{B}}^{-1}(\mathbf{u}-{\mathbf{u}}_b),\end{array}$

where the first term penalizes the discrepancy between the actual measurement $\mathbf{w}$ and the state variable mapped into the observation space $h\left(\mathbf{u}\right)$ (also called the model predicted measurement). The second term aims at incorporating the prior information, weighted by the inverse of the covariance matrix to reflect our confidence in this background. We highlight that all terms in Equation (4) are evaluated at the same time, and thus the 3DVAR can be referred to as a stationary case.

The minimizer of $J\left(\mathbf{u}\right)$ (i.e., the analysis) can be obtained by setting the gradient of the cost functional to zero as follows,

(5)$\begin{array}{c}\hfill \mathsf{\nabla }J\left(\mathbf{u}\right)=-{\mathbf{D}}_h^T\left({\mathbf{u}}_a\right){\mathbf{R}}^{-1}(\mathbf{w}-h\left({\mathbf{u}}_a\right))+{\mathbf{B}}^{-1}({\mathbf{u}}_a-{\mathbf{u}}_b)=0,\end{array}$

where ${\mathbf{D}}_h\left(\mathbf{u}\right)\in {\mathbb{R}}^{m\times n}$ is the Jacobian matrix of the operator $h\left(\mathbf{u}\right)$. The difficulty of solving Equation (5) depends on the form of $h\left({\mathbf{u}}_{\mathbf{a}}\right)$ as it can either be linear, or highly nonlinear.

3.1. Linear Case

For linear observation operator (i.e., $h\left(\mathbf{u}\right)=\mathbf{H}\mathbf{u}$, and ${\mathbf{D}}_h\left(\mathbf{u}\right)=\mathbf{H}$, where $\mathbf{H}$ is an $m\times n$ matrix), the evaluation of the analysis ${\mathbf{u}}_a$ in Equation (5) reduces to solving the following linear system of equations

(6)$\begin{array}{c}\hfill ({\mathbf{B}}^{-1}+{\mathbf{H}}^T{\mathbf{R}}^{-1}\mathbf{H}){\mathbf{u}}_a=({\mathbf{B}}^{-1}{\mathbf{u}}_b+{\mathbf{H}}^T{\mathbf{R}}^{-1}\mathbf{w}).\end{array}$

We note that $({\mathbf{B}}^{-1}+{\mathbf{H}}^T{\mathbf{R}}^{-1}\mathbf{H})$ on the left-hand side is an $n\times n$ matrix, and hence this is called the model-space approach to 3DVAR. Furthermore, a popular incremental form can be derived from Equation (6) by adding and subtracting ${\mathbf{H}}^T{\mathbf{R}}^{-1}\mathbf{H}{\mathbf{u}}_a$ to/from the right-hand side and rearranging to get the following form,

(7)$\begin{array}{c}\hfill {\mathbf{u}}_a={\mathbf{u}}_b+{({\mathbf{B}}^{-1}+{\mathbf{H}}^T{\mathbf{R}}^{-1}\mathbf{H})}^{-1}{\mathbf{H}}^T{\mathbf{R}}^{-1}(\mathbf{w}-\mathbf{H}{\mathbf{u}}_b).\end{array}$

Moreover, the Sherman–Morrison–Woodbury inversion formula can be used to derive an observation-space solution to the 3DVAR problem (for details, see [27], page 327) as follows,

(8)$\begin{array}{c}\hfill {\mathbf{u}}_a={\mathbf{u}}_b+\mathbf{B}{\mathbf{H}}^T{(\mathbf{R}+\mathbf{H}\mathbf{B}{\mathbf{H}}^T)}^{-1}(\mathbf{w}-\mathbf{H}{\mathbf{u}}_b).\end{array}$

Note that $(\mathbf{R}+\mathbf{H}\mathbf{B}{\mathbf{H}}^T)$ is an $m\times m$ matrix, compared to $({\mathbf{B}}^{-1}+{\mathbf{H}}^T{\mathbf{R}}^{-1}\mathbf{H})$ being an $n\times n$ matrix. Thus, Equation (7) or Equation (8) might be computationally favored based on the values of n and m. We also highlight that in either cases, matrix inversion is rarely (almost never) computed directly, and efficient linear system solvers should be utilized, instead. An example of a Python function for the implementation of the 3DVAR algorithm with a linear operator is shown in Listing 1.

3.2. Nonlinear Case

On the other hand, if $h\left(\mathbf{u}\right)$ is a nonlinear function, Equation (5) implies the solution of a system of nonlinear equations. Unlike linear systems, few algorithms are available to directly solve nonlinear systems and their convergence and stability are usually questionable. Alternatively, we can use Taylor series to expand $h\left(\mathbf{u}\right)$ around an initial estimate of ${\mathbf{u}}_a$, denoted as ${\mathbf{u}}_c$, where ${\mathbf{u}}_a={\mathbf{u}}_c+\mathsf{\Delta }\mathbf{u}$. The first-order approximation of $h\left({\mathbf{u}}_a\right)$ can be written as

(9)$\begin{array}{c}\hfill h\left({\mathbf{u}}_a\right)\approx h\left({\mathbf{u}}_c\right)+{\mathbf{D}}_h\left({\mathbf{u}}_c\right)\mathsf{\Delta }\mathbf{u},\end{array}$

and Equation (5) can be approximated as

(10)$\begin{array}{c}\hfill {\mathbf{D}}_h^T\left({\mathbf{u}}_c\right){\mathbf{R}}^{-1}(\mathbf{w}-h\left({\mathbf{u}}_c\right)-{\mathbf{D}}_h\left({\mathbf{u}}_c\right)\mathsf{\Delta }\mathbf{u})={\mathbf{B}}^{-1}({\mathbf{u}}_c+\mathsf{\Delta }\mathbf{u}-{\mathbf{u}}_b).\end{array}$

Thus, the correction to the initial guess of ${\mathbf{u}}_a$ can be computed by solving the following system of linear equations

(11)$\begin{array}{c}\hfill \left({\mathbf{B}}^{-1}+{\mathbf{D}}_h^T\left({\mathbf{u}}_c\right){\mathbf{R}}^{-1}{\mathbf{D}}_h\left({\mathbf{u}}_c\right)\right)\mathsf{\Delta }\mathbf{u}=\left({\mathbf{B}}^{-1}({\mathbf{u}}_b-{\mathbf{u}}_c)+{\mathbf{D}}_h{\left({\mathbf{u}}_c\right)}^T{\mathbf{R}}^{-1}(\mathbf{w}-h\left({\mathbf{u}}_c\right))\right),\end{array}$

and a new guess of ${\mathbf{u}}_a$ is estimated as ${\mathbf{u}}_c+\mathsf{\Delta }\mathbf{u}$, which is then plugged back into Equation (11) and the computations are repeated until convergence is reached. Python implementation of the 3DVAR in nonlinear observation operator is presented in Listing 2. Although we only present the first order approximation of $h\left(\mathbf{u}\right)$, higher order expansions can be utilized for increased accuracy [27].

3.3. Example: Lorenz 63 System

The Lorenz 63 equations have been utilized as a toy problem in data assimilation studies, capturing some of the interesting mechanisms of weather systems. The three-equation model can be written as

(12)$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}& =\sigma (y-x),\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}& =x(\rho -z)-y,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}}& =xy-\beta z,\hfill \end{array}$

where the values of $\sigma =10,\beta =8/4,\rho =28$ are usually used to exhibit a chaotic behavior. If we like to put Equation (12) with the notations introduced in Section 2, we can write $\mathbf{u}={[x,y,z]}^T$ with $n=3$, and $\theta ={[\sigma ,\beta ,\rho ]}^T$ with $p=3$. A Python function describing the dynamics of the Lorenz 63 system is given in Listing 3.

Equation (12) describe the continuous-time evolution of the Lorenz system. In order to obtain the discrete-time mapping $M(·;·)$, a temporal integration scheme has to be applied. In Listing 4, one-step time integration functions are provided in Python using the first-order Euler and the fourth-order Runge–KuKutta schemes. Note that these functions requires a right-hand side function as input, this is mainly the continuous-time model $f\left(\mathbf{u}\right)$ (e.g., Listing 3).

For twin experiment testing, we suppose a true initial condition of ${\mathbf{u}}_t\left(0\right)={[1,1,1]}^T$ and measurements are collected each $0.2$ time units for a total time of 2. We suppose that we measure the full system state (i.e., $h\left(\mathbf{u}\right)=\mathbf{u}$, $m=3$, and $\mathbf{H}={\mathbf{I}}_3$, where ${\mathbf{I}}_3$ is the $3\times 3$ identity matrix). Measurements are considered to be contaminated by a white Gaussian noise with a zero mean and a covariance matrix $\mathbf{R}=\mathrm{Diag}({\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2)$. For simplicity, we let ${\sigma }_1={\sigma }_2={\sigma }_3=0.15$. For data assimilation testing, we assume that we begin with a perturbed initial condition of $\mathbf{u}\left(0\right)={[2,3,4]}^T$. Then, background state values are computed at $t=0.2$ by time integration of Equation (12) starting from this false initial condition. Observations at $t=0.2$ are assimilated to provide the analysis at $t=0.2$. After that, background state values are computed at $t=0.4$ by time integration of Equation (12) starting from the analysis at $t=0.2$, and so on. A sample implementation of the 3DVAR framework is presented in Listing 5, where a fixed background covariance matrix $\mathbf{B}=\mathrm{Diag}(0.01,0.01,0.01)$ is assumed. Solution trajectories are presented in Figure 1 for a total time of 10, where observations are only available up to $t=2$.

4. Four Dimensional Variational Data Assimilation

We highlighted in Section 3 that the 3DVAR can be referred to as a stationary case since the observations, background, and analysis all correspond to a fixed time instant. In other words, the optimization problem that minimizes Equation (4) takes place in the spatial state-space only. As an extension, the four dimensional variational data assimilation (4DVAR) aims to solve the optimization in both space and time, proving a non-stationary framework. In particular, the model’s dynamics are incorporated into the optimization problem to relate different points in time to each other. The cost functional for the 4DVAR Can be written as follows,

(13)$\begin{array}{c}\hfill J\left(\mathbf{u}\left(t_0\right)\right)=\sum _{t_k\in \mathcal{T}}{\displaystyle \frac{1}{2}}{(\mathbf{w}\left(t_k\right)-h\left(\mathbf{u}\left(t_k\right)\right))}^T{\mathbf{R}}^{-1}\left(t_k\right)(\mathbf{w}\left(t_k\right)-h\left(\mathbf{u}\left(t_k\right)\right)),\end{array}$

where $\mathbf{w}\left(t_k\right)$ is the measurement at time $t_k$ and $\mathcal{T}$ defines the set of time instants where observations are available. Note that the argument of this cost functional is the initial condition $\mathbf{u}\left(t_0\right)$. In other words, the purpose of the 4DVAR algorithm is to evaluate an initial state estimate, which if evolved in time, would produce a trajectory that is as close to the collected measurements as possible (weighted by the inverse of the covariance matrix of interfering noise). This is the place where the model’s dynamics comes into play to relate initial condition to future predictions when measurements are accessible. In other words, the values of $\mathbf{u}\left(t_k\right))$ are constrained by the underlying model. Instead of presenting the linear and nonlinear mappings separately, we will focus on the general case of both nonlinear model mapping and nonlinear observation operator, where simplification to linear cases should be straightforward.

In Equation (2), we introduced the one-step transition map and here, we can extend it to the k-step transition case by applying Equation (2) recursively as

(14)$\begin{array}{c}\hfill \mathbf{u}\left(t_k\right)=M^{\left(k\right)}(\mathbf{u}\left(t_0\right);\theta )=M(M^{(k-1)}\left(\mathbf{u}\left(t_0\right);\theta );\theta \right),\end{array}$

where $M^{\left(1\right)}(\mathbf{u}\left(t_0\right);\theta )=M(\mathbf{u}\left(t_0\right);\theta )$. Now, we consider a base trajectory given by $\mathbf{ū}\left(t_k\right)$ for $k=1,2,\dots$ generated from an initial condition of $\mathbf{ū}\left(t_0\right)$. A perturbed trajectory ($\mathbf{u}\left(t_k\right)$ for $k=1,2,\dots ,$) can be obtained by correcting the initial condition as $\mathbf{u}\left(t_0\right)=\mathbf{ū}\left(t_0\right)+\mathsf{\Delta }{\mathbf{u}}_0$ and the difference between the perturbed and based trajectories can be written as

(15)$\begin{array}{c}\hfill \mathbf{u}\left(t_k\right)-\mathbf{ū}\left(t_k\right)=M^{\left(k\right)}(\mathbf{ū}\left(t_0\right)+\mathsf{\Delta }{\mathbf{u}}_0;\theta )-M^{\left(k\right)}(\mathbf{ū}\left(t_0\right);\theta ).\end{array}$

A first-order Taylor expansion of $M^{}(\mathbf{ū}\left(t_0\right)+\mathsf{\Delta }{\mathbf{u}}_0;\theta )$ around $\mathbf{ū}\left(t_0\right)$ can be given as follows

(16)$\begin{array}{c}\hfill M^{}(\mathbf{ū}\left(t_0\right)+\mathsf{\Delta }{\mathbf{u}}_0;\theta )\approx M^{}(\mathbf{ū}\left(t_0\right);\theta )+{\mathbf{D}}_M\left(\mathbf{ū}\left(t_0\right)\right)\mathsf{\Delta }{\mathbf{u}}_0,\end{array}$

where ${\mathbf{D}}_M\left(\mathbf{u}\left(t_k\right)\right)$ is the Jacobian of the model $M(\mathbf{u};\theta )$, evaluated at $\mathbf{u}\left(t_k\right)$, also known as the tangent linear operator. Note that $M(\mathbf{ū}\left(t_0\right);\theta )=\mathbf{ū}\left(t_1\right)$ and $M(\mathbf{ū}\left(t_0\right)+\mathsf{\Delta }{\mathbf{u}}_0;\theta )=\mathbf{u}\left(t_1\right)$, thus $\mathsf{\Delta }{\mathbf{u}}_1=\mathbf{u}\left(t_1\right)-\mathbf{ū}\left(t_1\right)\approx {\mathbf{D}}_M\left(\mathbf{ū}\left(t_0\right)\right)\mathsf{\Delta }{\mathbf{u}}_0$. Similarly, we can expand $M^{}(\mathbf{ū}\left(t_1\right)+\mathsf{\Delta }{\mathbf{u}}_1;\theta )$ around $\mathbf{ū}\left(t_1\right)$ as follows,

(17)$\begin{array}{c}\hfill M^{}(\mathbf{ū}\left(t_1\right)+\mathsf{\Delta }{\mathbf{u}}_1;\theta )\approx M^{}(\mathbf{ū}\left(t_1\right);\theta )+{\mathbf{D}}_M\left(\mathbf{ū}\left(t_1\right)\right)\mathsf{\Delta }{\mathbf{u}}_1,\end{array}$

where $\mathbf{u}(t_2=M^{}(\mathbf{u}\left(t_1\right);\theta )\approx M^{}(\mathbf{ū}\left(t_1\right)+\mathsf{\Delta }{\mathbf{u}}_1;\theta )$ and $\mathbf{u}\left(t_2\right)=M^{}(\mathbf{ū}\left(t_1\right);\theta )$. Consequently, $\mathsf{\Delta }{\mathbf{u}}_2=\mathbf{u}\left(t_2\right)-\mathbf{ū}\left(t_2\right)\approx {\mathbf{D}}_M\left(\mathbf{ū}\left(t_1\right)\right)\mathsf{\Delta }{\mathbf{u}}_1$, which can be generalized as,

(18)$\begin{array}{c}\hfill \mathsf{\Delta }{\mathbf{u}}_{k+1}\approx {\mathbf{D}}_M\left(\mathbf{ū}\left(t_k\right)\right)\mathsf{\Delta }{\mathbf{u}}_k,\end{array}$

with $\mathbf{u}\left(t_k\right)\approx \mathsf{\Delta }{\mathbf{u}}_k+\mathbf{ū}\left(t_k\right)$. It is customary to call Equation (18) as the perturbation equation, or the tangent linear system (TLS). Equation (18) can be related to $\mathsf{\Delta }{\mathbf{u}}_0$ by recursion as follows,

$\begin{array}{cc}\hfill \mathsf{\Delta }{\mathbf{u}}_{k+1}& \approx {\mathbf{D}}_M\left(\mathbf{ū}\left(t_k\right)\right)\mathsf{\Delta }{\mathbf{u}}_k\hfill \\ \hfill & \approx {\mathbf{D}}_M\left(\mathbf{ū}\left(t_k\right)\right){\mathbf{D}}_M\left(\mathbf{ū}\left(t_{k-1}\right)\right)\mathsf{\Delta }{\mathbf{u}}_{k-1}\hfill \\ \hfill & \approx {\mathbf{D}}_M\left(\mathbf{ū}\left(t_k\right)\right){\mathbf{D}}_M\left(\mathbf{ū}\left(t_{k-1}\right)\right){\mathbf{D}}_M\left(\mathbf{ū}\left(t_{k-2}\right)\right)\mathsf{\Delta }{\mathbf{u}}_{k-2}\hfill \\ \hfill & \approx {\mathbf{D}}_M\left(\mathbf{ū}\left(t_k\right)\right){\mathbf{D}}_M\left(\mathbf{ū}\left(t_{k-1}\right)\right){\mathbf{D}}_M\left(\mathbf{ū}\left(t_{k-2}\right)\right){\mathbf{D}}_M\left(\mathbf{ū}\left(t_{k-3}\right)\right)\dots {\mathbf{D}}_M\left(\mathbf{ū}\left(t_0\right)\right)\mathsf{\Delta }{\mathbf{u}}_0,\hfill \end{array}$

which can be short-handed as $\mathsf{\Delta }{\mathbf{u}}_{k+1}\approx {\mathbf{D}}_M\left(\mathbf{ū}\left(t_{k:0}\right)\right)\mathsf{\Delta }{\mathbf{u}}_0$ (please, notice the order of matrix multiplication and the subscript “$k:0$”).

Now, we investigate the first order variation $\mathsf{\Delta }J$ of the cost functional $J\left(\mathbf{ū}\left(t_0\right)\right)$ induced by the perturbation $\mathsf{\Delta }{\mathbf{u}}_0$ in the initial condition. This can be approximated as below,

(19)$\begin{array}{cc}\hfill \mathsf{\Delta }J& =\mathsf{\Delta }{\mathbf{u}}_0^T\mathsf{\nabla }J\left(\mathbf{ū}\left(t_0\right)\right)\hfill \end{array}$

(20)$\begin{array}{cc}\hfill & =-\sum _{t_k\in \mathcal{T}}\mathsf{\Delta }{\mathbf{u}}_k^T{\mathbf{D}}_h^T\left(\mathbf{ū}\left(t_k\right)\right){\mathbf{R}}^{-1}\left(t_k\right)\left(\mathbf{w}\left(t_k\right)-h\left(\mathbf{ū}\left(t_k\right)\right)\right).\hfill \end{array}$

Given that $\mathsf{\Delta }{\mathbf{u}}_k\approx {\mathbf{D}}_M\left(\mathbf{ū}\left(t_{k-1:0}\right)\right)\mathsf{\Delta }{\mathbf{u}}_0$, then $\mathsf{\Delta }{\mathbf{u}}_k^T\approx \mathsf{\Delta }{\mathbf{u}}_0^T{\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_0\right)\right){\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_1\right)\right)\dots {\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_{k-1}\right)\right)=\mathsf{\Delta }{\mathbf{u}}_0^T{\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_{0:k-1}\right)\right)$ and Equation (20) can be rewritten as

(21)$\begin{array}{cc}\hfill \mathsf{\Delta }J& =-\sum _{t_k\in \mathcal{T}}\mathsf{\Delta }{\mathbf{u}}_0^T{\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_{0:k-1}\right)\right){\mathbf{D}}_h^T\left(\mathbf{ū}\left(t_k\right)\right){\mathbf{R}}^{-1}\left(t_k\right)\left(\mathbf{w}\left(t_k\right)-h\left(\mathbf{ū}\left(t_k\right)\right)\right).\hfill \end{array}$

By comparing Equations (19) and (21), the gradient of the cost functional can be approximated as

(22)$\begin{array}{cc}\hfill \mathsf{\nabla }J\left(\mathbf{ū}\left(t_0\right)\right)& =-\sum _{t_k\in \mathcal{T}}{\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_{0:k-1}\right)\right){\mathbf{D}}_h^T\left(\mathbf{ū}\left(t_k\right)\right){\mathbf{R}}^{-1}\left(t_k\right)\left(\mathbf{w}\left(t_k\right)-h\left(\mathbf{ū}\left(t_k\right)\right)\right)\hfill \end{array}$

(23)$\begin{array}{cc}& =-\sum _{t_k\in \mathcal{T}}{\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_{0:k-1}\right)\right)\mathbf{f}\left(t_k\right),\hfill \end{array}$

where $\mathbf{f}\left(t_k\right)={\mathbf{D}}_h^T\left(\mathbf{ū}\left(t_k\right)\right){\mathbf{R}}^{-1}\left(t_k\right)\left(\mathbf{w}\left(t_k\right)-h\left(\mathbf{ū}\left(t_k\right)\right)\right)$. If we denote the time instants at which measurements are available as $\mathcal{T}=\{t_{O1},t_{O2},\dots ,t_{ON}\}$, Equation (23) can be expanded as

(24)$\begin{array}{cc}\hfill \mathsf{\nabla }J\left(\mathbf{ū}\left(t_0\right)\right)& =-\left\{{\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_{0:O1-1}\right)\right)\mathbf{f}\left(t_{O1}\right)+{\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_{0:O2-1}\right)\right)\mathbf{f}\left(t_{O2}\right)+\dots +{\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_{0:ON-1}\right)\right)\mathbf{f}\left(t_{ON}\right)\right\}.\hfill \end{array}$

Now, defining a sequence of ${\lambda }_k\in {\mathbb{R}}^n$ as below,

(25)${\lambda }_k=\left\{\begin{array}{cc}{\mathbf{f}}_k,\hfill & \mathrm{if}{t}_k=t_{ON}\hfill \\ {\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_k\right)\right){\lambda }_{k+1}+{\mathbf{f}}_k,\hfill & \mathrm{if}{t}_k\in \{t_{O1},t_{O2},\dots ,t_{ON-1}\}\hfill \\ {\mathbf{D}}_M^T\left(\mathbf{ū}\left(t_k\right)\right){\lambda }_{k+1},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$

It can be verified that $\mathsf{\nabla }J\left(\mathbf{ū}\left(t_0\right)\right)=-{\lambda }_0$ (assume some numbers and you can see this relation holds!). Therefore, in order to obtain the gradient of the cost functional, ${\lambda }_0$ has to be computed, which depends on the evaluation of ${\lambda }_1$. In turn, the computation of ${\lambda }_1$ requires ${\lambda }_2$ and so on. Equation (25) is known as the first-order adjoint equation, as it implies the evaluation of ${\lambda }_k$ sequence from $t_{k+1}$ to $t_k$ (i.e., reverse order).

Therefore, the first-order approximation of the 4DVAR works as follows. Starting from a prior guess of the initial condition, the base trajectory is computed by solving the model forward in time until the final time corresponding the last observation point (i.e., $k=0,1,2,\dots ,ON$). Then, the value of ${\lambda }_{ON}={\mathbf{f}}_{ON}$ is evaluated at final time. After that, ${\lambda }_k$ is evolved backward in time using Equation (25) until $\mathsf{\nabla }J\left(\mathbf{ū}\left(t_0\right)\right)=-{\lambda }_0$ is obtained. This value of the gradient is thus utilized to update the initial condition and a new base trajectory is generated. The solution of the 4DVAR problem requires the solution of the model dynamics forward in time and adjoint problem backward in time until to compute the gradient of the cost functional and update the initial condition. The process is thus repeated until convergence takes place. Listing 6 shows a sample function to compute the gradient of the cost functional $\mathsf{\nabla }J\left(\mathbf{ū}\left(t_0\right)\right)$ corresponding to a base trajectory generated from a guess of the initial condition $\mathbf{ū}\left(t_0\right)$. We highlight that, in practice, the storage of the base trajectory as well as the $\lambda$ sequence at every time instant might be overwhelming. However, we are not addressing such issues in these introductory tutorials to data assimilation techniques.

The gradient $\mathsf{\nabla }J\left(\mathbf{ū}\left(t_0\right)\right)$ should be used in a minimization algorithm to update the initial condition for the next iteration. One simple algorithm is the simple gradient descent where an updated value of the initial state is computed as $\mathbf{ū}\left(t_0\right){{)}^{new}=\mathbf{ū}\left(t_0\right))}^{old}-{\beta }_n\mathsf{\nabla }J\left(\mathbf{ū}{\left(t_0\right)}^{old}\right)$, where ${\beta }_n$ is some step parameter. This can be normalized as $\mathbf{ū}\left(t_0\right){{)}^{new}=\mathbf{ū}\left(t_0\right))}^{old}-\beta {\displaystyle \frac{\mathsf{\nabla }J\left(\mathbf{ū}{\left(t_0\right)}^{old}\right)}{\parallel \mathsf{\nabla }J\left(\mathbf{ū}{\left(t_0\right)}^{old}\right)\parallel }}$. The value of $\beta$ might be predefined, or more efficiently updated at each iteration using an additional optimization algorithm (e.g., line-search). For the sake of completeness, we present a line-search routine in Listing 7 using the Golden search algorithm. This is based on the definition of the cost functional in Listing 8.

Example: Lorenz 63 System

Similar to the 3DVAR demonstration, we apply the described 4DVAR using the first-order adjoint method on the Lorenz 63 system. We also begin with the same erroneous initial condition of $\mathbf{u}\left(0\right)={[2,3,4]}^T$ and observations are collected each $0.2$ time units, contaminated with a Gaussian noise with diagonal covariance matrix defined as $\mathbf{R}={\sigma }_m^2{\mathbf{I}}_3$, where ${\sigma }_m=0.15$ is the standard deviation for measurement noise. Moreover, we simply define a linear observation operator defined as $h\left(\mathbf{u}\right)=\mathbf{u}$, with a Jacobian of identity matrix. We utilize the simple gradient descent for minimizing the cost functional, equipped by a Golden search method for learning rate optimization. A maximum number of iterations is set to 1000, but we highlight that this is highly dependent on the adopted minimization algorithm as well as the line-search technique. In practice, the evaluation of each iteration might be too computationally expensive, so the number of iterations need to be as low as possible. We define two criteria for convergence, and iterations stop whenever any one of them is achieved. The first one is based on the change in the value of the cost or loss functional and the second one is based on the magnitude of its gradient. Extra criteria might be supplied as well.

Results for running Listing 9 is shown in Figure 2, where we can notice the significant improvement of predictions, compared to the background trajectories. Moreover, we highlight the correction to the initial conditions in Figure 2 which resulted in the analysis trajectory. This is opposed to the 3DVAR implementation, where correction is applied locally at measurements instants only as seen in Figure 1.

Before we move to other data assimilation techniques, we highlight a few remarks regarding our presentation of the 4DVAR

• In Listing 9, we utilize the gradient descent approach to minimize the cost function. Readers are encouraged to apply other optimization techniques (e.g., conjugate gradient) that achieve higher convergence rate.

• The determination of the learning rate can be further optimized using more efficient line-search methods, rather than the simple Golden search.

• The Lagrangian multiplier method can be applied to solve the 4DVAR problem instead of the adjoint method, similar results should be obtained.

• The presented algorithm relies on the definition of the cost functional given in Equation (13), based on the discrepancy between measurements and model’s predictions. When extra information is available, it can be incorporated into the cost functional. For instance, similar to Equation (4), a term that penalizes the correction magnitude can be added, weighted by the background covariance matrix. Furthermore, symmetries or other physical knowledge can be enforced as hard or weak constraints.

• The first-order adjoint algorithm requires the computation of the Jacobian ${\mathbf{D}}_M\left(\mathbf{u}\right)$ of the discrete-time model map $M(\mathbf{u};\theta )$. This can be computed by plugging the model $f(\mathbf{u};\theta )$ in a time integration scheme and rearranging everything to rewrite $M(\mathbf{u};\theta )$ as explicit function of $\mathbf{u}$ and differentiating with respect to components of $\mathbf{u}$. For Lorenz 63 and 1st Euler scheme, this can be an easy task. However, for a higher dimensional system and more accurate time integrators, this would be cumbersome. Instead, the chain rule can be utilized to compute ${\mathbf{D}}_M\left(\mathbf{u}\right)$ as presented in Listing 10, which takes as input the right-hand side of the continuous-time dynamics $f(·;·)$ (described in Listing 3 for Lorenz 63 system) as well as its Jacobian (given in Listing 11 for Lorenz 63).

5. Forward Sensitivity Method

We have seen in Section 4 that the minimization of the cost functional via the 4DVAR algorithm requires the solution of the adjoint problem at each iteration, which incurs a significant computational burden for high dimensional systems. Alternatively, Lakshmivarahan and Lewis [13] proposed the forward sensitivity method (FSM) to derive an expression for the correction vector in terms of the forward sensitivity matrices [14]. In their development, simultaneous correction to the initial condition $\mathbf{u}\left(t_0\right)$ and the model parameters $\theta$ is treated. For conciseness and consistency with the methods introduced here, we only consider erroneous initial conditions and assume model parameters are perfectly known. Given the discrete-time model map $M(\mathbf{u};\theta )$ in Equation (2), the forecast sensitivity at time $t_{k+a}$ to the initial conditions $\mathbf{u}\left(t_0\right)$ can be defined as follows,

(26)${\displaystyle \frac{\partial u_i\left(t_{k+1}\right)}{\partial u_j\left(t_0\right)}}=\sum _{q=1}^n\left({\displaystyle \frac{\partial M_i(\mathbf{u}\left(t_k\right);\theta )}{\partial u_q\left(t_k\right)}}\right)\left({\displaystyle \frac{\partial u_q\left(t_k\right)}{\partial u_j\left(t_0\right)}}\right),1⩽i,j⩽n,$

where $M(\mathbf{u}\left(t_k\right);\theta )={[M_1(\mathbf{u}\left(t_k\right);\theta ),M_2(\mathbf{u}\left(t_k\right);\theta ),\dots ,M_n(\mathbf{u}\left(t_k\right);\theta )]}^T$. Recall that the Jacobian of the model $M(\mathbf{u}\left(t_k\right);\theta )$ is defined by the matrix ${\mathbf{D}}_M\left(\mathbf{u}\left(t_k\right)\right)\in {\mathbb{R}}^{n\times n}$ whose ${(i,j)}^{th}$ entry is defined as ${\displaystyle \frac{\partial M_i(\mathbf{u}\left(t_k\right);\theta )}{\partial u_j\left(t_k\right)}}$. We also define $\mathbf{U}\left(t_k\right)$ as the forward sensitivity matrix of $\mathbf{u}\left(t_k\right)\in {\mathbb{R}}^{n\times n}$ with respect to initial state $\mathbf{u}\left(t_0\right)$, where ${\left[\mathbf{U}\left(t_k\right)\right]}_{i,j}={\displaystyle \frac{\partial u_i\left(t_k\right)}{\partial u_j\left(t_0\right)}}$ for $1⩽i,j⩽n$. Thus, Equation (26) can be rewritten in matrix form as,

(27)$\begin{array}{cc}\hfill \mathbf{U}\left(t_{k+1}\right)& ={\mathbf{D}}_M\left(\mathbf{u}\left(t_k\right)\right)\mathbf{U}\left(t_k\right).\hfill \end{array}$

Equation (27) provides the dynamic evolution of the forward sensitivity matrix in a recursive manner, initialized by $\mathbf{U}\left(t_0\right)={\mathbf{I}}_n$, that can be used to relate the prediction error at any time step to the initial condition.

Given the measurement $\mathbf{w}\left(t_k\right)$ at time $t_k\in \mathcal{T}$, the forecast error $\mathbf{e}\left(t_k\right)$ is defined as the difference between the model forecast and measurements as

(28)$\mathbf{e}\left(t_k\right)=\mathbf{w}\left(t_k\right)-\mathbf{h}\left(\mathbf{u}\left(t_k\right)\right).$

This is commonly called the innovation in DA terminology. The cost functional in Equation (13) can be rewritten as

(29)$J\left(\mathbf{u}\left(t_0\right)\right)=\sum _{t_k\in \mathcal{T}}{\displaystyle \frac{1}{2}}{\parallel \mathbf{e}\left(t_k\right)\parallel }_{{\mathbf{R}}^{-1}\left(t_k\right)}^2=\sum _{t_k\in \mathcal{T}}{\displaystyle \frac{1}{2}}\mathbf{e}{\left(t_k\right)}^T{\mathbf{R}}^{-1}\left(t_k\right)\mathbf{e}\left(t_k\right).$

With the assumption that the dynamical model is perfect (i.e., correctly encapsulates all the relevant processes) and the model parameters are known, the deterministic part of the forecast error can be attributed to the inaccuracy in the initial condition $\mathbf{u}\left(t_0\right)$, defined as $\mathsf{\Delta }{\mathbf{u}}_0={\mathbf{u}}_t\left(t_0\right)-{\mathbf{u}}_b\left(t_0\right)$, where ${\mathbf{u}}_t\left(t_0\right)$ denotes the true initial conditions.

Considering a base trajectory $\mathbf{ū}\left(t_k\right)$ for $k=1,2,\dots$ generated from the initial condition of $\mathbf{ū}\left(t_0\right)$, related to the corrected trajectory ($\mathbf{u}\left(t_k\right)$ for $k=1,2,\dots ,$) obtained by correcting the initial condition as $\mathbf{u}\left(t_0\right)=\mathbf{ū}\left(t_0\right)+\mathsf{\Delta }{\mathbf{u}}_0$, we define the difference between both trajectories at any time $t_k$ as $\mathsf{\Delta }{\mathbf{u}}_k=\mathbf{u}\left(t_k\right)-\mathbf{ū}\left(t_k\right)$. We highlight that $\mathbf{u}\left(t_k\right)$ is a function of both the initial condition (with the model parameters being known), the first-order Taylor expansion of $\mathbf{u}\left(t_k\right)$ around the base trajectory can be written as $\mathbf{u}\left(t_k\right)\approx \mathbf{ū}\left(t_k\right)+\mathbf{U}\left(t_k\right)\mathsf{\Delta }{\mathbf{u}}_0$, leading to the following relation

(30)$\mathsf{\Delta }{\mathbf{u}}_k\approx \mathbf{U}\left(t_k\right)\mathsf{\Delta }{\mathbf{u}}_0.$

If we let the perturbed (corrected) trajectory to be the sought true trajectory, Equation (3) can be rewritten as

(31)$\begin{array}{c}\hfill \mathbf{w}\left(t_k\right)=h(\mathbf{ū}\left(t_k\right)+\mathsf{\Delta }{\mathbf{u}}_k)+{\xi }_m,\end{array}$

and a first order expansion of $\mathbf{w}\left(t_k\right)$ (neglecting the measurement noise) will be as follows,

(32)$\begin{array}{c}\hfill \mathbf{w}\left(t_k\right)\approx h\left(\mathbf{ū}\left(t_k\right)\right)+{\mathbf{D}}_h\left(\mathbf{ū}\left(t_k\right)\right)\mathsf{\Delta }{\mathbf{u}}_k,\end{array}$

and the forecast error at the base trajectory can be approximated as

(33)$\mathbf{e}\left(t_k\right)={\mathbf{D}}_h\left(\mathbf{ū}\left(t_k\right)\right)\mathsf{\Delta }{\mathbf{u}}_k.$

Equations (30) and (33) can be combined to yield the following,

(34)$\mathbf{e}\left(t_k\right)={\mathbf{D}}_h\left(\mathbf{ū}\left(t_k\right)\right)\mathbf{U}\left(t_k\right)\mathsf{\Delta }{\mathbf{u}}_0,$

which relates the forecast error at any time $t_k$ and the discrepancy between the true and erroneous initial condition in a linear relationship. In order to account for all the time instants at which observations are available ($\mathcal{T}=\{t_{O1},t_{O2},\dots ,t_{ON}\}$), Equation (34) can be concatenated at different times and written as a linear system of equations as follows,

(35)$\begin{array}{c}\hfill \mathbf{Q}\mathsf{\Delta }{\mathbf{u}}_0={\mathbf{e}}_F,\end{array}$

where the matrix ${\mathbf{Q}}^{Nm\times n}$ and the vector ${\mathbf{e}}_F\in {\mathbb{R}}^{Nm}$ are computed as,

(36)$\begin{array}{c}\hfill \mathbf{Q}=\left[\begin{array}{c}{\mathbf{D}}_h\left(\mathbf{ū}\left(t_{O1}\right)\right)\mathbf{U}\left(t_{O1}\right)\\ {\mathbf{D}}_h\left(\mathbf{ū}\left(t_{O2}\right)\right)\mathbf{U}{t}_{O2})\\ ⋮\\ {\mathbf{D}}_h\left(\mathbf{ū}\left(t_{ON}\right)\right)\mathbf{U}\left(t_{ON}\right)\end{array}\right],{\mathbf{e}}_F=\left[\begin{array}{c}\mathbf{e}\left(t_{O1}\right)\\ \mathbf{e}\left(t_{O2}\right)\\ ⋮\\ \mathbf{e}\left(t_{ON}\right)\end{array}\right].\end{array}$

Depending on the value of $Nm$ relative to n, Equation (35) can give rise to either an over-determined or an under-determined linear inverse problem. In either case, the inverse problem can be solved in a weighted least squares sense to find a first-order estimation of the optimal correction or perturbation to the initial condition $\mathsf{\Delta }{\mathbf{u}}_0$, with ${\mathbf{R}}^{-1}$ being the weighting matrix, where $\mathbf{R}$ is a block-diagonal matrix constructed as follows,

(37)$\mathbf{R}=\left[\begin{array}{cccc}\mathbf{R}\left(t_{O1}\right)& & & \\ & \mathbf{R}\left(t_{O2}\right)& & \\ & & \ddots & \\ & & & \mathbf{R}\left(t_{ON}\right)\end{array}\right],$

and the solution of Equation (35) can be written as $\mathsf{\Delta }{\mathbf{u}}_0={\left({\mathbf{Q}}^T{\mathbf{R}}^{-1}\mathbf{Q}\right)}^{-1}{\mathbf{Q}}^T{\mathbf{R}}^{-1}{\mathbf{e}}_F$ for the over-determined case. This first order approximation progressively yields better results by repeating the entire process for multiple iterations until convergence with certain tolerance [13]. In essence, the FSM is an alternative to the 4DVAR algorithm, replacing the solution of the adjoint problem backward in time (i.e., Equation (25)), by the successive matrix evaluation in Equation (27). However, we highlight that in the 4DVAR approach, the actual forecast error is computed as ${\mathbf{e}}_k=\mathbf{w}\left(t_k\right)-h\left(\mathbf{ū}\left(t_k\right)\right)$ while its first-order approximation is utilized in the FSM development. The duality between the two approaches is further discussed in [13].

A Python implementation of the FSM approach is presented in Listing 12. We note that we solve the Equation (35) using the built-in numpy least-squares function. However, more efficient iterative schemes can be adopted in practice.

Example: Lorenz 63 System

We apply the described FSM to estimate the initial conditions for the Lorenz 63 system, using the same parameters and setup as described in Section 4. Sample code is presented in Listing 13. We highlight that instead of adding the correction vector $\mathsf{\Delta }{\mathbf{u}}_0$ directly to the base value $\mathbf{ū}\left(t_0\right)$, we multiply it with a learning rate to mitigate the effects of first-order approximations. We utilize the golden search method to update this learning rate at each iteration.

Prediction results are provided in Figure 3, where we notice a large discrepancy at the estimated initial conditions. However, the predicted trajectory perfectly match the true one for the rest of the testing time window. This is largely affected by the nature of the Lorenz system itself and the attachment to its attractor. Furthermore, this can be partially attributed to the lack of background information and its contribution to the cost functional. Moreover, this can be highly improved by adding more observations close to the initial time since the correction vector is estimated based on the forecast error computed at observation times. Anyhow, we see that the analysis trajectory is significantly more accurate than the background one, with iterative first-order approximations of the forward sensitivity method, even for long time predictions.

6. Kalman Filtering

The idea behind Kalman filtering techniques is to propagate the mean as well as the covariance matrix of the system’s state sequentially in time. That is, in addition to providing an improved state estimate (i.e., the analysis), it also gives some information about the statistical properties of this state estimate. This is one main difference between Kalman filtering and variational methods, which often assumes a fixed (stationary) background covariance matrices. Kalman filters are also very popular in systems engineering, robotics, navigation, and control. Almost all modern control systems use the Kalman filter. It assisted the guidance of the Apollo 11 lunar module to the moon’s surface, and most probably will do the same for next generations of aircraft as well.

Although most application in fluid dynamics involve nonlinear systems, we first describe the standard Kalman filter developed for the linear dynamical system case with linear observation operator described as

(38)$\begin{array}{cc}\hfill {\mathbf{u}}_t\left(t_{k+1}\right)& ={\mathbf{M}}_k{\mathbf{u}}_t\left(t_k\right)+{\xi }_p\left(t_{k+1}\right),\hfill \end{array}$

(39)$\begin{array}{cc}\hfill \mathbf{w}\left(t_k\right)& ={\mathbf{H}}_k{\mathbf{u}}_t\left(t_k\right)+{\xi }_m\left(t_k\right)\hfill \end{array}$

where $\mathbf{M}\in {\mathbb{R}}^{n\times n}$ is a non-singular system matrix defining the underlying governing processes and ${\xi }_p\in {\mathbb{R}}^n$ describes the process noise (or model error). $\mathbf{H}\in {\mathbb{R}}^{m\times n}$ represents the measurement system with a measurement noise of ${\xi }_m\in {\mathbb{R}}^m$.

As presented in Section 2, the true state ${\mathbf{u}}_t\left(t_k\right)$ is assumed to be a random variable with known mean $E\left[{\mathbf{u}}_t\left(t_k\right)\right]={\mathbf{u}}_b\left(t_k\right)$ and covariance matrix of $E\left[({\mathbf{u}}_t\left(t_k\right)-{\mathbf{u}}_b\left(t_k\right)){({\mathbf{u}}_t\left(t_k\right)-{\mathbf{u}}_b\left(t_k\right))}^T\right]={\mathbf{B}}_k$. In Kalman filtering, we note that the covariance matrix evolves in time, and thus appears the subscript. We also assume that the process noise is unbiased with zero mean and a covariance matrix $\mathbf{Q}$. That is $E\left({\xi }_p\left(t_k\right)\right)=0$ and $E\left({\xi }_p\left(t_k\right){\xi }_p{\left(t_k\right)}^T\right)={\mathbf{Q}}_k$.

Thus, the goal of the filtering problem is to find a good estimate (analysis) ${\mathbf{u}}_a\left(t_k\right)$ of the true system’s state ${\mathbf{u}}_t\left(t_k\right)$ given a dynamical model a set of noisy observation $\left\{\mathbf{w}\left(t_i\right)\right\}$ collected at some time instants $t_i\in (0,t_k]$. The optimality of the estimate ${\mathbf{u}}_a\left(t_k\right)$ is defined as the one which minimizes $E\left[{({\mathbf{u}}_t\left(t_k\right)-{\mathbf{u}}_a\left(t_k\right))}^T({\mathbf{u}}_t\left(t_k\right)-{\mathbf{u}}_a\left(t_k\right))\right]$. This filtering process generally consists of two steps: the forecast step and the data assimilation step.

The forecast step is performed using the predictable part of the given dynamical model starting from the best known information at time $t_k$ (denoted as ${\widehat{\mathbf{u}}}_b\left(t_k\right))$ to produce a forecast or background estimate ${\mathbf{u}}_b\left(t_{k+1}\right)={\mathbf{M}}_k{\widehat{\mathbf{u}}}_b\left(t_k\right)$. The difference between the background forecast and true state at $t_{k+1}$ can be written as follows,

$\begin{array}{cc}\hfill {\xi }_b\left(t_{k+1}\right)& ={\mathbf{u}}_t\left(t_{k+1}\right)-{\mathbf{u}}_b\left(t_{k+1}\right)\hfill \\ \hfill & =({\mathbf{M}}_k{\mathbf{u}}_t\left(t_k\right)+{\xi }_p\left(t_{k+1}\right))-{\mathbf{M}}_k{\widehat{\mathbf{u}}}_b\left(t_k\right)\hfill \\ \hfill & ={\mathbf{M}}_k({\mathbf{u}}_t\left(t_k\right)-{\widehat{\mathbf{u}}}_b\left(t_k\right))+{\xi }_p\left(t_{k+1}\right)\hfill \\ \hfill & ={\mathbf{M}}_k{\widehat{\xi }}_b\left(t_k\right)+{\xi }_p\left(t_{k+1}\right),\hfill \end{array}$

where ${\widehat{\xi }}_b\left(t_k\right)=({\mathbf{u}}_t\left(t_k\right)-{\widehat{\mathbf{u}}}_b\left(t_k\right))$ is the error estimate at $t_k$, with zero mean and covariance matrix of ${\widehat{\mathbf{B}}}_k$.

The covariance matrix of the background estimate at $t_{k+1}$ can be evaluated as ${\mathbf{B}}_{k+1}=E\left[{\xi }_b\left(t_{k+1}\right){\xi }_b{\left(t_{k+1}\right)}^T\right]=E\left[\left({\mathbf{M}}_k{\widehat{\xi }}_b\left(t_k\right)+{\xi }_p\left(t_{k+1}\right)\right){\left({\mathbf{M}}_k{\widehat{\xi }}_b\left(t_k\right)+{\xi }_p\left(t_{k+1}\right)\right)}^T\right]$. Since, ${\widehat{\xi }}_b\left(t_k\right)$ and ${\xi }_p\left(t_{k+1}\right)$ are assumed to be uncorrelated (i.e., $E\left[{\widehat{\xi }}_b\left(t_k\right){\xi }_p{\left(t_{k+1}\right)}^T\right]=\mathbf{0}$), the background covariance matrix at $t_{k+1}$ can be computed as follows,

(40)${\mathbf{B}}_{k+1}={\mathbf{M}}_k{\widehat{\mathbf{B}}}_k{\mathbf{M}}_k^T+{\mathbf{Q}}_{k+1}.$

Now, with the forecast step, we have a background estimate at $t_{k+1}$ defined as ${\mathbf{u}}_b\left(t_{k+1}\right)$ with a covariance matrix ${\mathbf{B}}_{k+1}$. Then, measurements $\mathbf{w}\left(t_{k+1}\right)$ are collected at $t_{k+1}$ with a linear operator ${\mathbf{H}}_{k+1}$ and measurement noise ${\xi }_m\left(t_{k+1}\right)$ with zero mean a covariance matrix of ${\mathbf{R}}_{k+1}$. Thus, we would like to fuse these pieces of information to create an optimal unbiased estimate (analysis) ${\mathbf{u}}_a\left(t_{k+1}\right)$ with a covariance matrix ${\mathbf{P}}_{k+1}$. This can be defines as linear function of ${\mathbf{u}}_b\left(t_{k+1}\right)$ and $\mathbf{w}\left(t_{k+1}\right)$ as follows,

(41)$\begin{array}{c}\hfill {\mathbf{u}}_a\left(t_{k+1}\right)={\mathbf{u}}_b\left(t_{k+1}\right)+{\mathbf{K}}_{k+1}\left(\mathbf{w}\left(t_{k+1}\right)-{\mathbf{H}}_{k+1}{\mathbf{u}}_b\left(t_{k+1}\right)\right),\end{array}$

where $\left(\mathbf{w}\left(t_{k+1}\right)-{\mathbf{H}}_{k+1}{\mathbf{u}}_b(t_{k+1})\right)$ is the innovation vector and $\mathbf{K}\in {\mathbb{R}}^{n\times m}$ is called the Kalman gain matrix. We highlight that Kalman gain matrix is defined in such a way to minimize $E\left[{({\mathbf{u}}_t\left(t_{k+1}\right)-{\mathbf{u}}_a\left(t_{k+1}\right))}^T({\mathbf{u}}_t\left(t_{k+1}\right)-{\mathbf{u}}_a\left(t_{k+1}\right))\right]=\mathrm{tr}\left({\mathbf{P}}_{k+1}\right)$. This can be written as [27]

(42)$\begin{array}{c}\hfill {\mathbf{K}}_{k+1}={\mathbf{B}}_{k+1}{\mathbf{H}}_{k+1}^T{({\mathbf{H}}_{k+1}{\mathbf{B}}_{k+1}{\mathbf{H}}_{k+1}^T+{\mathbf{R}}_{k+1})}^{-1},\end{array}$

resulting in an analysis covariance matrix defined as

(43)$\begin{array}{c}\hfill {\mathbf{P}}_{k+1}=({\mathbf{I}}_n-{\mathbf{K}}_{k+1}{\mathbf{H}}_{k+1}){\mathbf{B}}_{k+1},\end{array}$

where ${\mathbf{I}}_n$ is the $n\times n$ identity matrix. The resulting analysis ${\mathbf{u}}_a\left(t_{k+1}\right)$ is known as the best linear unbiased estimate (BLUE). We highlight that information at $t_k$ might correspond to the analysis (i.e., ${\widehat{\mathbf{u}}}_b\left(t_k\right)={\mathbf{u}}_a\left(t_k\right)$) obtained from the last data assimilation implementation, or just from previous forecast if no other information is available. Thus, the Kalman filtering process can be summarized as follows,

$\begin{array}{cc}\hfill \mathrm{Inputs}:& {\widehat{\mathbf{u}}}_b\left(t_k\right),{\widehat{\mathbf{B}}}_k,{\mathbf{M}}_k,{\mathbf{Q}}_{k+1},\mathbf{w}\left(t_{k+1}\right),{\mathbf{R}}_{k+1},{\mathbf{H}}_{k+1}\hfill \\ \hfill \mathrm{Forecast}:& {\mathbf{u}}_b\left(t_{k+1}\right)={\mathbf{M}}_k{\widehat{\mathbf{u}}}_b\left(t_k\right)\hfill \\ \hfill & {\mathbf{B}}_{k+1}={\mathbf{M}}_k{\widehat{\mathbf{B}}}_k{\mathbf{M}}_k^T+{\mathbf{Q}}_{k+1}\hfill \\ \hfill \mathrm{Kalman}\mathrm{gain}:& {\mathbf{K}}_{k+1}={\mathbf{B}}_{k+1}{\mathbf{H}}_{k+1}^T{\left({\mathbf{H}}_{k+1}{\mathbf{B}}_{k+1}{\mathbf{H}}_{k+1}^T+{\mathbf{R}}_{k+1}\right)}^{-1}\hfill \\ \hfill \mathrm{Analysis}:& {\mathbf{u}}_a\left(t_{k+1}\right)={\mathbf{u}}_b\left(t_{k+1}\right)+{\mathbf{K}}_{k+1}\left(\mathbf{w}\left(t_{k+1}\right)-{\mathbf{H}}_{k+1}{\mathbf{u}}_b\left(t_{k+1}\right)\right)\hfill \\ \hfill & {\mathbf{P}}_{k+1}=({\mathbf{I}}_n-{\mathbf{K}}_{k+1}{\mathbf{H}}_{k+1}){\mathbf{B}}_{k+1},\hfill \end{array}$

where the inputs at $t_k$ are defined as

$\begin{array}{c}\hfill ({\widehat{\mathbf{u}}}_b\left(t_k\right),{\widehat{\mathbf{B}}}_k)=\left\{\begin{array}{cc}({\mathbf{u}}_a\left(t_k\right),{\mathbf{P}}_k)\hfill & \mathrm{if}\mathbf{w}\left(t_k\right),\mathrm{is}\mathrm{available},\hfill \\ ({\mathbf{u}}_b\left(t_k\right),{\mathbf{B}}_k)\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$

Listing 14 describes a basic Python implementation of the data assimilation step using the KF algorithm described before. Although efficient matrix inversion routines that benefit from specific matrix properties can be utilized, we use the standard built-in Numpy matrix inversion function.

Different forms for evaluating the Kalman gain and the covariance matrices are presented in literature. Some of them are favored for computational cost aspects, while others maintain desirable properties (e.g., symmetry and positive definiteness) for numerically stable implementation [27]. Since we are more interested in nonlinear dynamical models, we shall discuss extensions for standard Kalman filters to account for nonlinearity in the following sections.

7. Extended Kalman Filter

Instead of dealing with linear stochastic dynamics, we look at the nonlinear case with general (nonlinear) observation operator written as

(44)$\begin{array}{cc}\hfill {\mathbf{u}}_t\left(t_{k+1}\right)& =M({\mathbf{u}}_t\left(t_k\right);\theta )+{\xi }_p\left(t_{k+1}\right),\hfill \end{array}$

(45)$\begin{array}{cc}\hfill \mathbf{w}\left(t_k\right)& =h\left({\mathbf{u}}_t\left(t_k\right)\right)+{\xi }_m\left(t_k\right).\hfill \end{array}$

The first challenge of applying Kalman filter for this system is the propagation of the background covariance matrix in the forecast step. The main clue behind the extended Kalman filter (EKF) to address this issue is to locally linearize $M\left(\mathbf{u}\left(t_k\right)\right)$ by expanding it around the estimate ${\widehat{\mathbf{u}}}_b\left(t_k\right)$ at $t_k$ using the first-order Taylor series as follows,

(46)$M({\mathbf{u}}_t\left(t_k\right);\theta )\approx M({\widehat{\mathbf{u}}}_b\left(t_k\right);\theta )+{\mathbf{D}}_M\left({\widehat{\mathbf{u}}}_b\left(t_k\right)\right){\widehat{\xi }}_b\left(t_k\right),$

where ${\mathbf{D}}_M\left({\widehat{\mathbf{u}}}_b\left(t_k\right)\right)$ is the Jacobian (also known as the tangent linear operator) of the forward model $M(·;·)$ evaluated at ${\widehat{\mathbf{u}}}_b\left(t_k\right)$ and ${\widehat{\xi }}_b\left(t_k\right)=({\mathbf{u}}_t\left(t_k\right)-{\widehat{\mathbf{u}}}_b\left(t_k\right))$ defining the error estimate at $t_k$, with zero mean and covariance matrix of ${\widehat{\mathbf{B}}}_k$. Thus, the difference between the background forecast and true state at $t_{k+1}$ can be written as follows,

$\begin{array}{cc}\hfill {\xi }_b\left(t_{k+1}\right)& ={\mathbf{u}}_t\left(t_{k+1}\right)-{\mathbf{u}}_b\left(t_{k+1}\right)\hfill \\ \hfill & =M({\mathbf{u}}_t\left(t_k\right);\theta )+{\xi }_p\left(t_{k+1}\right)-M({\widehat{\mathbf{u}}}_b\left(t_k\right);\theta )\hfill \\ \hfill & \approx M({\widehat{\mathbf{u}}}_b\left(t_k\right);\theta )+{\mathbf{D}}_M\left({\widehat{\mathbf{u}}}_b\left(t_k\right)\right){\widehat{\xi }}_b\left(t_k\right)+{\xi }_p\left(t_{k+1}\right)-M({\widehat{\mathbf{u}}}_b\left(t_k\right);\theta )\hfill \\ \hfill & \approx {\mathbf{D}}_M\left({\widehat{\mathbf{u}}}_b\left(t_k\right)\right){\widehat{\xi }}_b\left(t_k\right)+{\xi }_p\left(t_{k+1}\right).\hfill \end{array}$

Similar to the derivation in the linear case, with the assumption of uncorrelation between ${\widehat{\xi }}_b\left(t_k\right)$ and ${\xi }_p\left(t_{k+1}\right)$, the background covariance matrix at $t_{k+1}$ can be computed as follows,

(47)${\mathbf{B}}_{k+1}={\mathbf{D}}_M\left({\widehat{\mathbf{u}}}_b\left(t_k\right)\right){\widehat{\mathbf{B}}}_k{\mathbf{M}}_k^T+{\mathbf{Q}}_{k+1}.$

The next challenge regarding the analysis step is the computation of the Kalman gain in case of nonlinear observation operator. Again, $h\left({\mathbf{u}}_t\left(t_{k+1}\right)\right)$ is linearized using Talylor series expansion around ${\mathbf{u}}_b\left(t_{k+1}\right)$ (i.e., the background forecast) as follows,

(48)$h\left({\mathbf{u}}_t\left(t_{k+1}\right)\right)\approx h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)+{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right){\xi }_b\left(t_{k+1}\right),$

where ${\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)$ is the Jacobian of the observation operator h, computed with the forecast ${\mathbf{u}}_b\left(t_{k+1}\right)$. The Kalman gain is thus computed using this first-order approximation of h as follows,

(49)$\begin{array}{c}\hfill {\mathbf{K}}_{k+1}={\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T{\left({\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right){\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T+{\mathbf{R}}_{k+1}\right)}^{-1},\end{array}$

with an analysis estimate and analysis covariance matrix defined as

(50)$\begin{array}{cc}\hfill {\mathbf{u}}_a\left(t_{k+1}\right)& ={\mathbf{u}}_b\left(t_{k+1}\right)+{\mathbf{K}}_{k+1}\left(\mathbf{w}\left(t_{k+1}\right)-h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)\right),\hfill \end{array}$

(51)$\begin{array}{cc}\hfill {\mathbf{P}}_{k+1}& =\left({\mathbf{I}}_n-{\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)\right){\mathbf{B}}_{k+1}.\hfill \end{array}$

A summary of the EKF algorithm is described as follows,

$\begin{array}{cc}\hfill \mathrm{Inputs}:& {\widehat{\mathbf{u}}}_b\left(t_k\right),{\widehat{\mathbf{B}}}_k,M(·;·),{\mathbf{Q}}_{k+1},\mathbf{w}\left(t_{k+1}\right),{\mathbf{R}}_{k+1},h(·)\hfill \\ \hfill \mathrm{Forecast}:& {\mathbf{u}}_b\left(t_{k+1}\right)=M({\widehat{\mathbf{u}}}_b\left(t_k\right);\theta )\hfill \\ \hfill & {\mathbf{B}}_{k+1}={\mathbf{D}}_M\left({\widehat{\mathbf{u}}}_b\left(t_k\right)\right){\widehat{\mathbf{B}}}_k{\mathbf{D}}_M{\left({\widehat{\mathbf{u}}}_b\left(t_k\right)\right)}^T+{\mathbf{Q}}_{k+1}\hfill \\ \hfill \mathrm{Kalman}\mathrm{gain}:& {\mathbf{K}}_{k+1}={\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T{\left({\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}_{k+1}{\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T+{\mathbf{R}}_{k+1}\right)}^{-1}\hfill \\ \hfill \mathrm{Analysis}:& {\mathbf{u}}_a\left(t_{k+1}\right)={\mathbf{u}}_b\left(t_{k+1}\right)+{\mathbf{K}}_{k+1}\left(\mathbf{w}\left(t_{k+1}\right)-h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)\right)\hfill \\ \hfill & {\mathbf{P}}_{k+1}=({\mathbf{I}}_n-{\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)){\mathbf{B}}_{k+1},\hfill \end{array}$

and a Python implementation of the data assimilation step is presented in Listing 15.

Example: Lorenz 63 System

The first-order approximation of the Kalman filter in nonlinear case, known as extended Kalman filter, is applied for the test case of Lorenz 63 system. The computation of model Jacobian ${\mathbf{D}}_M(·)$ is presented in Listing 10 and Listing 11 in Section 4. We use the same parameters and initial conditions for the twin experiment framework as before. The sequential implementation of the forecast and analysis steps is shown in Listing 16 and results are illustrated in Figure 4. We adopt the 4th order Runge–Kutta scheme for time integration. For demonstration purposes, we consider zero process noise (i.e., $\mathbf{Q}=0$). However, we have found that assuming non-zero process noise (e.g., $\mathbf{Q}=0.01{\mathbf{I}}_3$) yields better performance.

8. Ensemble Kalman Filter

Despite the sound mathematical and theoretical foundation of Kalman filters (both linear and nonlinear cases), they are not widely utilized in geophysical sciences. The major bottleneck in the computational pipeline of Kalman filtering is the update of background covariance matrix. In typical implementation, the cost of this step is $O\left(n^3\right)$, where n is the size of the state vector. For systems governed by ordinary differential equations (ODES), n can be manageable (e.g., 3 in the Lorenz 63 model). However, fluid flows are often governed by partial differential equations. Thus, spatial discretization schemes (e.g., finite difference, finite volume, and finite element) are applied, resulting in a semi-discrete system of ODEs. In geophysical flow dynamics applications (e.g., weather forecast), a dimension of millions or even billions is not uncommon, which hinders the feasible implementation of standard Kalman filtering techniques.

Alternatively, reduced rank algorithms that provides low-order approximation of the covariance matrices are usually adopted. A very popular approach is the ensemble Kalman filter (EnKF), introduced by Evensen [17,28,29] based on the Monte Carlo estimation methods.

The main procedure for these methods is to create an ensemble of size N of the system state denoted as $\left\{\mathbf{u}{\left(t_k\right)}^{\left(i\right)}\right|1\le i\le N\}$ and apply the filtering algorithm to each member of the established ensemble. The statistical properties of the forecast and analysis are thus extracted from the ensemble using the standard Monte Carlo framework. In the previous discussions, the forecast (background) and analysis covariances are defined as follows,

$\begin{array}{cc}\hfill \mathbf{B}& =E\left[{\xi }_b{\xi }_b^T\right]=E\left[({\mathbf{u}}_t-{\mathbf{u}}_b){({\mathbf{u}}_t-{\mathbf{u}}_b)}^T\right],\hfill \\ \hfill \mathbf{P}& =E\left[{\xi }_a{\xi }_a^T\right]=E\left[({\mathbf{u}}_t-{\mathbf{u}}_a){({\mathbf{u}}_t-{\mathbf{u}}_a)}^T\right].\hfill \end{array}$

Alternatively, those can be approximated by the ensemble covariances, given as

$\begin{array}{cc}\hfill \mathbf{B}& \approx {\displaystyle \frac{1}{N-1}}\sum _{i=1}^N({\mathbf{u}}_b^{\left(i\right)}-{\mathbf{ū}}_b){({\mathbf{u}}_b^{\left(i\right)}-{\mathbf{ū}}_b)}^T,\hfill \\ \hfill \mathbf{P}& \approx {\displaystyle \frac{1}{N-1}}\sum _{i=1}^N({\mathbf{u}}_a^{\left(i\right)}-{\mathbf{ū}}_a){({\mathbf{u}}_a^{\left(i\right)}-{\mathbf{ū}}_a)}^T,\hfill \end{array}$

where the bar denotes the ensemble average defined as $\mathbf{ū}={\displaystyle \frac{1}{N}}{\sum }_{i=1}^N{\mathbf{u}}^{\left(i\right)}$. Thus, an interpretation of EnKF is that the ensemble mean is the best estimate of the state and the spreading of the ensemble around the mean is a definition of the error in this estimate. A larger ensemble size N yields a better approximation of the state estimate and its covariance matrix. In the following, we describe the typical steps for applying the EnKF algorithm.

We begin by creating an initial ensemble $\{{\widehat{\mathbf{u}}}_b^{\left(i\right)}\left(t_k\right)|1\le i\le N\}$ at time $t_k$ drawn from the distribution $\mathcal{N}({\widehat{\mathbf{u}}}_b\left(t_k\right),{\widehat{\mathbf{B}}}_k)$, where ${\widehat{\mathbf{u}}}_b\left(t_k\right)$ represents our best-known estimate at $t_k$. It can be verified that the ensemble mean and covariance converge to ${\widehat{\mathbf{u}}}_b\left(t_k\right)$ and ${\mathbf{B}}_k$ as $N\to \infty$. Then, the forecast step is applied to each member of the enesemble as follows,

(52)${\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)=M({\widehat{\mathbf{u}}}_b^{\left(i\right)}\left(t_k\right);\theta )+{\xi }_p^{\left(i\right)}\left(t_{k+1}\right),$

where ${\xi }_p^{\left(i\right)}\left(t_{k+1}\right)$ is drawn from the multivariate Gaussian distribution with zero mean and covariance matrix of ${\mathbf{Q}}_{k+1}$ representing the process noise applied to each member. The sample mean of the forecast ensemble can be thus computed, along with the corresponding covaiance matrix as,

(53)$\begin{array}{cc}\hfill {\mathbf{u}}_b\left(t_{k+1}\right)& \approx {\mathbf{ū}}_b\left(t_{k+1}\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right),\hfill \end{array}$

(54)$\begin{array}{cc}\hfill {\xi }_b^{\left(i\right)}\left(t_{k+1}\right)& ={\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)-{\mathbf{ū}}_b\left(t_{k+1}\right),\hfill \end{array}$

(55)$\begin{array}{cc}\hfill {\mathbf{B}}_{k+1}& \approx {\displaystyle \frac{1}{N-1}}\sum _{i=1}^N{\xi }_b^{\left(i\right)}\left(t_{k+1}\right){\xi }_b^{\left(i\right)}{\left(t_{k+1}\right)}^T,\hfill \end{array}$

which provides an approximation for the background covariance at $t_{k+1}$ without actually propagating the covariance matrix, as is the case in standard Kalman filtering.

An enesmble of observations $\left\{{\mathbf{w}}^{\left(i\right)}\left(t_{k+1}\right)\right|1\le i\le N\}$, also called virtual observations, is created assuming a Gaussian distribution with a mean equal to the actual observation $\mathbf{w}\left(t_{k+1}\right)$ and a covariance matrix ${\mathbf{R}}_{k+1}$. In other words, random Gaussian perturbations with zero mean and covariance matrix ${\mathbf{R}}_{k+1}$ are added to the actual measurements to create perturbed measurements. The Kalman gain matrix can be computed as before (repeated here for completeness),

(56)$\begin{array}{c}\hfill {\mathbf{K}}_{k+1}={\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T{\left({\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right){\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T+{\mathbf{R}}_{k+1}\right)}^{-1}.\end{array}$

Then, the analysis step is applied for each member in the ensemble cloud as below,

(57)$\begin{array}{cc}\hfill {\mathbf{u}}_a^{\left(i\right)}& ={\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)+{\mathbf{K}}_{k+1}\left({\mathbf{w}}^{\left(i\right)}\left(t_{k+1}\right)-h({\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right))\right),\hfill \end{array}$

and the analyzed estimate at $t_{k+1}$ is computed as the sample mean of the analysis ensemble along with its covariance matrix as follows,

(58)$\begin{array}{cc}\hfill {\mathbf{u}}_a\left(t_{k+1}\right)& \approx {\mathbf{ū}}_a\left(t_{k+1}\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\mathbf{u}}_a^{\left(i\right)}\left(t_{k+1}\right),\hfill \end{array}$

(59)$\begin{array}{cc}\hfill {\xi }_a^{\left(i\right)}\left(t_{k+1}\right)& ={\mathbf{u}}_a^{\left(i\right)}\left(t_{k+1}\right)-{\mathbf{ū}}_a\left(t_{k+1}\right),\hfill \end{array}$

(60)$\begin{array}{cc}\hfill {\mathbf{P}}_{k+1}& \approx {\displaystyle \frac{1}{N-1}}\sum _{i=1}^N{\xi }_a^{\left(i\right)}\left(t_{k+1}\right){\xi }_a^{\left(i\right)}{\left(t_{k+1}\right)}^T.\hfill \end{array}$

We observe that the EnKF algorithm provides approximations of the background and analysis covariance matrices, without the need to evaluate the computationally expensive propagation equations. This comes with the expense of having to evolve an ensemble of system’s states. However, the size of the ensemble N is usually smaller than the system’s dimension n. Moreover, with parallelization and high performance computing (HPC) frameworks, the forecast step can be distributed efficiently and computational speed-ups can be achieved. A summary of the EnKF algorithm is described as follows,   

$\begin{array}{cc}\hfill \mathrm{Inputs}:& {\widehat{\mathbf{u}}}_b\left(t_k\right),{\widehat{\mathbf{B}}}_k,M(·;·),{\mathbf{Q}}_{k+1},\mathbf{w}\left(t_{k+1}\right),{\mathbf{R}}_{k+1},h(·)\hfill \\ \hfill \mathrm{Ensemble}\mathrm{initialization}:& {\widehat{\mathbf{u}}}_b^{\left(i\right)}\left(t_k\right)={\widehat{\mathbf{u}}}_b\left(t_k\right)+{\mathbf{e}}_b^{\left(i\right)},{\mathbf{e}}_b^{\left(i\right)}\sim \mathcal{N}(\mathbf{0},{\widehat{\mathbf{B}}}_k)\hfill \\ \hfill \mathrm{Virtual}\mathrm{observations}:& {\mathbf{w}}^{\left(i\right)}\left(t_{k+1}\right)={\mathbf{w}}^{\left(i\right)}\left(t_{k+1}\right)+{\mathbf{e}}_m^{\left(i\right)},{\mathbf{e}}_m^{\left(i\right)}\sim \mathcal{N}(\mathbf{0},{\mathbf{R}}_{k+1})\hfill \\ \hfill \mathrm{Forecast}:& {\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)=M({\widehat{\mathbf{u}}}_b^{\left(i\right)}\left(t_k\right);\theta )+{\xi }_p^{\left(i\right)}\left(t_{k+1}\right)\hfill \\ \hfill & {\xi }_b^{\left(i\right)}\left(t_{k+1}\right)={\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)-{\mathbf{ū}}_b\left(t_{k+1}\right)\hfill \\ \hfill & {\mathbf{B}}_{k+1}\approx {\displaystyle \frac{1}{N-1}}\sum _{i=1}^N{\xi }_b^{\left(i\right)}\left(t_{k+1}\right){\xi }_b^{\left(i\right)}{\left(t_{k+1}\right)}^T\hfill \\ \hfill \mathrm{Kalman}\mathrm{gain}:& {\mathbf{K}}_{k+1}={\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T{\left({\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}_{k+1}{\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T+{\mathbf{R}}_{k+1}\right)}^{-1}\hfill \\ \hfill \mathrm{Analysis}:& {\mathbf{u}}_a^{\left(i\right)}\left(t_{k+1}\right)={\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)+{\mathbf{K}}_{k+1}\left({\mathbf{w}}^{\left(i\right)}\left(t_{k+1}\right)-h\left({\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)\right)\right)\hfill \\ \hfill & {\mathbf{u}}_a\left(t_{k+1}\right)\approx {\mathbf{ū}}_a\left(t_{k+1}\right)\hfill \\ \hfill & {\xi }_a^{\left(i\right)}\left(t_{k+1}\right)={\mathbf{u}}_a^{\left(i\right)}\left(t_{k+1}\right)-{\mathbf{ū}}_a\left(t_{k+1}\right)\hfill \\ \hfill & {\mathbf{P}}_{k+1}\approx {\displaystyle \frac{1}{N-1}}\sum _{i=1}^N{\xi }_a^{\left(i\right)}\left(t_{k+1}\right){\xi }_a^{\left(i\right)}{\left(t_{k+1}\right)}^T\hfill \end{array}$

and Listing 17 provides a Python execution of the presented EnKF approach.

We highlight a few remarks regarding the EnKF as below,

• Ensemble methods have gained significant popularity because of their simple conceptual formulation and relative ease of implementation. No optimization problem is required to be solved. They are considered non-intrusive in the sense that current solvers can be easily incorporated with minimal modification, as there is no need to derive model Jacobians or adjoint equations.

• The analysis ensemble can be used as initial ensemble for the next assimilation cycle (in which case, we need not compute ${\mathbf{P}}_{k+1}$). Alternatively, new ensemble can be built, by sampling from multivariate Gaussian distribution with a mean of ${\mathbf{u}}_a\left(t_{k+1}\right)$ and covariance matrix of ${\mathbf{P}}_{k+1}$ (i.e., using ${\mathbf{u}}_a\left(t_{k+1}\right)$ and ${\mathbf{P}}_{k+1}$ in lieu of ${\widehat{\mathbf{u}}}_b\left(t_{k+1}\right)$ and ${\widehat{\mathbf{B}}}_{k+1}$, respectively).

• After virtual observations are made-up, an ensemble measurement error covariance matrix can be arbitrarily computed as an alternative to the actual one [17]. This is especially valuable when the actual measurement noise covariance matrix is poorly known.

• Perturbed observations are needed in EnKF derivation and guarantees that the posterior (analysis) covariance is not underestimated. For instance, in case of small corrections to the forecast, the traditional EnKF without virtual observations yields a error covariance that is about twice smaller than that is needed to match Kalman filter [30]. In other words, the use of virtual observations forces the ensemble posterior covariance to be the same as that of the standard Kalman filter in the limit of very large N. Thus, the same Kalman gain matrix relation is borrowed from standard Kalman filter.

• Instead of assuming virtual observations, alternative formulations of ensemble Kalman filters have been proposed in literature, giving a family of deterministic ensemble Kalman filter (DEnKF), as opposed to the aforementioned (stochastic) ensemble Kalman filer (EnKF). One such variant is briefly discussed in Section 8.1.

8.1. Deterministic Ensemble Kalman Filter

The use of an ensemble of perturbed observations in the EnKF leads to a match between the analysis error covariance and its theoretical value given by Kalman filter. However, this is in a statistical sense only when the ensemble size is large. Unfortunately, this perturbation introduces sampling error, which renders the filter suboptimal, particularly for small ensembles [31]. Alternative formulations that do not require virtual observations can be found in literature, including ensemble square root filters [31,32]. We focus here on a simple formulation proposed by Sakov and Oke [30] that maintains the numerical effectiveness and simplicity EnKF without the need to virtual observations, denoted as deterministic ensemble Kalman filter (DEnKF).

Without measurement perturbation, it can be derived that the resulting analysis error covariance matrix is given as follows,

$\begin{array}{cc}\hfill {\mathbf{P}}_{k+1}& =({\mathbf{I}}_n-{\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)){\mathbf{B}}_{k+1}{(\mathbf{I}-\mathbf{K}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right))}^T\hfill \\ \hfill & ={\mathbf{B}}_{k+1}-{\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right){\mathbf{B}}_{k+1}-{\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T{\mathbf{K}}_{k+1}^T\hfill \\ \hfill & +{\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right){\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T{\mathbf{K}}_{k+1}^T.\hfill \end{array}$

With the definition of the Kalman gain, it can be seen that ${\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right){\mathbf{B}}_{k+1}={\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T{\mathbf{K}}_{k+1}^T$. Thus,

$\begin{array}{cc}\hfill {\mathbf{P}}_{k+1}& ={\mathbf{B}}_{k+1}-2{\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right){\mathbf{B}}_{k+1}+{\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right){\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T{\mathbf{K}}_{k+1}^T.\hfill \end{array}$

For small values of ${\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)$, this form converges to ${\mathbf{P}}_{k+1}={\mathbf{B}}_{k+1}-2{\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right){\mathbf{B}}_{k+1}$ up to the quadratic term. It can be seen that this asymptotically match the theoretical value of analysis covariance matrix in standard Kalman filtering (i.e., ${\mathbf{P}}_{k+1}=({\mathbf{I}}_n-{\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)){\mathbf{B}}_{k+1}={\mathbf{B}}_{k+1}-{\mathbf{K}}_{k+1}{\mathbf{D}}_h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)$) by dividing the Kalman gain by two. Therefore, it can be argued that the DEnKF linearly recovers the theoretical analysis error covariance matrix. This is achieved by applying the analysis equation separately to the forecast mean ${\mathbf{u}}_b\left(t_{k+1}\right)\approx {\mathbf{ū}}_b\left(t_{k+1}\right)$ with the Kalman gain matrix and ensemble of anomalies ${\xi }_b^{\left(i\right)}\left(t_{k+1}\right)={\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)-{\mathbf{ū}}_b\left(t_{k+1}\right)$ using half of the standard Kalman gain matrix. These steps are summarized as follows,

$\begin{array}{cc}\hfill \mathrm{Inputs}:& {\widehat{\mathbf{u}}}_b\left(t_k\right),{\widehat{\mathbf{B}}}_k,M(·;·),{\mathbf{Q}}_{k+1},\mathbf{w}\left(t_{k+1}\right),{\mathbf{R}}_{k+1},h(·)\hfill \\ \hfill \mathrm{Ensemble}\mathrm{initialization}:& {\widehat{\mathbf{u}}}_b^{\left(i\right)}\left(t_k\right)={\widehat{\mathbf{u}}}_b\left(t_k\right)+{\mathbf{e}}_b^{\left(i\right)},{\mathbf{e}}_b^{\left(i\right)}\sim \mathcal{N}(\mathbf{0},{\widehat{\mathbf{B}}}_k)\hfill \\ \hfill \mathrm{Forecast}:& {\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)=M({\widehat{\mathbf{u}}}_b^{\left(i\right)}\left(t_k\right);\theta )+{\xi }_p^{\left(i\right)}\left(t_{k+1}\right)\hfill \\ \hfill & {\mathbf{u}}_b\left(t_{k+1}\right)\approx {\mathbf{ū}}_b\left(t_{k+1}\right)\hfill \\ \hfill & {\xi }_b^{\left(i\right)}\left(t_{k+1}\right)={\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)-{\mathbf{ū}}_b\left(t_{k+1}\right)\hfill \\ \hfill & {\mathbf{B}}_{k+1}\approx {\displaystyle \frac{1}{N-1}}\sum _{i=1}^N{\xi }_b^{\left(i\right)}\left(t_{k+1}\right){\xi }_b^{\left(i\right)}{\left(t_{k+1}\right)}^T\hfill \\ \hfill \mathrm{Kalman}\mathrm{gain}:& {\mathbf{K}}_{k+1}={\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T{\left({\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}_{k+1}{\mathbf{B}}_{k+1}{\mathbf{D}}_h{\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)}^T+{\mathbf{R}}_{k+1}\right)}^{-1}\hfill \\ \hfill \mathrm{Analysis}:& {\mathbf{u}}_a\left(t_{k+1}\right)={\mathbf{u}}_b\left(t_{k+1}\right)+{\mathbf{K}}_{k+1}\left(\mathbf{w}\left(t_{k+1}\right)-h\left({\mathbf{u}}_b\left(t_{k+1}\right)\right)\right)\hfill \\ \hfill & {\xi }_a^{\left(i\right)}\left(t_{k+1}\right)={\xi }_b^{\left(i\right)}\left(t_{k+1}\right)-{\displaystyle \frac{1}{2}}{\mathbf{K}}_{k+1}\left(h\left({\xi }_b^{\left(i\right)}\left(t_{k+1}\right)\right)\right)\hfill \\ \hfill & {\mathbf{u}}_a^{\left(i\right)}\left(t_{k+1}\right)={\mathbf{u}}_a\left(t_{k+1}\right)+{\xi }_a^{\left(i\right)}\left(t_{k+1}\right)\hfill \\ \hfill & {\mathbf{P}}_{k+1}\approx {\displaystyle \frac{1}{N-1}}\sum _{i=1}^N{\xi }_a^{\left(i\right)}\left(t_{k+1}\right){\xi }_a^{\left(i\right)}{\left(t_{k+1}\right)}^T\hfill \end{array}$

and Listing 18 provides a Python execution of the presented DEnKF approach. Note that the ensemble of observations is not created in this case, compared to the EnKF.

8.2. Example: Lorenz 63 System

The same Lorenz 63 system is used to showcase the performance of both the EnKF and DEnKF. In Listing 19, we show the Python application of the EnKF algorithm. In general, the size of ensemble N is much smaller than the state dimension n for the implementation of EnKF to be computationally feasible. However, the state dimension in the Lorenz 63 is 3, and an ensemble of size 3 or less is trivial. The uncertainty in the covariance approximation via the Monte Carlo framework with such small ensemble becomes very high and resulting predictions are unreliable. There exists various approaches that help to increase the fidelity of small ensembles, including localization and inflation. In Section 9.2, we describe simple application of inflation factor and its impact with small ensembles. Here, we stick with the basic implementation with an ensemble size of 10 for both EnKF and DEnKF.

Figure 5 shows the EnKF results for the Lorenz 63 system. We see that the analysis trajectory is close to the true one and more accurate than the background. Readers are encouraged to play with the codes to explore the effect of increasing or decreasing the ensemble size with different levels of noise. Furthermore, different observation operators can be defined (for instance, observe only 1 or 2 variables, or assume some nonlinear function $h(·)$).

For the sake of completeness, we also sketch the DEnKF predictions in Figure 6. The same implementation in Listing 19 can be adopted, but calling the DEnKF function framed in Listing 18 instead of EnKF. Although different approaches might give slightly dissimilar results, we are not trying to benchmark them in this introductory presentation since we are only showing very simple implementation, with idealized twin experiments.

9. Applications

In this section, we gradually increase the complexity of the test cases using fluid dynamics applications. In Section 9.1, we slightly increase the dimensionality of the system from 3 (as in Lorenz 63 system) to 36 using the Lorenz 96 system and demonstrate the capability of DA algorithms to treat uncertainty in initial conditions. This is further extended in Section 9.2, where we show that DA can recover the hidden underlying processes and provide closure effects using the two-level variant of the Lorenz 96 model. We also introduce the utilization of an inflation factor and its impact to mitigate ensemble collapse and account for a slight under-representation of covariance due to the use of a small ensemble in EnKF and DEnKF. In Section 9.3, we illustrate the application of DA on systems governed by partial differential equations (PDEs) using the Kuramoto–Sivashinsky equation. We highlight that for each application, we only show results for a few selected algorithms and extensions to other approaches covered in this tutorial are left to readers as computer projects. We also emphasize that we are demonstrating the implementation and capabilities of the presented DA algorithms, not assessing their performance nor benchmarking different approaches against each other.

9.1. Lorenz 96 System

The Lorenz 96 model [33] is a system of ordinary differential equations that describes an arbitrary atmospheric quantity as it evolves on a circular array of sites, undergoing forcing, dissipation, and rotation invariant advection [34]. The Lorenz 96 dynamical model can be written as

(61)${\displaystyle \frac{\mathrm{d}{X}_i}{\mathrm{d}t}}=(X_{i+1}-X_{i-2})X_{i-1}-X_i+F,i=1,2,\dots ,n,$

where $X_i$ is the state of the system at the $i^{th}$ location and F represents a forcing constant. Periodicity is enforced by assuming that $X_{-1}=X_{n-1}$, $X_0=X_n$, and $X_{n+1}=X_1$. In the present study, we use $n=36$, and $F=8$ defining a forcing term. In order to obtain a valid initial condition, we begin at $t=-5$ using equilibrium conditions defined as ($X_i=F$ for $i=1,2,\dots ,n$) and adding a small perturbation to the 20th state variable as $X_{20}=F+0.01$. Then, ODE integrator is run up to $t=0$ and solution at $t=0$ is treated at the true initial conditions for our twin experiment. We assume a background erroneous initial condition by contaminating the true one with Gaussian noise with zero mean and standard deviation of 1.

A total time window of 20 time units is considered, with a time step of $\mathsf{\Delta }t=0.01$ and the RK4 schemes is adopted for time integration. Synthetic measurements are collected at every $0.2$ time unit (i.e., each 20 time integration steps) sampled at 9 equidistant locations (i.e., at $i\in \{4,8,12,\dots ,36\}$) from true trajectory assuming that sensors add a white noise with zero mean and a standard deviation of $0.1$. We also assume a process noise drawn from a multivariate Gaussian distribution with zero mean and covariance matrix $\mathbf{Q}$ defined as $\mathbf{Q}=0.1^2{\mathbf{I}}_{36}$. We first apply the EKF approach to correct the solution trajectory, which yields very good results as shown in Figure 7. For visualization, we only plot the time evolution of $X_9$, $X_{18}$, and $X_{36}$. We see that the Lorenz 96 is sensitive to the initial conditions and small perturbation is sufficient to produce a very different trajectory (e.g., background solution).

The second approach to test is the stochastic version of EnKF. We create an ensemble of 50 members to approximate the covariance matrices. Results are depicted in Figure 8 for $X_9$, $X_{18}$, and $X_{36}$. We highlight that observations appear only in the $X_{36}$ plot because observations are collected at $i=4,8,12,\dots ,36$ (neither $X_9$ nor $X_{18}$ are measured). We highlight that, generally speaking, increasing the size of ensemble improves the predictions. However, this comes on the expense of the solution of the forward nonlinear model for each added member. Thus, a compromise between the accuracy and computational burden is in place.

9.2. Two-Level Lorenz 96 System

In this section, we describe the two-level variant of the Lorenz 96 model proposed by Lorenz [33]. The two-level Lorenz 96 model can be written as

(62)$\begin{array}{cc}\hfill \frac{d{X}_i}{dt}& =-X_{i-1}(X_{i-2}-X_{i+1})-X_i-\frac{hc}{b}\sum _{j=1}^J{Y}_{j,i}+F,\hfill \end{array}$

(63)$\begin{array}{cc}\hfill \frac{d{Y}_{j,i}}{dt}& =-cb{Y}_{j+1,i}(Y_{j+2,i}-Y_{j-1,i})-c{Y}_{j,i}+\frac{hc}{b}{X}_i,\hfill \end{array}$

where Equation (62) represents the evolution of slow, high-amplitude variables $X_i(i=1,\dots ,n)$, and Equation (63) provides the evolution of a coupled fast, low-amplitude variable $Y_{j,i}(j=1,\dots ,J)$. We use $n=36$ and $J=10$ in our computational experiments. We utilize $c=10$ and $b=10$, which implies that the small scales fluctuate 10 times faster than the larger scales. Furthermore, the coupling coefficient h between two scales is equal to 1 and the forcing is set at $F=10$ to make both variables exhibit the chaotic behavior.

We utilize the fourth-order Runge–Kutta numerical scheme with a time step $\mathsf{\Delta }t=0.001$ for temporal integration of the Lorenz 96 model. We apply the periodic boundary condition for the slow variables, i.e., $X_{i-n}=X_{i+n}=X_i$. The fast variables are extended by letting $Y_{j,i-n}=Y_{j,i+n}=Y_{j,i}$, $Y_{j-J,i}=Y_{j,i-1}$, and $Y_{j+J,i}=Y_{j,i+1}$. The physical initial condition is computed by starting with an equilibrium condition at time $t=-5$ for slow variables. The equilibrium condition for slow variables is $X_i=F$ for $i\in 1,2,\dots ,n$. We perturb the equilibrium solution for the 18th state variable as $X_{18}=F+0.01$. At the time $t=-5$, the fast variables are assigned with random numbers between $-F/10$ to $F/10$. We integrate a two-level Lorenz 96 model by solving both Equations (62) and (63) in a coupled manner up to time $t=0$. The solution at time $t=0$ represent the true initial condition. For our twin experiment, we obtain observations by adding noise drawn from the Gaussian distribution with zero mean and ${\sigma }_o^2=1.0$. We assume that observations are sparse in space and are collected at every 10th time step.

The motivation behind this example is to demonstrate how covariance inflation can be utilized to account for the model error. Usually the imperfections in the forecast model is taken into account by adding a Gaussian noise to the forecast model. Another method to account for model error is covarinace inflation. It also helps in alleviating the effect finite number of ensemble members in practical data assimilation and addresses the problem of covariance underestimation in the EnKF algorithm. We use the multiplicative inflation [35] where the ensemble members are pushed away from the ensemble mean by a given inflation factor and mathematically it can be expressed as

(64)$\begin{array}{c}\hfill {\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)\leftarrow {\mathbf{u}}_a\left(t_{k+1}\right)+\lambda ·({\mathbf{u}}_b^{\left(i\right)}\left(t_{k+1}\right)-{\mathbf{u}}_a^{\left(i\right)}\left(t_{k+1}\right)),\end{array}$

where $\lambda$ is the inflation factor. The inflation factor can be a constant scalar over the entire domain at all time step or it can space and time dependent.

In this example, we discard parameterizations of fast variables in the forecast model. The forecast model for two-level Lorenz system with no parameterizations is equivalent to setting the coupling coefficient $h=0$ in Equation (62) and it reduces to one-level Lorenz 96 model as presented in Section 9.1. We note here that the observations used for data assimilation are obtained by solving a two-level Lorenz 96 model in a coupled manner (i.e., without discarding fast-variables). Therefore, the effect of unresolved scales is embedded in observations. The parameterization of fast variables (i.e.,$\frac{hc}{b}{\sum }_{j=1}^J{Y}_{j,i}$ term in Equation (62)) can be considered as an added noise to the true state of the system for a one-level Lorenz 96 model. Figure 9 displays the RMSE for a two-level Lorenz system when 18 observations are used for DA with different number of ensemble members and inflation factors for EnKF and DEnKF algorithms. Figure 10 shows the full state trajectory of two-level Lorenz system corresponding to minimum RMSE, which is obtained with 50 ensemble members and inflation factor $\lambda =1.04$ for the EnKF algorithm. The parameters corresponding to minimum RMSE for the DEnKF algorithm are 45 ensemble members and $\lambda =1.05$.

9.3. Kuramato Sivashinsky

In this section, we describe the Kuramoto–Sivashinsky (K-S) equation derived by Kuramoto [36], which is used as a turbulence model for different flows. The one-dimensional K-S equation can be written as

(65)$\frac{\partial u}{\partial t}=-\nu \frac{{\partial }^{4}u}{\partial x^4}-\frac{{\partial }^{2}u}{\partial x^2}-u\frac{\partial u}{\partial x},$

where $\nu$ is the viscosity coefficient. The K-S equation is characterized by the second-order unstable diffusion term responsible for an instability at large scales, the fourth-order stabilizing viscous term that provides damping at small scales, and a quadratic nonlinear term which transfers energy between large and small scales. We use the computational domain extending from 0 to L, i.e., $x\in [0,L]$ and time $t\in [0,\infty ]$. We impose the Dirichlet and Neumann boundary conditions as given below

(66)$\begin{array}{c}\hfill u(0,t)=u(L,t)=0,\end{array}$

(67)$\begin{array}{c}\hfill \frac{\partial u}{\partial x}{|}_{x=0}=\frac{\partial u}{\partial x}{|}_{x=L}=0.\end{array}$

We spatially discretize the domain with the grid size $\mathsf{\Delta }x=L/(n-1)$, where n is the degrees of freedom. We set $L=50$ and $n=129$ for our numerical experiments. The state of the system at discretized grid is denoted as $u_i=u\left((i-1)\mathsf{\Delta }x\right)$ for $i=1,\dots ,n$. Using the second-order finite difference discretization, the discretized K-S equation can be written as

(68)$\begin{array}{c}\hfill \frac{d{u}_i}{dt}=-\nu \frac{u_{i+2}-4u_{i+1}+6u_i-4u_{i-1}+u_{i-2}}{\mathsf{\Delta }{x}^4}-\frac{u_{i+1}-2u_i+u_{i-1}}{\mathsf{\Delta }{x}^2}-\frac{1}{2}\frac{u_{i+1}^2-u_{i-1}^2}{2\mathsf{\Delta }x}.\end{array}$

The first term on the right hand side is computed by utilizing ghost nodes and the Neumann boundary condition is assigned for ghost points. We impose $u_0=u_2$ and $u_{n+1}=u_{n-1}$ using the second-order discretization for Equation (67) at boundary points $u_1$ and $u_n$, respectively.

We use the fourth-order Runge–Kutta scheme for time integration with a time step $\mathsf{\Delta }t=0.25$. To generate an initial condition for the forward run we start with an equilibrium condition at time $t=-50$ and integrate up to time $t=0$. The equilibrium condition for the model is $u_i=0.1$ for $j\in \{1,\dots ,n\}$. Once the true initial condition is generated, we run the forward solver up to time $t=50$. We test the prediction capability of sequential data assimilation algorithms for forecast up to $t=50$. The K-S equation exhibits different levels of chaotic behavior depending on the value of viscosity coefficient $\nu$. The chaos depends upon the bifurcation parameter $\tilde{L}=L/2\pi \sqrt{\nu }$. We utilize $\nu =1/2$ which represent the less chaotic behaviour. The observations for twin experiments are obtained by adding some noise to the true state of the system to account for experimental uncertainties and measurement errors. The observations are also sparse in time, meaning that the time interval between two observations can be different from the time step of the forecast model. For our twin experiments, we assume that observations are recorded at every 10th time step of the model for $\nu =1/2$. Therefore, the time difference between two observations is $\delta t=2.5$ for the K-S equation. We present results for the EnKF algorithm with three sets of observations. The first set of observations is very sparse with only 12.5% of the full state of the system. The first set utilizes observations for states $[u_8,u_{16},\dots ,u_{128}]\in R^{16}$. In a second set of observations we employ observations at $[u_4,u_8,\dots ,u_{128}]\in R^{32}$ for the assimilation. The third set of observations consists of 50% of the full state of the system, i.e., observations at states $[u_2,u_4,\dots ,u_{128}]\in R^{64}$ for the assimilation. We apply ${\sigma }_o^2=1.0\times {10}^{-2}$ and ${\sigma }_i^2=1.0\times {10}^{-2}$ as the variance of observation noise and initial condition uncertainty, respectively.

In Figure 11, we present the time evolution of selected states for three different number of observations included in the assimilation of the EnKF algorithm. There is an excellent agreement between true and assimilated states $u_{51}$ and $u_{101}$, for which observations are not present. We also provide the full state trajectory of the K-S equation in Figure 12. The results obtained clearly indicate that the EnKF algorithm is able to determine the correct state trajectory even when the observation data are very sparse, i.e., $m=16$. With an increase in the number of observations, the prediction of the full state trajectory gets smoother, and almost the exact state is recovered with 50% observations.

9.4. Quasi-Geostrophic (QG) Ocean Circulation Model

We consider a simple single-layer QG model to illustrate the application of sequential data assimilation for two-dimensional flows. Specifically, we use the deterministic ensemble Kalman filter (DenKF) algorithm discussed in Section 8.1 to improve the prediction of the single-layer QG model. The wind-driven oceanic flows exhibit a vast range of spatio-temporal scales and modeling of these scales with all the relevant physics has always been challenging. The barotropic vorticity equation (BVE) with various dissipative and forcing terms is one of the most commonly used models for geostrophic flows [37,38]. The dimensionless vorticity-streamfunction formulation for the BVE [39] with forcing and dissipative terms can be written as

(69)$\frac{\partial \omega }{\partial t}+J(\omega ,\psi )-\frac{1}{Ro}\frac{\partial \psi }{\partial x}=\frac{1}{Re}{\mathsf{\nabla }}^2\omega +\frac{1}{Ro}sin\left(\pi y\right),$

where $\omega$ is the vorticity, $\psi$ is the streamfunction, ${\mathsf{\nabla }}^2$ is the standard two-dimensional Laplacian operator, Re is the Reynolds number, and $Ro$ is the Rossby number. The kinematic relation between vorticity and streamfunction is given by the following Poisson equation

(70)${\mathsf{\nabla }}^2\psi =-\omega .$

The nonlinear convection term is given by the Jacobian as follows

(71)$J(\omega ,\psi )=\frac{\partial \psi }{\partial y}\frac{\partial \omega }{\partial x}-\frac{\partial \psi }{\partial x}\frac{\partial \omega }{\partial y}.$

The computational domain for the QG model is $(x,y)\in [0,1]\times [-1,1]$ and is discretized using $128\times 256$ grid resolution. Therefore, the QG model has the dimension of about $3.2\times {10}^4$. We utilize the homogeneous Dirichlet boundary condition for the vorticity and streamfunction at all boundaries. The vorticity and streamfunction is initialized from quiescent state, i.e., ${\omega }_{t=0}{=\psi |}_{t=0}=0$. The QG model is numerically solved by discretizing Equation (69) using second-order finite difference scheme. The nonlinear Jacobian term is discretized with the energy-conserving Arakawa [40] numerical scheme. A third-order total-variation-diminishing Runge–Kutta scheme is used for the temporal integration and a fast sine transform Poisson solver is utilized to update streamfunction from the vorticity [41]. For the physical parameters, we use values of $\mathrm{Re}=100$ and $\mathrm{Ro}=1.75\times {10}^{-3}$.

The QG model is integrated with a constant time step of $5\times {10}^{-5}$ from time $t=0$ to $t=0.25$ to generate the true initial condition at final time $t=0.25$. Then the data assimilation is conducted from time $t=0.25$ to $t=0.4$ with observations getting assimilated at every tenth time step. The synthetic observations are generated by sampling vorticity field on 16 equidistant points in x and y directions respectively and then adding the Gaussian noise, i.e., ${\mathbf{v}}_k\sim \mathcal{N}(0,{\mathbf{R}}_k)$, where ${\mathbf{R}}_k={\sigma }_b^2\mathbf{I}$. We set the observation noise variance at ${\sigma }_b^2=5$. The typical vorticity and streamfunction field along with the locations of measurements are shown in Figure 13.

We employ 20 ensemble members for the DEnKF algorithm. The initialization of the ensemble members is an important step to get accurate prediction with any type of the EnKF algorithm. We initialize different ensemble members by randomly selecting the vorticity field snapshots between time $t=0.24$ to $t=0.25$. The other methods such as adding a random perturbation from the Gaussian distribution to the true initial condition can also be adopted. Figure 14 displays the vorticity field and the predicted vorticity field at three different time instances along with the difference (error) between the two. We can see that the true and analysis field are similar at all time instances and the magnitude of error is also small. We recall that we observe only around 2% of the system (i.e., observations at $16\times 32$ locations). With more observations, the quality of the results can be further improved.

10. Concluding Remarks

In this tutorial paper, we provided a 101 introduction to common data assimilation techniques. In particular, we briefly covered the relevant mathematical foundation and the algorithmic steps for three dimensional variational (3DVAR), four dimensional variational (4DVAR), forward sensitivity method (FSM), and Kalman filtering approaches. Since it is considered as a first exposure to DA, we focused on the simplest implementations that anybody can easily follow. For example, to treat nonlinearity (e.g., in 3DVAR, 4DVAR, FSM and EKF), we only presented the first order Taylor expansions. We demonstrated the execution of the covered approaches with a series of Python modules that can be linked to each other easily. Again, we preferred to keep our codes as concise and simple as possible, even if it comes on the expense of computational efficiency. The Python codes used to generate this tutorial are publicly available through our GitHub repository https://github.com/Shady-Ahmed/PyDA.

Since it is introductory exploration, we should admit that we have bypassed a few important analyses and shortcut some key derivations. Interested readers are referred to well-established textbooks that offer in-depth discussions about various DA techniques [14,27,29,42,43,44,45,46,47]. Likewise, more advanced topics such as the particle filters [48], maximum likelihood ensemble filters [49,50], optimal sensor placement [51,52], higher-order analysis of variational methods [53], or hybrid methods [54,55,56,57,58,59] are omitted in our current presentation.

Author Contributions

Data curation, S.E.A., S.P., and O.S.; Supervision, O.S.; Writing—original draft, S.E.A., and S.P.; and Writing—review and editing, S.E.A., S.P., and O.S. All authors have read and agreed to the published version of the manuscript.

Funding

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under Award Number DE-SC0019290. Omer San gratefully acknowledges their support. Disclaimer: This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Acknowledgments

We thank Sivaramakrishnan Lakshmivarahan for his insightful comments as well as his archival NPTEL lectures [60] on dynamic data assimilation that greatly helped us in developing the PyDA module. Special thanks go to Ionel Michael Navon for providing his lecture notes on data assimilation.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Listings

Figure 1. Results of 3DVAR implementation for the Lorenz 63 system.

Figure 2. Results of 4DVAR implementation for the Lorenz 63 system.

Figure 3. Results of FSM implementation for the Lorenz 63 system.

Figure 4. EKF results for the Lorenz 63 system with the assumption of zero process noise.

Figure 5. EnKF results for the Lorenz 63 system with virtual observations.

Figure 6. DEnKF results for the Lorenz 63 system without virtual observations.

Figure 7. EKF results for the Lorenz 96 system. The trajectories of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] are shown.

Figure 8. EnKF results for the Lorenz 96 system using an ensemble of 50 members.

Figure 9. RMSE for a two-level Lorenz model for different combinations of number of ensembles and inflation factor.

Figure 10. Full state trajectory of the multiscale Lorenz 96 model with no parameterizations in the forecast model. The EnKF algorithm uses the inflation factor [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] and the DEnKF uses the inflation factor [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. The observation data for both EnKF an DEnKF algorithm is obtained by adding measurement noise to the exact solution of the two-level Lorenz 96 system.

Figure 11. Selected trajectories of the Kuramoto–Sivashinsky model ([Forumla omitted. See PDF.]) with the analysis performed by the ensemble Kalman filter (EnKF) using observations from [Forumla omitted. See PDF.] (left), [Forumla omitted. See PDF.] (middle), and [Forumla omitted. See PDF.] (right) state variables at every 10 time steps.

Figure 12. Full state trajectory of the Kuramoto–Sivashinsky model ([Forumla omitted. See PDF.]) with the analysis performed by the ensemble Kalman filter (EnKF) using observations from [Forumla omitted. See PDF.] (left), [Forumla omitted. See PDF.] (middle), and [Forumla omitted. See PDF.] (right) state variables at every 10 time steps.

Figure 13. A Typical vorticity (left) and streamfunction (right) field for the single-layer QG model. The dots shows the locations of observations.

Figure 14. Snapshots of the true vorticity field (left), analysis estimate of the DEnKF algorithm (middle), and the difference (error) between the two fields (right) obtained for a particular run of the single-layer QG model. The snapshots of vorticity field are plotted at [Forumla omitted. See PDF.] (from top to bottom).

Listing 1. Implementation of 3DVAR for with a linear observation operator.

Listing 2. Implementation of the 3DVAR for with a nonlinear observation operator, usingfirst-order approximation.

Listing 3. A Python function for the Lorenz 63 dynamics.

Listing 4. Python functions for the time integration using the 1st Euler and the 4thRunge–Kutta schemes.

Listing 5. Implementation of the 3DVAR for the Lorenz 63 system.

Listing 5. Implementation of the 3DVAR for the Lorenz 63 system.

Listing 6. Computation of the gradient of the cost functional with the 4DVAR using the first-orderadjoint algorithm.

Listing 6. Computation of the gradient of the cost functional with the 4DVAR using the first-orderadjoint algorithm.

Listing 6. Computation of the gradient of the cost functional with the 4DVAR using the first-orderadjoint algorithm.

Listing 7. A line-search Python function using the Golden search method.

Listing 8. Computation of the cost functional defined in Equation (13).

Listing 8. Computation of the cost functional defined in Equation (13).

Listing 9. Implementation of 4DVAR using the first-order adjoint method for the Lorenz 63 system.

Listing 9. Implementation of 4DVAR using the first-order adjoint method for the Lorenz 63 system.

Listing 9. Implementation of 4DVAR using the first-order adjoint method for the Lorenz 63 system.

Listing 10. Python functions for computing the Jacobian [Forumla omitted. See PDF.] of the discrete-time model map [Forumla omitted. See PDF.] using the 1st Euler and the 4th Runge–Kutta schemes with chain rule.

Listing 11. A Python function for the Jacobian of the continuous-time Lorenz 63 dynamics.

Listing 12. Python function for computing the correction vector [Forumla omitted. See PDF.] sing the forward sensitivity method with first-order approximation.

Listing 12. Python function for computing the correction vector [Forumla omitted. See PDF.] sing the forward sensitivity method with first-order approximation.

Listing 13. Implementation of the FSM for the Lorenz 63 system.

Listing 13. Implementation of the FSM for the Lorenz 63 system.

Listing 14. Implementation of the KF with linear dynamics and observation operator.

Listing 14. Implementation of the KF with linear dynamics and observation operator.

Listing 15. Implementation of the (first-order) EKF with nonlinear dynamics and nonlinearobservation operator.

Listing 16. Implementation of the EKF for the Lorenz 63 system.

Listing 17. Implementation of the EnKF with virtual observations.

Listing 18. Implementation of DEnKF without virtual observations.

Listing 18. Implementation of DEnKF without virtual observations.

Listing 19. Implementation of EnKF for the Lorenz 63 system.

Listing 19. Implementation of EnKF for the Lorenz 63 system.

References

1. Navon, I.M. Data assimilation for numerical weather prediction: A review. Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications; Springer: Berlin/Heidelberg, Germany, 2009; pp. 21-65.

2. Blum, J.; Le Dimet, F.X.; Navon, I.M. Data assimilation for geophysical fluids. Handb. Numer. Anal.; 2009; 14, pp. 385-441.

3. Le Dimet, F.X.; Navon, I.M.; Ştefănescu, R. Variational data assimilation: Optimization and optimal control. Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications (Vol. III); Springer: Berlin/Heidelberg, Germany, 2017; pp. 1-53.

4. Attia, A.; Sandu, A. DATeS: A highly extensible data assimilation testing suite v1.0. Geosci. Model Dev.; 2019; 12, pp. 629-649. [DOI: https://dx.doi.org/10.5194/gmd-12-629-2019]

5. Lorenc, A.C. Analysis methods for numerical weather prediction. Q. J. R. Meteorol. Soc.; 1986; 112, pp. 1177-1194. [DOI: https://dx.doi.org/10.1002/qj.49711247414]

6. Parrish, D.F.; Derber, J.C. The National Meteorological Center’s spectral statistical-interpolation analysis system. Mon. Weather Rev.; 1992; 120, pp. 1747-1763. [DOI: https://dx.doi.org/10.1175/1520-0493(1992)120<1747:TNMCSS>2.0.CO;2]

7. Courtier, P. Dual formulation of four-dimensional variational assimilation. Q. J. R. Meteorol. Soc.; 1997; 123, pp. 2449-2461. [DOI: https://dx.doi.org/10.1002/qj.49712354414]

8. Rabier, F.; Järvinen, H.; Klinker, E.; Mahfouf, J.F.; Simmons, A. The ECMWF operational implementation of four-dimensional variational assimilation. I: Experimental results with simplified physics. Q. J. R. Meteorol. Soc.; 2000; 126, pp. 1143-1170. [DOI: https://dx.doi.org/10.1002/qj.49712656415]

9. Elbern, H.; Schmidt, H.; Talagrand, O.; Ebel, A. 4D-variational data assimilation with an adjoint air quality model for emission analysis. Environ. Model. Softw.; 2000; 15, pp. 539-548. [DOI: https://dx.doi.org/10.1016/S1364-8152(00)00049-9]

10. Courtier, P.; Thépaut, J.N.; Hollingsworth, A. A strategy for operational implementation of 4D-Var, using an incremental approach. Q. J. R. Meteorol. Soc.; 1994; 120, pp. 1367-1387. [DOI: https://dx.doi.org/10.1002/qj.49712051912]

11. Lorenc, A.C.; Rawlins, F. Why does 4D-Var beat 3D-Var?. Q. J. R. Meteorol. Soc.; 2005; 131, pp. 3247-3257. [DOI: https://dx.doi.org/10.1256/qj.05.85]

12. Gauthier, P.; Tanguay, M.; Laroche, S.; Pellerin, S.; Morneau, J. Extension of 3DVAR to 4DVAR: Implementation of 4DVAR at the Meteorological Service of Canada. Mon. Weather Rev.; 2007; 135, pp. 2339-2354. [DOI: https://dx.doi.org/10.1175/MWR3394.1]

13. Lakshmivarahan, S.; Lewis, J.M. Forward sensitivity approach to dynamic data assimilation. Adv. Meteorol.; 2010; 2010, 375615. [DOI: https://dx.doi.org/10.1155/2010/375615]

14. Lakshmivarahan, S.; Lewis, J.M.; Jabrzemski, R. Forecast Error Correction Using Dynamic Data Assimilation; Springer: Cham, Switzerland, 2017.

15. Houtekamer, P.L.; Mitchell, H.L. Data assimilation using an ensemble Kalman filter technique. Mon. Weather Rev.; 1998; 126, pp. 796-811. [DOI: https://dx.doi.org/10.1175/1520-0493(1998)126<0796:DAUAEK>2.0.CO;2]

16. Burgers, G.; Jan van Leeuwen, P.; Evensen, G. Analysis scheme in the ensemble Kalman filter. Mon. Weather Rev.; 1998; 126, pp. 1719-1724. [DOI: https://dx.doi.org/10.1175/1520-0493(1998)126<1719:ASITEK>2.0.CO;2]

17. Evensen, G. The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean. Dyn.; 2003; 53, pp. 343-367. [DOI: https://dx.doi.org/10.1007/s10236-003-0036-9]

18. Houtekamer, P.L.; Mitchell, H.L. A sequential ensemble Kalman filter for atmospheric data assimilation. Mon. Weather Rev.; 2001; 129, pp. 123-137. [DOI: https://dx.doi.org/10.1175/1520-0493(2001)129<0123:ASEKFF>2.0.CO;2]

19. Houtekamer, P.L.; Mitchell, H.L. Ensemble kalman filtering. Q. J. R. Meteorol. Soc.; 2005; 131, pp. 3269-3289. [DOI: https://dx.doi.org/10.1256/qj.05.135]

20. Treebushny, D.; Madsen, H. A new reduced rank square root Kalman filter for data assimilation in mathematical models. Proceedings of the International Conference on Computational Science; Melbourne, Australia, 2 June 2003; Springer: Berlin/Heidelberg, Germany, 2003; pp. 482-491.

21. Buehner, M.; Malanotte-Rizzoli, P. Reduced-rank Kalman filters applied to an idealized model of the wind-driven ocean circulation. J. Geophys. Res. Ocean.; 2003; 108, [DOI: https://dx.doi.org/10.1029/2001JC000873]

22. Lakshmivarahan, S.; Stensrud, D.J. Ensemble Kalman filter. IEEE Control. Syst. Mag.; 2009; 29, pp. 34-46.

23. Apte, A.; Hairer, M.; Stuart, A.; Voss, J. Sampling the posterior: An approach to non-Gaussian data assimilation. Phys. Nonlinear Phenom.; 2007; 230, pp. 50-64. [DOI: https://dx.doi.org/10.1016/j.physd.2006.06.009]

24. Bocquet, M.; Pires, C.A.; Wu, L. Beyond Gaussian statistical modeling in geophysical data assimilation. Mon. Weather Rev.; 2010; 138, pp. 2997-3023. [DOI: https://dx.doi.org/10.1175/2010MWR3164.1]

25. Vetra-Carvalho, S.; Van Leeuwen, P.J.; Nerger, L.; Barth, A.; Altaf, M.U.; Brasseur, P.; Kirchgessner, P.; Beckers, J.M. State-of-the-art stochastic data assimilation methods for high-dimensional non-Gaussian problems. Tellus Dyn. Meteorol. Oceanogr.; 2018; 70, pp. 1-43. [DOI: https://dx.doi.org/10.1080/16000870.2018.1445364]

26. Attia, A.; Moosavi, A.; Sandu, A. Cluster sampling filters for non-Gaussian data assimilation. Atmosphere; 2018; 9, 213. [DOI: https://dx.doi.org/10.3390/atmos9060213]

27. Lewis, J.M.; Lakshmivarahan, S.; Dhall, S. Dynamic Data Assimilation: A Least Squares Approach; Cambridge University Press: Cambridge, UK, 2006; Volume 104.

28. Evensen, G. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res. Ocean.; 1994; 99, pp. 10143-10162. [DOI: https://dx.doi.org/10.1029/94JC00572]

29. Evensen, G. Data Assimilation: The Ensemble Kalman Filter; Springer: Berlin/Heidelberg, Germany, 2009.

30. Sakov, P.; Oke, P.R. A deterministic formulation of the ensemble Kalman filter: An alternative to ensemble square root filters. Tellus Dyn. Meteorol. Oceanogr.; 2008; 60, pp. 361-371. [DOI: https://dx.doi.org/10.1111/j.1600-0870.2007.00299.x]

31. Whitaker, J.S.; Hamill, T.M. Ensemble data assimilation without perturbed observations. Mon. Weather Rev.; 2002; 130, pp. 1913-1924. [DOI: https://dx.doi.org/10.1175/1520-0493(2002)130<1913:EDAWPO>2.0.CO;2]

32. Tippett, M.K.; Anderson, J.L.; Bishop, C.H.; Hamill, T.M.; Whitaker, J.S. Ensemble square root filters. Mon. Weather Rev.; 2003; 131, pp. 1485-1490. [DOI: https://dx.doi.org/10.1175/1520-0493(2003)131<1485:ESRF>2.0.CO;2]

33. Lorenz, E.N. Predictability: A problem partly solved. Proceedings of the Seminar on Predictability; Reading, UK, 9–11 September 1996; Volume 1.

34. Kerin, J.; Engler, H. On the Lorenz’96 Model and Some Generalizations. arXiv; 2020; arXiv: 2005.07767

35. Anderson, J.L.; Anderson, S.L. A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts. Mon. Weather Rev.; 1999; 127, pp. 2741-2758. [DOI: https://dx.doi.org/10.1175/1520-0493(1999)127<2741:AMCIOT>2.0.CO;2]

36. Kuramoto, Y. Diffusion-induced chaos in reaction systems. Prog. Theor. Phys. Suppl.; 1978; 64, pp. 346-367. [DOI: https://dx.doi.org/10.1143/PTPS.64.346]

37. Majda, A.; Wang, X. Nonlinear Dynamics and Statistical Theories for basic Geophysical Flows; Cambridge University Press: New York, NY, USA, 2006.

38. Greatbatch, R.J.; Nadiga, B.T. Four-gyre circulation in a barotropic model with double-gyre wind forcing. J. Phys. Oceanogr.; 2000; 30, pp. 1461-1471. [DOI: https://dx.doi.org/10.1175/1520-0485(2000)030<1461:FGCIAB>2.0.CO;2]

39. San, O.; Staples, A.E.; Wang, Z.; Iliescu, T. Approximate deconvolution large eddy simulation of a barotropic ocean circulation model. Ocean. Model.; 2011; 40, pp. 120-132. [DOI: https://dx.doi.org/10.1016/j.ocemod.2011.08.003]

40. Arakawa, A. Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I. J. Comput. Phys.; 1997; 135, pp. 103-114. [DOI: https://dx.doi.org/10.1006/jcph.1997.5697]

41. Press, W.H.; Flannery, B.P.; Teukolsky, S.A.; Vetterling, W.T. Numerical Recipes; Cambridge University Press: New York, NY, USA, 1989.

42. Cacuci, D.G.; Navon, I.M.; Ionescu-Bujor, M. Computational Methods for Data Evaluation and Assimilation; CRC Press: New York, NY, USA, 2013.

43. Kalnay, E. Atmospheric Modeling, Data Assimilation and Predictability; Cambridge University Press: New York, NY, USA, 2003.

44. Law, K.; Stuart, A.; Zygalakis, K. Data Assimilation: A Mathematical Introduction; Springer: Cham, Switzerland, 2015.

45. Asch, M.; Bocquet, M.; Nodet, M. Data Assimilation: Methods, Algorithms, and Applications; SIAM: Philadelphia, PA, USA, 2016.

46. Simon, D. Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches; John Wiley & Sons: Hoboken, NJ, USA, 2006.

47. Labbe, R. Kalman and bayesian filters in Python. Chap; 2014; 7, 246.

48. Van Leeuwen, P.J.; Künsch, H.R.; Nerger, L.; Potthast, R.; Reich, S. Particle filters for high-dimensional geoscience applications: A review. Q. J. R. Meteorol. Soc.; 2019; 145, pp. 2335-2365. [DOI: https://dx.doi.org/10.1002/qj.3551] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31598012]

49. Zupanski, M. Maximum likelihood ensemble filter: Theoretical aspects. Mon. Weather Rev.; 2005; 133, pp. 1710-1726. [DOI: https://dx.doi.org/10.1175/MWR2946.1]

50. Zupanski, M.; Navon, I.M.; Zupanski, D. The Maximum Likelihood Ensemble Filter as a non-differentiable minimization algorithm. Q. J. R. Meteorol. Soc.; 2008; 134, pp. 1039-1050. [DOI: https://dx.doi.org/10.1002/qj.251]

51. Kang, W.; Xu, L. Optimal placement of mobile sensors for data assimilations. Tellus Dyn. Meteorol. Oceanogr.; 2012; 64, 17133. [DOI: https://dx.doi.org/10.3402/tellusa.v64i0.17133]

52. Mons, V.; Chassaing, J.C.; Sagaut, P. Optimal sensor placement for variational data assimilation of unsteady flows past a rotationally oscillating cylinder. J. Fluid Mech.; 2017; 823, pp. 230-277. [DOI: https://dx.doi.org/10.1017/jfm.2017.313]

53. Le Dimet, F.X.; Navon, I.M.; Daescu, D.N. Second-order information in data assimilation. Mon. Weather Rev.; 2002; 130, pp. 629-648. [DOI: https://dx.doi.org/10.1175/1520-0493(2002)130<0629:SOIIDA>2.0.CO;2]

54. Lorenc, A.C.; Bowler, N.E.; Clayton, A.M.; Pring, S.R.; Fairbairn, D. Comparison of hybrid-4DEnVar and hybrid-4DVar data assimilation methods for global NWP. Mon. Weather Rev.; 2015; 143, pp. 212-229. [DOI: https://dx.doi.org/10.1175/MWR-D-14-00195.1]

55. Desroziers, G.; Camino, J.T.; Berre, L. 4DEnVar: Link with 4D state formulation of variational assimilation and different possible implementations. Q. J. R. Meteorol. Soc.; 2014; 140, pp. 2097-2110. [DOI: https://dx.doi.org/10.1002/qj.2325]

56. Wang, X.; Barker, D.M.; Snyder, C.; Hamill, T.M. A hybrid ETKF–3DVAR data assimilation scheme for the WRF model. Part I: Observing system simulation experiment. Mon. Weather Rev.; 2008; 136, pp. 5116-5131. [DOI: https://dx.doi.org/10.1175/2008MWR2444.1]

57. Buehner, M.; Morneau, J.; Charette, C. Four-dimensional ensemble-variational data assimilation for global deterministic weather prediction. Nonlinear Process. Geophys.; 2013; 20, pp. 669-682. [DOI: https://dx.doi.org/10.5194/npg-20-669-2013]

58. Kleist, D.T.; Ide, K. An OSSE-based evaluation of hybrid variational-ensemble data assimilation for the NCEP GFS. Part I: System description and 3D-hybrid results. Mon. Weather Rev.; 2015; 143, pp. 433-451. [DOI: https://dx.doi.org/10.1175/MWR-D-13-00351.1]

59. Kleist, D.T.; Ide, K. An OSSE-based evaluation of hybrid variational-ensemble data assimilation for the NCEP GFS. Part II: 4DEnVar and hybrid variants. Mon. Weather Rev.; 2015; 143, pp. 452-470. [DOI: https://dx.doi.org/10.1175/MWR-D-13-00350.1]

60. Lakshmivarahan, S. Video Lectures on Dynamic Data Assimilation; NPTEL Program IIT Madras: Chennai, India, 2016; Available online: https://nptel.ac.in/courses/111/106/111106082/ (accessed on 28 November 2020).