Analysis Step
ua[:,k+1],B = EKF(ua[:,k+1],w[:,km],h,Dh,R,B)
km = km+1
✝ ✆
```
**Figure 4.** EKF results for the Lorenz 63 system with the assumption of zero process noise.

#### **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*(*n*3), 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 {**u**(*tk*)(*i*)|<sup>1</sup> <sup>≤</sup> *<sup>i</sup>* <sup>≤</sup> *<sup>N</sup>*} 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{aligned} \mathbf{B} &= E[\mathbf{\mathcal{J}}\_b \mathbf{\mathcal{J}}\_b^T] = E[(\mathbf{u}\_t - \mathbf{u}\_b)(\mathbf{u}\_t - \mathbf{u}\_b)^T], \\ \mathbf{P} &= E[\mathbf{\mathcal{J}}\_a \mathbf{\mathcal{J}}\_a^T] = E[(\mathbf{u}\_t - \mathbf{u}\_d)(\mathbf{u}\_t - \mathbf{u}\_d)^T]. \end{aligned}$$

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

$$\begin{aligned} \mathbf{B} & \approx \frac{1}{N-1} \sum\_{i=1}^{N} (\mathbf{u}\_b^{(i)} - \overline{\mathbf{u}\_\mathbf{b}})(\mathbf{u}\_b^{(i)} - \overline{\mathbf{u}\_\mathbf{b}})^T, \\ \mathbf{P} & \approx \frac{1}{N-1} \sum\_{i=1}^{N} (\mathbf{u}\_a^{(i)} - \overline{\mathbf{u}\_\mathbf{a}})(\mathbf{u}\_a^{(i)} - \overline{\mathbf{u}\_\mathbf{a}})^T, \end{aligned}$$

where the bar denotes the ensemble average defined as **<sup>u</sup>** <sup>=</sup> <sup>1</sup> *<sup>N</sup>* <sup>∑</sup>*<sup>N</sup> <sup>i</sup>*=<sup>1</sup> **<sup>u</sup>**(*i*). 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 {**<sup>u</sup>**(*i*) *<sup>b</sup>* (*tk*)|1 ≤ *i* ≤ *N*} at time *tk* drawn from the distribution <sup>N</sup> (**<sup>u</sup>***<sup>b</sup>*(*tk*), **<sup>B</sup>***<sup>k</sup>*), where **<sup>u</sup>***<sup>b</sup>*(*tk*) represents our best-known estimate at *tk*. It can be verified that the ensemble mean and covariance converge to **<sup>u</sup>***<sup>b</sup>*(*tk*) and **<sup>B</sup>***<sup>k</sup>* as *<sup>N</sup>* <sup>→</sup> <sup>∞</sup>. Then, the forecast step is applied to each member of the enesemble as follows,

$$\mathbf{u}\_{b}^{(i)}(t\_{k+1}) = M(\hat{\mathbf{u}}\_{b}^{(i)}(t\_{k}); \boldsymbol{\theta}) + \mathfrak{J}\_{p}^{(i)}(t\_{k+1}),\tag{52}$$

#### *Fluids* **2020**, *5*, 225

where *ξ* (*i*) *<sup>p</sup>* (*tk*<sup>+</sup>1) is drawn from the multivariate Gaussian distribution with zero mean and covariance matrix of **Q***k*+<sup>1</sup> 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,

$$\mathbf{u}\_b(t\_{k+1}) \approx \overline{\mathbf{u}\_\mathbf{b}}(t\_{k+1}) = \frac{1}{N} \sum\_{i=1}^N \mathbf{u}\_b^{(i)}(t\_{k+1}),\tag{53}$$

$$\mathbf{f}\_{b}^{(i)}(t\_{k+1}) = \mathbf{u}\_{b}^{(i)}(t\_{k+1}) - \overline{\mathbf{u}\_{\mathbf{b}}}(t\_{k+1}),\tag{54}$$

$$\mathbf{B}\_{k+1} \approx \frac{1}{N-1} \sum\_{i=1}^{N} \mathfrak{F}\_{b^\*}^{(i)}(t\_{k+1}) \mathfrak{F}\_{b^\*}^{(i)}(t\_{k+1})^T,\tag{55}$$

which provides an approximation for the background covariance at *tk*<sup>+</sup><sup>1</sup> without actually propagating the covariance matrix, as is the case in standard Kalman filtering.

An enesmble of observations {**w**(*i*)(*tk*<sup>+</sup>1)|<sup>1</sup> <sup>≤</sup> *<sup>i</sup>* <sup>≤</sup> *<sup>N</sup>*}, also called virtual observations, is created assuming a Gaussian distribution with a mean equal to the actual observation **w**(*tk*<sup>+</sup>1) and a covariance matrix **R***k*<sup>+</sup>1. In other words, random Gaussian perturbations with zero mean and covariance matrix **R***k*+<sup>1</sup> are added to the actual measurements to create perturbed measurements. The Kalman gain matrix can be computed as before (repeated here for completeness),

$$\mathbf{K}\_{k+1} = \mathbf{B}\_{k+1} \mathbf{D}\_{h} (\mathbf{u}\_{b}(t\_{k+1}))^{T} \left( \mathbf{D}\_{h} (\mathbf{u}\_{b}(t\_{k+1})) \mathbf{B}\_{k+1} \mathbf{D}\_{h} (\mathbf{u}\_{b}(t\_{k+1}))^{T} + \mathbf{R}\_{k+1} \right)^{-1}. \tag{56}$$

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

$$\mathbf{u}\_{a}^{(i)} = \mathbf{u}\_{b}^{(i)}(t\_{k+1}) + \mathbf{K}\_{k+1} \Big( \mathbf{w}^{(i)}(t\_{k+1}) - h(\mathbf{u}\_{b}^{(i)}(t\_{k+1})) \Big),\tag{57}$$

and the analyzed estimate at *tk*<sup>+</sup><sup>1</sup> is computed as the sample mean of the analysis ensemble along with its covariance matrix as follows,

$$\mathbf{u}\_{\mathbf{a}}(t\_{k+1}) \approx \overline{\mathbf{u}\_{\mathbf{a}}}(t\_{k+1}) = \frac{1}{N} \sum\_{i=1}^{N} \mathbf{u}\_{\mathbf{a}}^{(i)}(t\_{k+1}),\tag{58}$$

$$\mathfrak{J}\_{\mathfrak{a}}^{(i)}(t\_{k+1}) = \mathfrak{u}\_{\mathfrak{a}}^{(i)}(t\_{k+1}) - \overline{\mathfrak{u}\_{\mathfrak{a}}}(t\_{k+1}),\tag{59}$$

$$\mathbf{P}\_{k+1} \approx \frac{1}{N-1} \sum\_{i=1}^{N} \mathfrak{F}\_{a}^{(i)}(t\_{k+1}) \mathfrak{F}\_{a}^{(i)}(t\_{k+1})^T. \tag{60}$$

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,

Ensemble initialization: **<sup>u</sup>**(*i*) Virtual observations: **w**(*i*)

Inputs: **<sup>u</sup>***<sup>b</sup>*(*tk*), **<sup>B</sup>***<sup>k</sup>*, *<sup>M</sup>*(·; ·), **<sup>Q</sup>***k*<sup>+</sup>1, **<sup>w</sup>**(*tk*<sup>+</sup>1), **<sup>R</sup>***k*<sup>+</sup>1, *<sup>h</sup>*(·) *<sup>b</sup>* (*tk*) = **<sup>u</sup>***<sup>b</sup>*(*tk*) + **<sup>e</sup>** (*i*) *<sup>b</sup>* , **e** (*i*) *<sup>b</sup>* ∼ N (**0**, **B***<sup>k</sup>*) (*tk*<sup>+</sup>1) = **<sup>w</sup>**(*i*) (*tk*<sup>+</sup>1) + **e** (*i*) *<sup>m</sup>* , **e** (*i*) *<sup>m</sup>* ∼ N (**0**, **R***k*+1) Forecast: **<sup>u</sup>**(*i*) *<sup>b</sup>* (*tk*<sup>+</sup>1) = *<sup>M</sup>*(**<sup>u</sup>**(*i*) *<sup>b</sup>* (*tk*); *θ*) + *ξ* (*i*) *<sup>p</sup>* (*tk*<sup>+</sup>1) *ξ* (*i*) *<sup>b</sup>* (*tk*<sup>+</sup>1) = **<sup>u</sup>**(*i*) *<sup>b</sup>* (*tk*<sup>+</sup>1) − **ub**(*tk*<sup>+</sup>1) **<sup>B</sup>***k*+<sup>1</sup> <sup>≈</sup> <sup>1</sup> *N* − 1 *N* ∑ *i*=1 *ξ* (*i*) *<sup>b</sup>* (*tk*<sup>+</sup>1)*ξ* (*i*) *<sup>b</sup>* (*tk*<sup>+</sup>1)*<sup>T</sup>*

$$\mathbf{B}\_{k+1} \approx \frac{1}{N-1} \sum\_{l=1}^{N} \mathbf{\tilde{z}}\_{b}^{(l)} (t\_{k+1}) \mathbf{\tilde{z}}\_{b}^{(l)} (t\_{k+1})^T$$

$$\text{Kalman gain:} \qquad \mathbf{K}\_{k+1} = \mathbf{B}\_{k+1} \mathbf{D}\_{h} (\mathbf{u}\_{b}(t\_{k+1}))^T \left( \mathbf{D}\_{h} (\mathbf{u}\_{b}(t\_{k+1}))\_{k+1} \mathbf{B}\_{k+1} \mathbf{D}\_{h} (\mathbf{u}\_{b}(t\_{k+1}))^T + \mathbf{R}\_{k+1} \right)^{-1}$$

$$\text{Analysis:} \qquad \mathbf{u}\_{a}^{(l)}(t\_{k+1}) = \mathbf{u}\_{b}^{(l)}(t\_{k+1}) + \mathbf{K}\_{k+1} \left( \mathbf{w}^{(l)}(t\_{k+1}) - h(\mathbf{u}\_{b}^{(l)}(t\_{k+1})) \right)$$

$$\begin{aligned} \mathbf{u}\_{a}(t\_{k+1}) &\approx \overline{\mathbf{u}\_{a}}(t\_{k+1}) \\ \mathbf{\tilde{z}}\_{a}^{(l)}(t\_{k+1}) &= \mathbf{u}\_{a}^{(l)}(t\_{k+1}) - \overline{\mathbf{u}\_{a}}(t\_{k+1}) \\ \mathbf{P}\_{k+1} &\approx \frac{1}{N-1} \sum\_{l=1}^{N} \underline{\mathbf{z}}\_{b}^{(l)}(t\_{k+1}) \mathbf{\tilde{z}}\_{a}^{(l)} (t\_{k+1})^T \end{aligned}$$

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

**Listing 17.** Implementation of the EnKF with virtual observations. ✞ ☎

```
import numpy as np
def EnKF(ubi,w,ObsOp,JObsOp,R,B):