Continuous data assimilation of large eddy simulation by lattice Boltzmann method and local ensemble transform Kalman filter (LBM-LETKF)

The data assimilation (DA) refers to the methodology which interpolates the simulation and the experiment to achieve more consistent state estimation. The continuous DA is one of the major DA fields, in which observation and DA are repeated in time (cf figures 1(a) and (b)). The continuous DA is used to estimate and/or predict the time-dependent state of chaotic systems, which have been the major interest in the community of the numerical weather prediction (NWP) (Kalnay 2002). Recently, the continuous DA, also known as the reconstruction, regeneration, or state estimation of flow fields, has also been applied to the field of computational fluid dynamics (CFD) to predict turbulent flows. Fundamental properties of DA in turbulent flows have been investigated in direct numerical simulations (DNSs) of homogeneous and isotropic turbulence. It was shown that by imposing a reference solution to larger-scale Fourier modes below a certain critical wavenumber, 3D DNSs of isotropic turbulence were fully reproduced or synchronized up to smaller-scale Fourier modes (Yoshida et al 2005, Lalescu et al 2013, Wang et al 2022). The simplest continuous DA approach called the nudging (Hoke and Anthes 1976, Lakshmivarahan and Lewis 2013) was applied to 2D isotropic turbulence, and conditions for synchronization, in which the assimilated solution converges to the reference solution in time, were shown for nudging in Fourier space (Olson and Titi 2003) and in configuration space (Azouani et al 2014), and were examined in 2D DNSs of isotropic turbulence (Olson and Titi 2008, Gesho et al 2016). The DA accuracy of the nudging was investigated also in 3D DNSs of isotropic turbulence, and it was shown that to obtain synchronized solutions, ${\sim} 20\%$ of volume or grid points have to be assimilated in configuration space, while nudging in Fourier space leads to synchronized solutions with an order of magnitude lower number of nudged Fourier modes (Leoni et al 2020). Li et al (2022) also studied the nudging DA of 3D isotropic turbulence with the large eddy simulation (LES).

Zoom In Zoom Out Reset image size

The nudging based on flow field data was also applied to turbulent flows around a square cylinder, where 2D unsteady Reynolds-averaged Navier–Stokes (URANS) simulations were assimilated with the reference flow field data from 3D DNSs (Zauner et al 2022). Although many of these works addressed an impact of observation resolution in configuration space or in Fourier space on the DA accuracy, an influence of observation noise was studied only in a few works. Wang and Zaki (2022) investigated the nudging DA for the turbulent channel flow with observation noises; however, in the latter work, the gain of the nudging was set to 1, i.e. the simulation state was totally replaced by the observation data where they exist. Theoretically, it is difficult to determine the optimal gain of the nudging against given observation/simulation errors, and thus, the gain of the nudging is determined empirically in general.

On the other hand, the issue of observation error was addressed also by using Kalman filter approaches. Suzuki (2012) proposed a hybrid simulation combining particle tracking velocimetry (PTV) and DNS, in which 2D DNSs of a planar-jet flow were assimilated with PTV data from the experiment using the reduced-order Kalman filter. Colburn et al (2011) showed the capability of the ensemble Kalman filter (EnKF) (Evensen 1994, 2003) applied to the flow reconstruction in the turbulent channel flow with the measurement of the pressure on the walls. With increasing computing power, continuous DAs were also applied to engineering problems in recent years. Bauweraerts and Meyers (2021) showed a turbulent flow reconstruction in the atmospheric boundary layer by the 4D-Var DA (Le Dimet and Talagrand 1986) with the (virtual) lidar measurement. Labahn et al (2020) applied the EnKF to 3D LES of a turbulent jet, which was assimilated with flow fields reconstructed based on tomographic particle image velocimetry (TPIV) data from the experiment. They examined its DA accuracy by changing DA conditions such as the number of ensemble members, the localization parameter, the observation error, the DA frequency, and the resolution of reconstructed flow fields.

Although the latter study (Labahn et al 2020) demonstrated promising features of the EnKF in engineering problems, the analysis was still posteriori, and the observation resolution was relatively high (several experimental data points per computational cell). As the goal of DA is to produce accurate forecasts based on observation data, real-time or even faster CFD simulations are needed. Although the space and time resolution of the experimental data is dramatically improving with PIV techniques in the laboratory experiment, it may still be limited in the open fields such as atmospheric boundary layer flows. Therefore, one needs to examine the applicability of the continuous DA to more sparse observation conditions. To address these two issues, in this work, we apply the continuous DA to the lattice Boltzmann method (LBM), and address its DA accuracy with noisy and spatially-sparse observation conditions. Here, the LBM is an essential tool for the large-scale LES, in which over billion degree-of-freedoms (DOFs) can be computed on state-of-the-art supercomputers. In our previous works (Onodera et al 2018, 2021), for instance, the LBM was dramatically accelerated on GPUs, and real-time ensemble LESs for plume dispersion analysis in urban areas were enabled by combining the LBM and the block-structured adaptive mesh refinement (AMR) method. Regarding the continuous DA, we used the nudging in the latter work. However, in the nudging, spatial interpolation is needed when the simulation and the observation have different resolutions. In addition, the nudging requires the conversion of the variables, when the variables of the simulation and the observation are different; namely, the LBM solves the velocity distribution function, while the macroscopic variables such as the velocity and the pressure are observed. Here, the conversion from the macroscopic variables to the velocity distribution cannot be uniquely determined, and an approximation such as assuming the equilibrium distribution function is needed. Such approximations can be naturally avoided by using the EnKF, in which the state vector and the observation vector can be designed by using different variables at different resolutions and dimensions. However, the applicability of the EnKF to the LBM has not been tested to our best knowledge. Thus, in this work, we clarify this point by comparing with our previous method, the nudging.

Here, to consider the candidate of the continuous DA approaches applied to the LBM, we briefly review methodologies developed in the NWP community. The continuous DA can be classified into two major approaches: one is the variational approach and the other is the statistical approach. Historically, in the former, the 3D or 4D variational method (3D/4D-Var) has been the de facto standard DA approach. However, the 3D/4D-Var has difficulty in determining the background error (i.e. the error and the covariance of the numerical state are prescribed by the model hyperparameters). The determination of the background error is often based on an empirical rule, which varies depending on the numerical model. One needs to redesign the background error whenever the model is changed, e.g. adjusting the domain size, the mesh resolution, the time step width, boundary conditions, and/or switching numerical schemes such as space discretization and time integration. Such frequent redesign of the background error is very expensive and is a serious problem in mesh refinement approaches, in which the mesh resolution changes in time. In addition, especially in the 4D-Var, one needs to implement the so-called adjoint code, which is another critical issue for the adaptive mesh refinement approaches.

On the other hand, another major paradigm, which can avoid the redesign of the background error, is the EnKF. The EnKF computes the ensemble simulation, in which the statistics of the numerical model are estimated by multiple simulations with slightly different simulation conditions (cf figure 1(c)). The ensemble simulation is regarded as the Monte Carlo sampling, and thus, its statistical error trivially corresponds to the background error. This is a more straightforward approach than the 3D/4D-Var because an empirical design of the background error is no longer needed. Therefore, the applicability of the EnKF is drastically enhanced from the 3D/4D-Var, and we address the EnKF in this work.

Since the original EnKF (Evensen 1994) had the problem of too small ensemble variation (Burgers et al 1998), the EnKF has been improved in former studies. They are classified mostly into two types: One is the perturbed observation method (EnKF-PO) (Burgers et al 1998), which superimposes the additional ensemble perturbation onto the observation to implicitly maintain the ensemble variation. The other is the ensemble squared root filter (EnSRF) (Tippett et al 2003), which explicitly updates the separation of each ensemble member by solving the squared root of the ensemble covariance matrix. Although the relative merits of the two approaches are still under debate, in this study, we choose one of the EnSRF approaches, the local ensemble transform Kalman filter (LETKF) (Hunt et al 2007). The reason for the decision is that the LETKF can be computed in an embarrassingly parallel manner at each grid point, and thus, is suitable for high-performance computing (HPC). In fact, the LETKF has been applied to extreme-scale DA simulations with large ensemble sizes and/or DOFs on state-of-the-art supercomputers, e.g. on the CPU-based supercomputers in the former works (Miyoshi et al 2014, 2016, Yashiro et al 2020), and on the GPU-based supercomputer in our previous study (Hasegawa et al 2022).

The LETKF and the other EnKFs have been well-studied in the field of CFD, however, a few works of those were conducted with the LBM. Roussel et al (2013) estimated unknown boundary conditions in the 2D lid-driven cavity flow problem by using the EnKF-PO with the velocity measurements. Salman et al (2022) recovered the spatial profile of the temperature in the simulated indoor air quality with unknown boundary temperature by using the 3D-Var or the EnKF-PO with the temperature measurements within the room. Those studies, however, utilized the DA for the estimation of the unknown boundary condition rather than the flow state estimation. Thus, in this study, we apply the LETKF to the LBM, and examine its capability in turbulent flow reconstructions.

As the first test case of the LBM-LETKF, we consider a minimal benchmark problem designed to minimize the requirements on physics modelings and computational costs while keeping chaotic behavior. For example, in the NWP field, the NCEP GFS model (Szunyogh et al 2005) and the SPEEDY model (Molteni 2003) are the minimal 3D global atmospheric models whose DOFs are typically $\mathcal{O}(10^4)$. We also consider a small system with similar DOFs. In this work, we choose a 2D forced isotropic turbulence (Lilly 1969), in which properties of the nudging have been well understood in former works. We also apply an LES model to suppress the grid number to $256\times256$ for computing isotropic turbulence at a relatively high Reynolds number. Under these conditions, the ensemble simulation with up to 64 ensembles can be computed within a few minutes on GPUs, and thus, we can conduct thousands of trial-and-error experiments with reasonable computing resources.

The remainder of this paper is organized as follows. The numerical method of the LBM and the LETKF are described in sections 2 and 3, respectively. While the validation of the LETKF is the major objective of this study, a simpler DA model, the nudging, is also introduced for comparison. In the LBM, the state and observable variables are different, i.e. the former is the velocity distribution function and the latter is the macroscopic variables of the flow field. We also describe how to treat such state and observable variables in the LETKF or the nudging. Section 4 presents the DA experiment of 2D forced isotropic turbulence with the LBM-LETKF. We first design and validate the numerical configuration, and then, the DA experiments using the LBM-LETKF and the LBM-nudging are compared. As the observation data, we assume spatially-sparse and noisy observation data of velocity and pressure field, which may be considered for real measurements of the atmospheric boundary layer environment. Especially, data sparsity will be studied in detail to investigate the accuracy and robustness of the LETKF and to show the advantage of the LETKF compared to the simple nudging method. Finally, we conclude this study in section 5.

We start with the definition of the CFD model without the DA scheme. Since our future target is to implement the ensemble DA into the LBM-based real-time urban wind simulation code (Onodera et al 2021), we employ the LBM as the CFD model for this study. Although our previous study utilized the 3D cumulant model (Geier et al 2015), in this study we use the traditional 2D lattice BGK model (Qian et al 1992) for simplicity. In the LBM with D2Q9 (two-dimensional nine-velocities) model, the discrete velocity $\boldsymbol{c}_\alpha$ is defined as

$$\begin{align} \boldsymbol{c}_\alpha = \left\{\begin{array}{ll} \left(0, 0\right)^\top, & \alpha = 0 \\ \left(\pm c, 0\right)^\top, & \alpha = 1, 2 \\ \left(0, \pm c\right)^\top, & \alpha = 3, 4 \\ \left(\pm c, \pm c\right)^\top, & \alpha = 5, 6, 7, 8 \end{array}\right., \end{align} \tag{ 1 }$$

where the subscript α is the index of the discrete velocity. c is the model constant called the lattice speed, which is determined from the given grid spacing $\Delta x$ and time step width $\Delta t$, namely, $c = \Delta x / \Delta t$. The governing equation is defined as

$$\begin{align} f_\alpha\left(\boldsymbol{\xi}+\boldsymbol{c}_\alpha\Delta t, t+\Delta t\right) - f_\alpha\left(\boldsymbol{\xi}, t\right) = \Delta t\left[ - \frac{1}{\tau}\left( \,f_\alpha\left(\boldsymbol{\xi},t\right) - f\,\,_\alpha^\mathrm{eq}\left(\boldsymbol{\xi},t\right) \right) + F_\alpha \right], \end{align} \tag{ 2 }$$

where $\boldsymbol{\xi} = (x,y)$ and t are the coordinate and the time, respectively. f_α is the velocity distribution function, whose moments correspond to the macroscopic variables as

$$\begin{align} \rho & = \sum_\alpha f_\alpha, \end{align} \tag{ 3 }$$

$$\begin{align} \boldsymbol{u} & = \frac{1}{\rho} \sum_\alpha f_\alpha \boldsymbol{c}_\alpha, \end{align} \tag{ 4 }$$

where ρ is the fluid density, and $\boldsymbol{u} = (u, v)$ is the fluid velocity. $f\,\,^\mathrm{eq}_\alpha$ is the equilibrium distribution function, defined by

$$\begin{align} f\,\,^\mathrm{eq}_\alpha = \rho w_\alpha \left( 1 + 3\frac{\boldsymbol{c}_\alpha \cdot \boldsymbol{u}}{c^2} + \frac{9}{2}\frac{\left(\boldsymbol{c}_\alpha \cdot \boldsymbol{u}\right)^2}{c^4} - \frac{3}{2}\frac{\boldsymbol{u}\cdot\boldsymbol{u}}{c^2} \right), \end{align} \tag{ 5 }$$

where w_α is the weighting factor, which is given by

$$\begin{align} w_\alpha = \left\{ \begin{array}{ll} 4/9, & \alpha = 0 \\ 1/9, & \alpha = 1,2,3,4 \\ 1/36, & \alpha = 5,6,7,8 \end{array}\right.. \end{align} \tag{ 6 }$$

τ is the relaxation time. From Chapman-Enskog expansion (Krüger et al 2017), given the lattice speed c and the time step width $\Delta t$, the relaxation time is chosen to satisfy

$$\begin{align} \nu + \nu_\mathrm{SGS} = \frac{c^2}{3}\left(\tau - \frac{\Delta t}{2}\right), \end{align} \tag{ 7 }$$

where ν is the kinetic viscosity and $\nu_\mathrm{SGS}$ is the sub-grid scale (SGS) viscosity. The SGS viscosity is computed based on the standard Smagorinsky model (Smagorinsky 1963) because the benchmark problem in this study is the isotropic turbulence and thus the velocity field is statistically homogeneous. $\nu_\mathrm{SGS}$ is then defined as

$$\begin{align} \nu_\mathrm{SGS} & = \left(C_s\bar{\Delta}\right)^2|S|, \end{align} \tag{ 8 }$$

$$\begin{align} \notag\\|S| & = \sqrt{2 S_{ij} S_{ji}}, \end{align} \tag{ 9 }$$

$$\begin{align} \notag\\ S_{ij} & = \frac{1}{2}\left( \frac{\partial u_i}{\partial \xi_j} + \frac{\partial u_j}{\partial \xi_i} \right). \end{align} \tag{ 10 }$$

Here, $C_s = 0.2$ is the Smagorinsky constant, and $\bar{\Delta}$ is the model constant called the filter length, which is set to $\bar{\Delta} = \Delta x$ in this study. S_ij is the strain-rate tensor with the subscripts $i,j$ being the indices of spatial axes. F_α is the forcing term, which is defined with the given acceleration F by

$$\begin{align} F_\alpha = 3w_\alpha \rho \frac{\boldsymbol{F}\cdot \boldsymbol{c}_\alpha}{c^2}. \end{align} \tag{ 11 }$$

The continuous data assimilation (DA) corrects model variables in the numerical simulation based on the observation data iteratively in time. We introduce several vectors to treat the numerical and observation data in the context of the DA: let $\boldsymbol{x}(\boldsymbol{\xi})$ and $\boldsymbol{y}(\boldsymbol{\eta})$ be respectively the vector of numerical state (state vector) at the grid space ξ , and the vector of observable variables (observation vector) at the grid space η . The number of grid points in the observation space η is usually much smaller than that in the simulation state ξ . In this study, the state vector and observation vector are respectively chosen as the distribution function of the LBM and the macroscopic variables, namely,

$$\begin{align} \boldsymbol{x}\left(\boldsymbol{\xi}\right) = \left[\,f_0\left(\boldsymbol{\xi}\right), f_1\left(\boldsymbol{\xi}\right), \dots, f_8\left(\boldsymbol{\xi}\right)\right]^\top, \end{align} \tag{ 12 }$$

$$\begin{align} \boldsymbol{y}\left(\boldsymbol{\eta}\right) = \left[\rho\left(\boldsymbol{\eta}\right), u\left(\boldsymbol{\eta}\right), v\left(\boldsymbol{\eta}\right)\right]^\top. \end{align} \tag{ 13 }$$

Let $\boldsymbol{x}^\mathrm{f}$, $\boldsymbol{x}^\mathrm{a}$, and $\boldsymbol{y}^\mathrm{o}$ be respectively the forecast (prior) and the analysis (posterior) state vectors in the numerical model, and the observation vector (measurement data) from the true state at the time t. Hereafter, the superscripts 'f', 'a', and 'o' are used also for other variables in the same convention. The workflow of the continuous DA, which repeats the cycle of the time integration of the numerical model and the DA, consists of the following procedures:

• (i)  
  (Initial condition) Set the initial analysis state vector, $\boldsymbol{x}^\mathrm{a}(\boldsymbol{\xi})|_{t = 0}$. Set t = 0.
• (ii)  
  (Forecast) Compute the time integration of the numerical model and get the forecast state vector in the next time step, $\boldsymbol{x}^\mathrm{f}(\boldsymbol{\xi})|_{t = t+s\Delta t}$. Here, the number of time steps per the observation period is chosen as s = 200 in this study. The reason for this choice is based on the maximum Lyapunov exponent of the turbulence, which will be discussed in section 4.2.
• (iii)  
  (Observation) Measure the ground truth state and get the observation vector, $\boldsymbol{y}^\mathrm{o}(\boldsymbol{\eta})|_{t = t+s\Delta t}$.
• (iv)  
  (Analysis) Assimilate the forecast state vector with the observation vector and get the new analysis state vector based on the linear estimation theory (Talagrand 1997, Desroziers et al 2005):

  $$\begin{align} \boldsymbol{x}^\mathrm{a} = \boldsymbol{x}^\mathrm{f} - \boldsymbol{K}\left(\boldsymbol{y}^\mathrm{o} - \mathcal{H}\boldsymbol{x}^\mathrm{f}\right), \end{align} \tag{ 14 }$$

  K is the gain matrix, which varies depending on the DA method. $\mathcal{H}$ is the observation operator defined as $\boldsymbol{y} = \mathcal{H}\boldsymbol{x}$.
• (v)  
  Update the time step, $t \leftarrow t+s\Delta t$, and repeat the procedures (i)–(iv).

In practice, once the vectors x and y are designed, the observation operator $\mathcal{H}$ is designed to satisfy $\boldsymbol{y} = \mathcal{H}\boldsymbol{x}$, and the gain K is finally given. The design of $\boldsymbol{x}, \boldsymbol{y}$, and K in each DA method is described in the following subsections.

3.1. Nudging

Firstly, we introduce a traditional DA approach, which is called nudging (Hoke and Anthes 1976, Lakshmivarahan and Lewis 2013), as a reference in this work. The nudging gives the simplest linear interpolation between the state and the observation, i.e.

$$\begin{align} \boldsymbol{K} = \gamma\mathcal{H}^{-1}, \end{align} \tag{ 15 }$$

where $0\unicode{x2A7D}\gamma\unicode{x2A7D}1$ is an empirical scalar constant. In this study, we choose optimal γ in each case by the preliminary numerical experiment. Since the nudging method is an elementwise DA procedure, the state and the observation can be treated as scalars rather than vectors. $\mathcal{H}^{-1}$ is the inversion of the observation operator, namely, $\boldsymbol{x} = \mathcal{H}^{-1}\boldsymbol{y}$. To compute $\mathcal{H}^{-1}\boldsymbol{y}$, one needs a priori data processing which transforms from the macroscopic variables to the velocity distribution function. However, such transformation is not unique. In this study, we assume that the observation is given as the equilibrium distribution, and the non-equilibrium components are ignored.

In addition, one also needs spatial interpolation for the observation data if the resolution of the observation is coarser than that of the simulation. In this study, we apply simple linear interpolation to the sparse observation data. Although these approximations may ruin the DA accuracy, the nudging has merits in its computing cost and ease of use. Finally, the nudging for the LBM in this study is formed as

$$\begin{align} f\,^\mathrm{a}_\alpha\left(\boldsymbol{\xi}\right) = f\,^\mathrm{f}_\alpha\left(\boldsymbol{\xi}\right) + \gamma \left[f\,\,^\mathrm{eq}_\alpha|_{\rho = \hat{\rho}^\mathrm{o}\left(\boldsymbol{\xi}\right), u = \hat{u}^\mathrm{o}\left(\boldsymbol{\xi}\right), v = \hat{v}^\mathrm{o}\left(\boldsymbol{\xi}\right)} - f\,\,^\mathrm{f}_\alpha\left(\boldsymbol{\xi}\right)\right], \end{align} \tag{ 16 }$$

where $\hat{\rho}^\mathrm{o}(\boldsymbol{\xi})$, $\hat{u}^\mathrm{o}(\boldsymbol{\xi})$, and $\hat{v}^\mathrm{o}(\boldsymbol{\xi})$ are spatially interpolated variables of the observation $\rho^\mathrm{o}(\boldsymbol{\eta})$, $u^\mathrm{o}(\boldsymbol{\eta})$, and $v^\mathrm{o}(\boldsymbol{\eta})$, respectively.

3.2. Local ensemble transform Kalman filter (LETKF)

The LETKF (Hunt et al 2007) is one of the advanced models of the continuous DA, which is a derived form of the EnKF (Evensen 2003). In the LETKF, the DA is defined along with the ensemble simulation. One of its advantages is that the gain K is computed from the statistics of the ensemble simulation, and the adjoint of observation operator $\mathcal{H}^\dagger$ is thus not required. It is noted that the observation operator $\mathcal{H}$ and the gain K are not explicitly given for simplicity. To see their step-by-step deviations, readers can refer to, e.g. the book by Li et al (2012).

3.2.1. Basic model.

Let M be the number of ensemble members (i.e. the ensemble size), and $\boldsymbol{x}_m^\mathrm{f}$ be the forecast state vector of the mth ensemble member, where $m = 0,1,\dots, M-1$. Their ensemble mean and perturbation are then introduced as

$$\begin{align} \bar{\boldsymbol{x}}^\mathrm{f} & = \frac{1}{M}\sum_{m = 0}^{M-1} \boldsymbol{x}_m^\mathrm{f}, \end{align} \tag{ 17 }$$

$$\begin{align} \delta \boldsymbol{X}^\mathrm{f} & = \left[ \boldsymbol{x}_0^\mathrm{f} - \bar{\boldsymbol{x}}^\mathrm{f}, \boldsymbol{x}_1^\mathrm{f} - \bar{\boldsymbol{x}}^\mathrm{f}, \dots, \boldsymbol{x}_{M-1}^\mathrm{f} - \bar{\boldsymbol{x}}^\mathrm{f} \right]. \end{align} \tag{ 18 }$$

Similarly, those mapped to the observation space are defined by

$$\begin{align} \bar{\boldsymbol{y}}^\mathrm{f} & = \frac{1}{M}\sum_{m = 0}^{M-1} \boldsymbol{y}_m^\mathrm{f}, \end{align} \tag{ 19 }$$

$$\begin{align} \delta \boldsymbol{Y}^\mathrm{f} & = \left[ \boldsymbol{y}_0^\mathrm{f} - \bar{\boldsymbol{y}}^\mathrm{f}, \boldsymbol{y}_1^\mathrm{f} - \bar{\boldsymbol{y}}^\mathrm{f}, \dots, \boldsymbol{y}_{M-1}^\mathrm{f} - \bar{\boldsymbol{y}}^\mathrm{f} \right], \end{align} \tag{ 20 }$$

where $\boldsymbol{y}^\mathrm{f}_m = {\cal H}\boldsymbol{x}^\mathrm{f}_m$.

The analysis state vectors are solved in a latent ensemble space. Let $\boldsymbol{w}^\mathrm{a}_m\in\mathbb{R}^M$ be the analysis state vector in the ensemble space, which is given by the transform $\boldsymbol{x}^\mathrm{a} = \boldsymbol{X}^\mathrm{f}\boldsymbol{w}^\mathrm{a}$. Its ensemble mean $\bar{\boldsymbol{w}}^\mathrm{a}$ and perturbation matrix $\delta \boldsymbol{W}^\mathrm{a}$ are then computed by

$$\begin{align} \bar{\boldsymbol{w}}^\mathrm{a} & = \tilde{\boldsymbol{P}}^\mathrm{a} \, \delta \boldsymbol{Y}^\mathrm{f\top} \boldsymbol{R}^{-1} \left(\boldsymbol{y}^\mathrm{o} - \bar{\boldsymbol{y}}^\mathrm{f}\right), \end{align} \tag{ 21 }$$

$$\begin{align} \delta \boldsymbol{W}\,^\mathrm{a} & = \left[\left(M-1\right) \tilde{\boldsymbol{P}}^\mathrm{a}\right]^{1/2}. \end{align} \tag{ 22 }$$

Here, R is the matrix called the observation error covariance, which is a model constant determined by the observation noise. $\tilde{\boldsymbol{P}}^\mathrm{a}$ is the analysis covariance in the ensemble space, given by

$$\begin{align} \tilde{\boldsymbol{P}}^\mathrm{a} = \left[ \left(M-1\right)\boldsymbol{I} + \delta \boldsymbol{Y}^{\mathrm{f}\top} \boldsymbol{R}^{-1} \delta\boldsymbol{Y}^{\mathrm{f}} \right]^{-1}. \end{align} \tag{ 23 }$$

Finally, the analysis state vector in the ensemble space is transformed back to the model space, so that the analysis state vector of the mth ensemble member $\boldsymbol{x}_m^{\mathrm{a}}$ is computed as

$$\begin{align} \boldsymbol{x}_m^{\mathrm{a}} = \bar{\boldsymbol{x}}^\mathrm{f} + \delta \boldsymbol{X}^\mathrm{f} \left( \bar{\boldsymbol{w}}^\mathrm{a} + \delta \boldsymbol{W}\,^{\mathrm{a}}|_m \right), \end{align} \tag{ 24 }$$

where $\delta \boldsymbol{W}^{\mathrm{a}}|_m$ denotes the mth column of $\delta \boldsymbol{W}\,^\mathrm{a}$.

3.2.2. Covariance localization and inflation.

In most practical applications, the EnKF faces a lack of ensemble size. The ensemble size and the DOFs of the system are typically $M = \mathcal{O}(10)\sim \mathcal{O}(10^3)$ and $N = \mathcal{O}(10^4)\sim \mathcal{O}(10^{10})$, respectively. In such EnKF, spatial localization of the covariance (covariance localization) is generally considered to reduce the pseudo-correlations coming from the sampling error due to the small ensemble size. In this study, we utilize the covariance localization in the observation space based on the R -localization (Hunt et al 2007, Miyoshi 2011). The R -localization increases the observation error at the spatially distant observation points so that the influences from those points are weakened. The localized covariance is given by $\boldsymbol{R}_\mathrm{loc} = \boldsymbol{G}^{-1} \circ \boldsymbol{R}$. Here, G is the matrix of localization weight, whose elements have values between 0 to 1, and the values generally decay depending on the spatial distance between the grid point, ξ , and the observation points, namely, η . The symbol $\circ$ denotes the elementwise matrix–matrix multiplication (i.e. the Schur product). In the case that the observation error of each variable is uncorrelated, R and G are diagonal, and the above equation can be rewritten as $\boldsymbol{R}_\mathrm{loc}^{-1} = \boldsymbol{G} \circ \boldsymbol{R}^{-1}$. Here, we employ the Gaspari–Cohn function as the definition of each element of G . This function gives a nearly Gaussian decay, $e^{-1.5d^2}$, approximated by

$$\begin{align} G\left(d\right) & = \left\{\begin{array}{ll} -\frac{1}{4}d^5 + \frac{1}{2}d^4 + \frac{5}{8}d^3 - \frac{5}{3}d^2 + 1, & 0\unicode{x2A7D} d \unicode{x2A7D}1 \\ \frac{1}{12}d^5 - \frac{1}{2}d^4 + \frac{5}{8}d^3 + \frac{5}{3}d^2 - 5d + 4 - \frac{2}{3}d^{-1}, & 1\unicode{x2A7D} d\unicode{x2A7D}2 \\ 0, & d\unicode{x2A7E} 2, \end{array}\right., \end{align} \tag{ 25 }$$

where $d = |\boldsymbol{\xi} - \boldsymbol{\eta}|/d_\mathrm{cut}$ and $d_\mathrm{cut}$ is the cutoff parameter. The merit to use this function is that the distant observation points can be explicitly ignored because the approximated value is exactly zero for $d\unicode{x2A7E}2$ (see figure 5). By using R -localization, the observation error covariance $\boldsymbol{R}_\mathrm{loc}$ and thus the analysis vectors $\bar{\boldsymbol{w}}^\mathrm{a}, \delta \boldsymbol{W}^\mathrm{a}|_m$ have different values on each grid point. It leads to the batched computation of matrices, where each batch computes the analysis at each grid point. Such batched computation was parallelized and accelerated by GPU in our previous study (Hasegawa et al 2022).

It is also commonly known that the EnKF faces the so-called filter divergence problem, in which too small error covariances result in spuriously tiny Kalman gains. As a workaround for the filter divergence, we employ the multiplicative covariance inflation method (Anderson and Anderson 1999), namely, $\boldsymbol{P}^\mathrm{f}_\mathrm{inf} = \beta \boldsymbol{P}^\mathrm{f}$, where β is the inflation coefficient. The value of β is, practically, hand-tuned ad hoc (Hunt et al 2007, Miyoshi et al 2007) or automatically estimated by the adaptive inflation algorithm (Anderson 2007, Miyoshi 2011). In this study, the hand-tuned constant parameters were used. It is noted that the hand-tuned inflation coefficients tended to have larger values with smaller ensemble sizes or with fewer observation points.

In the LETKF, the R -localization and the multiplicative covariance inflation affect only $\bar{\boldsymbol{w}}^\mathrm{a}$ and $\tilde{\boldsymbol{P}}^\mathrm{a}$, and thus, equations (21) and (23) are rewritten as

$$\begin{align} \bar{\boldsymbol{w}}^\mathrm{a} & = \tilde{\boldsymbol{P}}^\mathrm{a} \, \delta \boldsymbol{Y}^\mathrm{f\top} \boldsymbol{R}_\mathrm{loc}^{-1} \left( \boldsymbol{y}^\mathrm{o} - \bar{\boldsymbol{y}}^\mathrm{f} \right), \end{align} \tag{ 26 }$$

$$\begin{align} \tilde{\boldsymbol{P}}^\mathrm{a} & = \left( \frac{M-1}{\beta}\boldsymbol{I} + \delta \boldsymbol{Y}^{\mathrm{f}\top} \boldsymbol{R}_\mathrm{loc}^{-1} \delta\boldsymbol{Y}^{\mathrm{f}} \right)^{-1}. \end{align} \tag{ 27 }$$

4.1. Design of numerical configuration

To clarify the properties of the ensemble DA, we need a numerical experiment, which shows chaotic dynamics and requires the minimum computational cost. The latter feature is essential to study the properties of the LETKF with various observation conditions up to large ensemble sizes. For this purpose, we employ a minimal 2D problem of the forced isotropic turbulence, which is well understood by the classical theory by Kraichnan (1967), Leith (1968), Batchelor (1969) (KLB theory).

We consider the 2D squared computational domain on the x − y plane with the system size of $2\pi\times2\pi$, double periodic boundaries, and the grid number of $256\times256$. The acceleration is given by the superposition of the energy injection and the friction, $\boldsymbol{F} = \boldsymbol{F}^i + \boldsymbol{F}^\mu$. Here, the energy injection F ⁱ is given by

$$\begin{align} \boldsymbol{F}^i & = \Delta x\nabla F_e\times \boldsymbol{i}_z, \end{align} \tag{ 28 }$$

$$\begin{align} F_e & = F_{e,0} \left( \sum_{||\boldsymbol{k}| - k_\mathrm{f}| \unicode{x2A7D} \Delta k} \boldsymbol{r}_1 \cdot \boldsymbol{k} \cos{\left(\boldsymbol{k} \cdot \boldsymbol{\xi}\right)} + \boldsymbol{r}_2 \cdot \boldsymbol{k} \sin{\left(\boldsymbol{k} \cdot \boldsymbol{\xi}\right)} \right), \end{align} \tag{ 29 }$$

where k is the wavenumber vector, i _z is the unit vector in the direction of z-axis, $k_\mathrm{f}$ and $\Delta k$ are respectively the center and width of the energy injection wavenumber regime, r ₁ and r ₂ are Gaussian random vectors with zero-mean and the standard deviation being $\Delta x$, and $F_{e,0}$ is the energy injection parameter. $F_{e,0}$ is determined so that the velocity becomes of the order of the reference velocity on average. The friction is defined as

$$\begin{align} \boldsymbol{F}^\mu = -\mu\boldsymbol{u}, \end{align} \tag{ 30 }$$

with the friction parameter µ. Here, the motivation to introduce the friction force is to avoid the immoderate energy increase in the low-wavenumber regime, which is known as the energy condensation phenomenon (Smith and Yakhot 1994, Chertkov et al 2007). This is essential for quasi-steady long-term simulations. The details of other parameters are summarized in table 1, where the normalization is chosen so that the reference velocity u₀ and the reference density ρ₀ are unity. The time-stepping width $\Delta t$ is adjusted so that the Courant number to the reference velocity is set to 0.05. In these parameters, the SGS viscosity $\nu_\mathrm{SGS}$ has the value of the same order as the kinetic viscosity ν.

Table 1. Detailed configuration of 2D forced isotropic turbulence.

Symbol	Description	Value
N	grid size	2562
u0	reference velocity (Courant number)	1
$\Delta t$	time-stepping width	0.001 232
ρ0	reference fluid density	1
ν	kinetic viscosity	10−4
µ	friction force coefficient	$5\times10^{-4}$
$k_\mathrm{f}$	center of energy injection wavenumber regime	4
$\Delta k$	width of energy injection wavenumber regime	2

Figure 2 shows a typical quasi-stationary state of the computed 2D turbulence. Here, the energy spectrum is defined by the shell average over the 2D wavenumber:

$$\begin{align} E\left(k\right) & = \sum_{||\boldsymbol{k}^{^{\prime}}| - k|\unicode{x2A7D}1/2} E\left(\boldsymbol{k}^{^{\prime}}\right), \end{align} \tag{ 31 }$$

$$\begin{align} E\left(\boldsymbol{k}\right) & = \frac{1}{2}\ \left|\hat{\boldsymbol{u}}\left(\boldsymbol{k}\right)\right|^2, \end{align} \tag{ 32 }$$

where k is the 2D wavenumber, k is the 1D wavenumber for the shell average, $\hat{\boldsymbol{u}}(\boldsymbol{k})$ is the Fourier transform of u . There exists the energy spectrum with the dual cascade, composed of the inverse cascade for $k \lt k_\mathrm{f}$, and the normal cascade for $k_\mathrm{f} \lt k \lt k_\mathrm{d}$. The dissipation regime is also observed in $k \gt k_\mathrm{d}$, where $k_\mathrm{d} \approx 90$. Although the dual cascade was observed, the normal cascade does not follow the power law of k⁻³ and is characterized by the steeper slope ${\approx} k^{-4}$. Although it seems to be inconsistent with the KLB theory, it is a well-known nature in isotropic turbulence simulations. The KLB theory assumes the limit of infinite resolution and/or Reynolds number, which cannot be satisfied in simulations with finite resolution, and thus, viscosity. Theoretically, the power law for the normal cascade is modified to the range between k⁻³ and k⁻⁵ for the moderate resolution and Reynolds number (Tran and Shepherd 2002, Tran and Bowman 2003). The power lows of the normal cascade in former numerical studies were also steeper than k⁻³ (Herring and Mcwilliams 1985, Legras et al 1988, Scott 2007, Tsang and Young 2009, Clark et al 2020, Xie 2020). It is noted that simulations at higher resolution may make the result close to the KLB theory (see the work by Boffetta (2007): the mesh resolution up to $32\,768^2$ resulted in the power low of ${\approx} k^{-3.35}$); however, such a larger scale simulation is not in the scope of the current study, i.e. the minimal benchmark of the ensemble DA with the LBM-LETKF. In addition, the LBM solves weakly-compressible flows rather than incompressible flows, and thus, flow characteristics may be slightly differ from the KLB theory. For these reasons, in this study, we accept the normal cascade with the power law steeper than k⁻³ and adopt this numerical setup as the minimal benchmark problem.

Zoom In Zoom Out Reset image size

4.2. Lyapunov exponent of the turbulence

Firstly, we investigated the error growth (also known as the uncertainty of the initial condition) at the designed numerical setup. The purpose is to estimate the maximum Lyapunov exponent, which indicates the predictability of turbulence (Nastac et al 2017) and gives guidance to design the time interval of the continuous DA (Labahn et al 2020). In this test, we prepared the following perturbed distribution function $f\,^\mathrm{pert}_\alpha$ at t = 0,

$$\begin{align} f\,_\alpha^\mathrm{pert} & = f_{\alpha,0} + \delta_0 \left( f\,\,_\alpha^\mathrm{eq}\left(\rho_0,u^\mathrm{rand},v^\mathrm{rand}\right) - f_{\alpha,0}\right), \end{align} \tag{ 33 }$$

where $u^\mathrm{rand}, v^\mathrm{rand}$ are given by Gaussian random numbers with zero-mean and the standard deviation being the reference velocity u₀. $f_\alpha^\mathrm{eq}(\rho_0,u^\mathrm{rand},v^\mathrm{rand})$ is the equilibrium distribution function with the constant density ρ₀ and the random velocity ($u^\mathrm{rand}, v^\mathrm{rand}$), and δ₀ is the perturbation magnitude. Here, the time t = 0 is not the initial step of the simulation; the original distribution function $f_{\alpha,0}$ is prepared by the spin-up simulation for $10^6\Delta t$, during which the system is fully-developed and is converged to the quasi-stationary state.

The separation of the simulations with and without the perturbation was evaluated by

$$\begin{align} |\delta\left(t\right)| = \sqrt{ \frac{1}{N} \sum_{\boldsymbol{\xi}}\left|\boldsymbol{u}^\mathrm{pert}\left(\boldsymbol{\xi}, t\right) - \boldsymbol{u}\left(\boldsymbol{\xi}, t\right)\right|^2 }, \end{align} \tag{ 34 }$$

where $\boldsymbol{u}^\mathrm{pert}$ and u are the velocities in perturbed and unperturbed simulations, respectively.

Following the chaotic nature of turbulence, the error grows exponentially; then, $|\delta(t)|$ may be approximated by

$$\begin{align} |\delta\left(t\right)| = e^{\lambda t} |\delta\left(0\right)|, \end{align} \tag{ 35 }$$

where λ is the error growth rate, namely, the maximum Lyapunov exponent.

We evaluated the maximum Lyapunov exponent against the various ranges of the initial perturbations from 10⁻⁵ to 10⁻². In figure 3, the time history of the error and the vorticity contour map indicate the growth of the error. As seen in the vorticity contour map, the initial states are quite similar since the initial perturbation is small enough compared to the turbulence; in contrast, the final states are completely different between perturbed and unperturbed simulations, which reflect the error growth following the chaotic dynamics of turbulence. The error grows exponentially and saturates at the level of $\delta(t)\sim\mathcal{O}(1)$ as seen in the time history. The maximum Lyapunov exponents, corresponding to the slopes of the time histories, were almost the same regardless of the magnitude of the initial perturbation. The value of the maximum Lyapunov exponent was $\lambda\sim4.64\times10^{-4}\Delta t^{-1}$, and then, the predictability time was estimated as $t_\mathrm{pred} = \lambda^{-1} \sim 2150\Delta t$.

Zoom In Zoom Out Reset image size

The time interval of the continuous DA for chaotic systems has to be frequent enough compared to the predictability time (Labahn et al 2020). Hence, we determined the design value of the sub-timesteps per DA cycle as s = 200, which refers to $s\Delta t \lessapprox t_\mathrm{pred}/10$.

4.3. Data assimilation experiments

As shown in figure 4, we carried out the DA experiment in the observing system simulation experiment (OSSE) manner. First, the ground truth was computed with a certain initial condition. To emulate the noisy observation, the macroscopic observable variables $(u^\mathrm{o}, v^\mathrm{o}, \rho^\mathrm{o})$ were measured with artificial noises on equidistant sampling points in the ground truth. Here, the density $\rho^\mathrm{o}$ was also measured because the density in the LBM reflects the pressure ($P = \rho/3$), which would be also observable in real experiments. The artificial observation noises were given by Gaussian random numbers with zero-mean, whose standard deviations for $u^\mathrm{o}$/$v^\mathrm{o}$ and $\rho^\mathrm{o}$ were 10% of u₀ and 1% of ρ₀, respectively. Here, the noise amplitude for the density differs from that for the velocity because the density corresponds to the pressure. The noise amplitude for the density is designed based on its fluctuation ${\lt}{\pm}5\%$. The spatiotemporal resolution of the observation was designed as follows:

• The observation points were sampled uniformly with the spatial interval p, which gave $256/p\times256/p$ observation points. The observation condition was varied from the dense observation case, p = 1, to the sparse observation cases, $p = 2, 4, \dots, 64$.
• The temporal interval of the observation was chosen to be the same as the interval of the continuous DA, s = 200.

Zoom In Zoom Out Reset image size

Next, the other simulations were prepared with the random initial conditions. The initial conditions for the ground truth and the ensemble simulations are all uncorrelated with each other. We conducted the simulations without the DA, with the nudging, and with the LETKF. In the nudging cases, the gain parameter was tuned by the preliminary numerical experiments. The optimal parameters were γ = 0.11 in the dense observation case (section 4.3.1) and $0.11\unicode{x2A7D}\gamma\unicode{x2A7D}0.31$ in the sparse observation cases (section 4.3.2). Here, more sparse observation leads to larger γ. In the LETKF cases, the ensemble sizes were chosen as $M = 4, 16, 64$ and the inflation parameter was optimized in each case. In the R -localization, the cutoff parameter was chosen to be $d_\mathrm{cut} = p\Delta x$, which guarantees that each computational grid is affected by ${\sim} {\cal O}(10)$ observation points within the distance of $2d_\mathrm{cut}$ (see figure 5).

Zoom In Zoom Out Reset image size

4.3.1. Accuracy evaluation at the dense observation case.

Firstly, we investigated the DA experiment with the dense observation case, p = 1. The DA accuracy was defined by the root means squared error (RMSE), namely,

$$\begin{align} |\delta\left(t\right)| = \sqrt{\frac{1}{N}\sum_{\boldsymbol{\xi}} |{\boldsymbol{u}}^\mathrm{DA}\left(\boldsymbol{\xi}, t\right) - \boldsymbol{u}^\mathrm{t}\left(\boldsymbol{\xi}, t\right)|^2}, \end{align} \tag{ 36 }$$

where $\boldsymbol{u}^\mathrm{t}$ is the velocity in the ground truth state, and $\boldsymbol{u}^\mathrm{DA}$ is the velocity obtained by each DA approach. In the LETKF, the ensemble mean, ${\boldsymbol{u}}^\mathrm{DA} = M^{-1}\sum_{m = 0}^{M-1} \boldsymbol{u}_m$, is considered as the best estimate. Figure 6 shows the time history of the RMSE for each DA approach. The errors in all cases were $\mathcal{O}(1)$ at t = 0 and became smaller in time except for the non-DA case. They converged into the values of 0.0190 (nudging), 0.0121 (LETKF, M = 4), 0.0110 (LETKF, M = 16), and 0.0109 (LETKF, M = 64), while the RMSE of the observation noise was 0.142. It is noted that the observation noise was 0.1 for each of $u^\mathrm{o}, v^\mathrm{o}$, thus, the observation noise of the velocity magnitude was $0.1\times\sqrt2 = 0.142$. The converged value is calculated by the time-averaging over the converged state for $70\,000\unicode{x2A7D} t/\Delta t\unicode{x2A7D}80\,000$. Although both the nudging and LETKF cases achieved good results, the LETKF cases showed lower RMSEs than the nudging case. The LETKF gave better results at larger ensemble sizes, and its accuracy was almost converged to the best value for $M\unicode{x2A7E} 16$. Besides, the duration to converge the RMSE was shorter in the LETKF cases than in the nudging cases. In the LETKF cases, the RMSE reached around the lower bound at $t/\Delta t = 10\,000\sim20\,000$, while the nudging case required ${\sim}70\,000\Delta t$ to reach the converged state. This feature is due to the difference in the gain: In the nudging, the gain parameter is a constant value tuned for the fully assimilated state ($t/\Delta t\gtrapprox 70\,000$); thus, the gain parameter might be too small for the initial transient period ($t/\Delta t \lt 70\,000$), where the calculation error was larger than the observation error. In contrast, in the LETKF, the gain is dynamically determined by the error covariances of the state vector and the observation vector, which result in the larger gain at the larger error covariance during the initial transient period. This feature of the LETKF would be convenient to reduce the spin-up time from the non-assimilated state to the fully assimilated state. Here, it should be noted that the faster convergence may be achieved also in the nudging if one introduces an adaptive optimization of γ, so that it becomes larger in the initial transient period. However, it is hard to optimize γ in an adaptive manner, because in reality, we cannot estimate the RMSE against the true state and the nudging does not have automatic optimization strategy to fit the model to the true state. On the other hand, the LETKF can automatically optimize the Kalman gain based on the ensemble statistics. Therefore, the adaptive nudging might require large number of test runs to try different time-varying parameters by hand-tuning. In addition, such a hand-tuned parameter might not be optimal for other initial conditions. Thus, we did not consider to introduce the adaptive nudging in this study.

Zoom In Zoom Out Reset image size

To see the details of the fully assimilated state in the nudging and LETKF cases, we examined the energy spectra. In the LETKF, the energy spectrum was evaluated by the mean of the energy spectra in each ensemble member $|E_m(k)|$. Figure 7 shows the energy spectra which are averaged over $t/\Delta t = 76\,000\sim 80\,000$. Since the observation error was given by a Gaussian random noise (i.e. a white noise), the energy spectra of the observation error were almost constant $\Delta E^\mathrm{o}(k) \sim O(10^{-4})$ as seen in the figure. In the low wavenumber regime, the observation error was small enough, and both the nudging and LETKF cases showed good agreement with the ground truth. The difference of the energy spectra can be seen only in the high wavenumber regime $k \gtrapprox 50$. The result of the nudging case was closer to the observation spectrum and deviated from the ground truth; in contrast, the LETKF case showed a good agreement with the ground truth up to the high wavenumber regime. The visible error in the LETKF case was seen only in the dissipation regime $k \gtrapprox 90$.

Zoom In Zoom Out Reset image size

These results might be caused due to the following reasons. Firstly, the nudging potentially induced larger errors because:

• A constant gain parameter γ is assumed, and there is no mechanism to compensate for the influence of the observation error.
• There is no exact transform from the macroscopic variables $(u^\mathrm{o}, v^\mathrm{o}, \rho^\mathrm{o})$ to the velocity distribution function $f_{\alpha}^\mathrm{o}$. In the nudging case, the equilibrium distribution $f_{\alpha}^\mathrm{(eq,o)}$ was employed for reconstructing a reference state vector, which might induce extra errors.

Secondly, in the LETKF, these issues were properly treated:

• The Kalman gain is optimized based on the ensemble statistics by taking account of the observation error, which improves the accuracy much better than the constant gain in the nudging.
• The observation data from neighboring observation points are considered, leading to a non-local DA process. This contains much more information than nudging based on the observation data from a single observation point.

Therefore, in applying the DA to the LBM, the LETKF is more straightforward and accurate than nudging.

4.3.2. Accuracy evaluation at the sparse observation cases.

Next, we investigated the sparse observation cases. The RMSEs observed in the experiments were summarized in figure 8. The numerical experiments were conducted up to p = 64, and points not shown in the figure mean numerically unstable cases. This kind of numerical instability was known as the catastrophic filter divergence problem (Harlim and Majda 2010, Gottwald and Majda 2013, Kelly et al 2015), where small ensemble sizes and sparse observation points produced poor statistics and uncorrelated observation data, leading to inaccurate estimations of the error covariance and thus the Kalman gain. The LETKF cases with $M = 4, p\unicode{x2A7E}32$ and with $M = 16, p = 64$ were unstable and diverged, while M = 64 was stable in all the cases. These results suggest that the shortage of the observation points could be avoided by increasing the ensemble size, and M = 64 was sufficient for the current numerical experiments. Here, it should be noted that the numerical instability of the EnKF in turbulent flows was discussed by Suzuki and Hasegawa (2017, appendix D). They proposed a relaxation factor to mitigate the improper injection energy at the initial few DA cycles, where the randomly distributed initial states of the ensemble give the large amplitude of the Kalman gain. This initial numerical instability appears to be different from the catastrophic filter divergence. In our case, however, the numerical instability occurred in the converged state where the ensemble variance and thus the background error were converged into small enough values. Except for the numerical instability, the DA accuracy in the LETKF was always better than that in the nudging. Here, we scored the DA accuracy against the observation noise $\epsilon^\mathrm{o} = 0.142$, and the result was classified as good when the RMSE was less than $\epsilon^\mathrm{o}$. The nudging showed good results for $p\unicode{x2A7D}8$, while the LETKF showed good results up to p = 32.

Zoom In Zoom Out Reset image size

To clarify the difference between the LETKF and nudging cases, we showed their energy spectra in figure 9 and vorticity maps in figure 10. In the nudging cases, the energy spectra with $p\unicode{x2A7E} 4$ indicated the normal cascade with a steeper power law than that in the ground truth, and this affected also the inverse cascade regime when the observation interval p was larger. This was attributed to the linear spatial interpolation of the observation data. The underestimation of the energy spectra lead to unphysical vorticity maps, in which eddy structures were smeared out. Once such unphysical flows were generated, the nudging could not recover the system in the correct direction. In the nudging cases, the above issues produced a worse result than the above estimate. In contrast, the LETKF cases showed better agreements with the ground truth. The energy spectra showed good agreement except for the energy dissipation regime. The vorticity maps also indicated reasonable agreements, while only the small-scale eddies were smoothed out due to the ensemble average of vorticity fields.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the current numerical experiments, the wavenumber regime of the energy injection due to the forced 2D turbulence was $k_\mathrm{f} = 4$, and larger scale flows were developed via the inverse energy cascade. The ratio of the energy contained below a given wavenumber k can be estimated by

$$\begin{align} R\left(k\right) = \frac{ \sum_{k^{^{\prime}} = 1}^{k} E^\mathrm{t}\left(k^{^{\prime}}\right) }{ \sum_{k^{^{\prime}} = 1}^{k_\mathrm{max}} E^\mathrm{t}\left(k^{^{\prime}}\right) }, \end{align} \tag{ 37 }$$

where $k_\mathrm{max} = 128$ and $E^\mathrm{t}(k)$ denotes the 1D energy spectrum in the ground truth. The value at the injection wavenumber was $R(k_\mathrm{f}) = 89.9\%$, therefore, the DA accuracy was mainly determined by the inverse cascade regime. According to the sampling theorem, the minimum required sampling interval for $k_\mathrm{f} = 4$ is estimated as $p\unicode{x2A7D} k_\mathrm{max}/k_\mathrm{f} = 32$. As seen in the vorticity maps in figure 10, the results in the nudging were smeared out for $p\unicode{x2A7E}32$. In the LETKF, such numerical diffusion of the vorticity map was largely suppressed, and the vortices at the energy containing-scale was well maintained up to p = 64, which exceeded the above criterion. This result suggests that the LETKF could assimilate the flow at finer scale than that of the observation. This might be attributed to the covariance matrices, which takes account of the spatial correlation at finer scale.

Consequently, the DA experiment in this section revealed that the LETKF is highly robust and accurate than the nudging, in particular, at the spatially-sparse observation condition.

4.3.3. Uncertainty and stability evaluation.

The previous section demonstrated the excellent DA capability of the LETKF against the nudging. In this section, we present further details of the ensemble statistics and their relation to the uncertainty and numerical instability in the LETKF.

Figure 11 shows the ensemble-averaged vorticity and that in each ensemble member. The ensemble average smeared out small-scale vortices in each ensemble member, while the individual members kept the small-scale eddies but their positions were slightly different from each other. This indicates that the LETKF did not work as the numerical viscosity; in other words, the amplitudes of the energy spectra were well maintained up to the high wavenumber regime, while only the phase information was ruined due to the shortage of observation points. The evidences are shown in figures 9(left) and 12: The former shows that the amplitudes of the energy spectra agreed well with the ground truth. The latter shows that the phase errors were significantly increased by decreasing the observation points. Here, the phase error was computed as follows:

$$\begin{align} \mathrm{err}\left(k\right) & = \left\langle|\mathrm{err}\left(\boldsymbol{k}^{^{\prime}}\right)|\right\rangle_{k - 1/2 \unicode{x2A7D} |\boldsymbol{k}^{^{\prime}}| \unicode{x2A7D} k + 1/2}, \end{align} \tag{ 38 }$$

$$\begin{align} \mathrm{err}\left(\boldsymbol{k}\right) & = \min |\arg \boldsymbol{E}^\mathrm{t}\left(\boldsymbol{k}\right) - \arg \boldsymbol{E}^\mathrm{DA}\left(\boldsymbol{k}\right) + 2n\pi| \ \forall n\in\left\{-1, 0, 1\right\}, \end{align} \tag{ 39 }$$

where $E^\mathrm{t}(\boldsymbol{k})$ and $E^\mathrm{DA}(\boldsymbol{k})$ are energy spectra of the ground truth and the result of the DA, respectively, and n is introduced for periodicity. The bracket $\langle val \rangle_{cond}$ computes the mean of the value val over the range where the condition cond is satisfied.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 13 is helpful to understand the uncertainty. The figure shows the spaghetti plot, which visually demonstrates the uncertainty of the positions of eddies in the ensemble simulation. The figure indicates the following insights: the results with $p = 1,2,4,8,16, 32$ showed relatively small uncertainty, while those with p = 64 showed the visibly large uncertainty in the small-scale vortices. These findings are also consistent with the quantitative evaluation of the RMSE in figure 8.

Zoom In Zoom Out Reset image size

From these results, it was shown that in the sparse observation case, the disappearance of the small-scale vortices was not due to the numerical viscosity but was due to the ensemble average; this can be seen as a kind of Reynolds averaging over the ensemble simulation. Overall, it was clarified that the LETKF is characterized by a weak numerical dissipation, and the amplitude of the energy spectrum is well maintained even with sparse observation, while the phase error or the time-dependent information is largely affected by the resolution of observation points.

Next, we discuss the cause of the numerical instability at small ensemble sizes, which was shown in figure 8. The origin of the instability might be the analysis increment, namely, the difference between the analysis and the forecast. We evaluated the energy spectrum of the analysis increment,

$$\begin{align} \Delta\left(k\right) = \frac{1}{M} \sum_{m = 0}^{M-1} \sum_{||\boldsymbol{k}^{^{\prime}}|-k|\unicode{x2A7D}1/2} \frac{1}{2}|\left\{\boldsymbol{u}_m^\mathrm{a} - \boldsymbol{u}_m^\mathrm{f}\right\}\left(\boldsymbol{k}^{^{\prime}}\right)|^2. \end{align} \tag{ 40 }$$

Figure 14 shows the results at the moderately sparse observation cases. At p = 32, only the case with M = 4 was numerically unstable, and the corresponding $\Delta(k)$, which was observed just before the divergence, was larger than that with $M = 16, 64$ in the whole wavenumber regime. Especially, it was enhanced in the high wavenumber regime ($k\gtrapprox10$); the amplitude of $\Delta(k)$ in the unstable case (M = 4) was a few orders of magnitudes larger than that in the stable cases ($M = 16, 64$). In addition, the amplitude of $\Delta(k)$ was comparable to the energy spectrum E(k) for $k\gtrapprox 30$. A similar trend was also seen at p = 64; the cases with $M = 4, 16$ were numerically unstable, and they showed larger $\Delta(k)$ than the stable case (M = 64) for $k\gtrapprox10$. In the unstable cases, $\Delta(k)$ was comparable to E(k) for $k\gtrapprox 30$, and was larger than E(k) in the dissipation regime ($k\gtrapprox100$). Hence, such large analysis increment in the high wavenumber regime was the cause of the noisy gain, leading to the spurious energy injection in the high wavenumber regime, and thus the numerical instability.

Zoom In Zoom Out Reset image size

This paper investigated the ensemble DA for the LES by the LBM coupled with the LETKF. We showed that the LETKF enables a straightforward implementation in the LBM, where the velocity distribution function is assimilated based on the macroscopic observable variables such as the density and the velocity. We compared two DA methods of the LETKF and the nudging by the OSSE of the 2D forced isotropic turbulence simulation. Firstly, in the spatially-dense observation case, both the nudging and the LETKF showed high DA accuracy, while the LETKF showed better DA accuracy than the nudging. The duration of the transient period from non-assimilated to fully-assimilated state was also much faster in the LETKF than in the nudging. Although the nudging resulted in a noisy field because of the observation error, the LETKF reproduced a quite accurate energy spectrum up to the high wavenumber regime, where the observation error was relatively large. Secondly, in the spatially-sparse observation case, the nudging gave good results up to the observation points of $32\times32$ against the numerical resolution of $256\times256$. The LETKF, in contrast, showed good results up to the observation points of $8\times8$, which is consistent with the Nyquist wavenumber needed to capture the energy-containing scale. The analysis of the energy spectra also revealed the important features of the LETKF:

• The lack of observation points affected mainly the time-dependent information (i.e. the phase), as shown in figure 12. The mean amplitude was well maintained, as seen in figure 9(left).
• The lack of observation points produced the spurious energy injection at the high wavenumber regime, leading to numerical instability. This phenomenon was known as catastrophic filter divergence (Harlim and Majda 2010, Gottwald and Majda 2013, Kelly et al 2015), which could be suppressed by increasing the ensemble size.

Overall, the LBM-LETKF is highly robust and accurate, provided that the ensemble size is large enough for the given configuration. From these results, we conclude that the LBM-LETKF is an encouraging DA approach for the LES of turbulent flows at least in 2D problems.

In future work, we will study the applicability of the LBM-LETKF in 3D LESs. For this purpose, we are going to address HPC techniques of the LBM-LETKF toward extreme-scale 3D LESs.

Computations were performed on the NVIDIA DGX-2 system at Japan Atomic Energy Agency (JAEA), the HPE SGI8600 supercomputer at JAEA, and the Wisteria-Aquarius supercomputer at The University of Tokyo.

The software that provided results of this study is available on https://github.com/hasegawa-yuta-jaea/LBM2D-LETKF. Output data are available upon reasonable requests.

This work was funded by the JSPS KAKENHI [Grant Numbers JP21K17755, JP22H03599], and the Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures (JHPCN) [Project IDs jh210003, jh210049, jh220030, jh230014].

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Yuta Hasegawa: Methodology, Software, Validation, Formal analysis, Investigation, Data Curation, Writing—Original draft preparation. Naoyuki Onodera: Software, Formal analysis. Yuuichi Asahi: Formal analysis. Takuya Ina: Software. Toshiyuki Imamura: Software. Yasuhiro Idomura: Conceptualization, Writing—Review & Editing, Formal analysis, Supervision.