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

We investigate the applicability of the data assimilation (DA) to large eddy simulations based on the lattice Boltzmann method (LBM). We carry out the observing system simulation experiment of a two-dimensional (2D) forced isotropic turbulence, and examine the DA accuracy of the nudging and the local ensemble transform Kalman filter (LETKF) with spatially sparse and noisy observation data of flow fields. The advantage of the LETKF is that it does not require computing spatial interpolation and/or an inverse problem between the macroscopic variables (the density and the pressure) and the velocity distribution function of the LBM, while the nudging introduces additional models for them. The numerical experiments with 256×256 grids and 10% observation noise in the velocity showed that the root mean square error of the velocity in the LETKF with 8×8 observation points ( ∼0.1% of the total grids) and 64 ensemble members becomes smaller than the observation noise, while the nudging requires an order of magnitude larger number of observation points to achieve the same accuracy. Another advantage of the LETKF is that it well keeps the amplitude of the energy spectrum, while only the phase error becomes larger with more sparse observation. We also see that a lack of observation data in the LETKF produces a spurious energy injection in high wavenumber regimes, leading to numerical instability. Such numerical instability is known as the catastrophic filter divergence problem, which can be suppressed by increasing the number of ensemble members. From these results, it was shown that the LETKF enables robust and accurate DA for the 2D LBM with sparse and noisy observation data.


Introduction
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. Figure 1 (a), (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 & Anthes, 1976;Lakshmivarahan & 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 & Titi, 2003) and in configuration space (Azouani et al., 2014), and were examined in 2D DNSs of isotropic turbulence (Olson & 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, ∼ 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).
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 & Zaki (2022) investigated Figure 1.Conceptual schematics of the continuous DA for a chaotic system with noisy observation.(a) Single-shot simulation without DA, which may predict the wrong state due to the chaotic nature of the system.f true and f calc are the quantities of the true state and the numerical simulation, respectively.(b) Single-shot simulation with the DA (typically, by the nudging method).f obs is the observed data from the true state with some noises.f DA is the predicted state improved by the DA.(c) Ensemble simulation with ensemble Kalman filter.
The filled circles show ensemble variation, where the thin and thick filling colors respectively indicate the range of ensemble data before and after the DA.
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(Evensen, , 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 & Meyers (2021) showed a turbulent flow reconstruction in the atmospheric boundary layer by the 4D-Var DA (Le Dimet & 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 spatiallysparse 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 stateof-the-art supercomputers.In our previous works (Onodera et al., 2018(Onodera et al., , 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 CPUbased supercomputers in the former works (Miyoshi et al., 2014(Miyoshi et al., , 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 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 × 256 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 section 2 and section 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.

Numerical model
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 c α is defined as 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 ∆x and time step width ∆t, namely, c = ∆x/∆t.The governing equation is defined as where ξ = (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 where ρ is the fluid density, and u = (u, v) is the fluid velocity.f eq α is the equilibrium distribution function, defined by where w α is the weighting factor, which is given by τ is the relaxation time.From Chapman-Enskog expansion (Krüger et al., 2017), given the lattice speed c and the time step width ∆t, the relaxation time is chosen to satisfy where ν is the kinetic viscosity and ν 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.ν SGS is then defined as Here, C s = 0.2 is the Smagorinsky constant, and ∆ is the model constant called the filter length, which is set to ∆ = ∆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

Continuous data assimilation models
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 x(ξ) and y(η) 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, Let x f , x a , and y 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, x a (ξ)| 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, x f (ξ)| t=t+s∆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, y o (η)| t=t+s∆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): K is the gain matrix, which varies depending on the DA method.H is the observation operator defined as y = Hx.(v) Update the time step, t ← t + s∆t, and repeat the procedures (i-iv).
In practice, once the vectors x and y are designed, the observation operator H is designed to satisfy y = Hx, and the gain K is finally given.The design of x, y, and K in each DA method is described in the following subsections.

Nudging
Firstly, we introduce a traditional DA approach, which is called nudging (Hoke & Anthes, 1976;Lakshmivarahan & Lewis, 2013), as a reference in this work.The nudging gives the simplest linear interpolation between the state and the observation, i.e., where 0 ≤ γ ≤ 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.H −1 is the inversion of the observation operator, namely, x = H −1 y.To compute H −1 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 where ρo (ξ), ûo (ξ), and vo (ξ) are spatially interpolated variables of the observation ρ o (η), u o (η), and v o (η), respectively.

Local ensemble transform Kalman filter (LETKF)
The local ensemble transform Kalman filter (LETKF) (Hunt et al., 2007) is one of the advanced models of the continuous DA, which is a derived form of the ensemble Kalman filter (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 H † is thus not required.It is noted that the observation operator 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).

Basic model
Let M be the number of ensemble members (i.e. the ensemble size), and x f m be the forecast state vector of the m-th ensemble member, where m = 0, 1, . . ., M − 1.Their ensemble mean and perturbation are then introduced as Similarly, those mapped to the observation space are defined by ȳf where y f m = Hx f m .The analysis state vectors are solved in a latent ensemble space.Let w a m ∈ R M be the analysis state vector in the ensemble space, which is given by the transform x a = X f w a .Its ensemble mean wa and perturbation matrix δW a are then computed by wa = P a δY f⊤ R −1 (y o − ȳf ), ( 21) Here, R is the matrix called the observation error covariance, which is a model constant determined by the observation noise.P a is the analysis covariance in the ensemble space, given by Finally, the analysis state vector in the ensemble space is transformed back to the model space, so that the analysis state vector of the m-th ensemble member x a m is computed as where δW a | m denotes the m-th column of δW a .

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 = O(10) ∼ O(10 3 ) and N = O(10 4 ) ∼ 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 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 • 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 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 where d = |ξ−η|/d cut and d 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 ≥ 2 (see Figure 5).By using R-localization, the observation error covariance R loc and thus the analysis vectors wa , δW 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 & Anderson, 1999), namely, P f inf = βP 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 wa and P a , and thus, Eqs. ( 21) and ( 23) are rewritten as 4. Experiment with 2D forced isotropic turbulence

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π × 2π, double periodic boundaries, and the grid number of 256 × 256.The acceleration is given by the superposition of the energy injection and the friction, F = F i + F µ .Here, the energy injection F i is given by where k is the wavenumber vector, i z is the unit vector in the direction of z-axis, k f and ∆k are respectively the center and width of the energy injection wavenumber regime, r 1 and r 2 are Gaussian random vectors with zero-mean and the standard deviation being ∆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 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 & 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 0 and the reference density ρ 0 are unity.The time-stepping width ∆t is adjusted so that the Courant number to the reference velocity is set to 0.05.In these parameters, the SGS viscosity ν SGS has the value of the same order as the kinetic viscosity ν.
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: where k is the 2D wavenumber, k is the 1D wavenumber for the shell average, û(k) is the Fourier transform of u.There exists the energy spectrum with the dual cascade, composed of the inverse cascade for k < k f , and the normal cascade for k f < k < k d .The dissipation regime is also observed in k > k d , where k d ≈ 90.Although the dual cascade was observed, the normal cascade does not follow the power law of k −3 and is characterized by the steeper slope ≈ 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 −3 and k −5 for the moderate resolution and Reynolds number (Tran & Shepherd, 2002;Tran & Bowman, 2003).The power lows of the normal cascade in former numerical studies were also steeper than k −3 (Herring & Mcwilliams, 1985;Legras et al., 1988;Scott, 2007;Tsang & 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 (cf. the work by Boffetta (2007): the mesh resolution up to 32, 768 2 resulted in the power low of ≈ 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 −3 and adopt this numerical setup as the minimal benchmark problem.

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 pert α at t = 0, where u rand , v rand are given by Gaussian random numbers with zero-mean and the standard deviation being the reference velocity u 0 .f eq α (ρ 0 , u rand , v rand ) is the equilibrium distribution function with the constant density ρ 0 and the random velocity (u rand , v rand ), and δ 0 is the perturbation magnitude.Here, the time t = 0 is not the initial step of the simulation; the original distribution function f α,0 is prepared by the spin-up simulation for 10 6 ∆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 where u pert and u are the velocities in perturbed and unperturbed simulations, respectively.
Following the chaotic nature of turbulence, the error grows exponentially; then, |δ(t)| may be approximated by 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 −5 to 10 −2 .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 δ(t) ∼ 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 λ ∼ 4.64 × 10 −4 ∆t −1 , and then, the predictability time was estimated as 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∆t ⪅ t pred /10.

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 o , v o , ρ o ) were measured with artificial noises on equidistant sampling points in the ground truth.Here, the density ρ o was also measured because the density in the LBM reflects the pressure (P = ρ/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 o /v o and ρ o were 10% of u 0 and 1% of ρ 0 , 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 < ±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×256/p observation points.The observation condition was varied from the dense observation case, p = 1, to the sparse observation cases, p = 2, 4, . . ., 64.
• The temporal interval of the observation was chosen to be the same as the interval of the continuous DA, s = 200.

Forecast ensemble Analysis ensemble
Data assimilation (DA)  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 ≤ γ ≤ 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 cut = p∆x, which guarantees that each computational grid is affected by ∼ O(10) observation points within the distance of 2d cut (see Figure 5).

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, where u t is the velocity in the ground truth state, and u DA is the velocity obtained by each DA approach.In the LETKF, the ensemble mean, u DA = M −1 M −1 m=0 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 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 o , v o , thus, the observation noise of the velocity magnitude was 0.1 × √ 2 = 0.142.The converged value is calculated by the time-averaging over the converged state for 70, 000 ≤ t/∆t ≤ 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 ≥ 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/∆t = 10, 000 ∼ 20, 000, while the nudging case required ∼ 70, 000∆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/∆t ⪆ 70, 000); thus, the gain parameter might be too small for the initial transient period (t/∆t < 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.
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/∆t = 76, 000 ∼ 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 ∆E o (k) ∼ 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 ⪆ 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 ⪆ 90.
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 o , v o , ρ o ) to the velocity distribution function f o α .In the nudging case, the equilibrium distribution f (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 & Majda, 2010;Gottwald & 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 ≥ 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 & 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 ϵ o = 0.142, and the result was classified as good when the RMSE was less than ϵ o .The nudging showed good results for p ≤ 8, while the LETKF showed good results up to p = 32.
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 ≥ 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.In the current numerical experiments, the wavenumber regime of the energy injection due to the forced 2D turbulence was k 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 where k max = 128 and E t (k) denotes the 1D energy spectrum in the ground truth.The value at the injection wavenumber was R(k 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 f = 4 is estimated as p ≤ k max /k f = 32.As seen in the vorticity maps in Figure 10, the results in the nudging were smeared out for p ≥ 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.

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 (k) (k) . Energy spectra of the analysis increment ∆(k), compared with E(k).
Left shows the cases at p = 32, where ∆(k) was observed at t/∆t = 5400, just before the divergence of the case with M = 4. Right shows the cases at p = 64, where ∆(k) was observed at t/∆t = 3200, just before the divergence of the cases with M = 4, 16.
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 Figure 9 (left) and Figure 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: where E t (k) and E DA (k) are energy spectra of the ground truth and the result of the DA, respectively, and n is introduced for periodicity.The bracket ⟨val⟩ cond computes the mean of the value val over the range where the condition cond is satisfied.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.
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, 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 ∆(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 ⪆ 10); the amplitude of ∆(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 ∆(k) was comparable to the energy spectrum E(k) for k ⪆ 30.A similar trend was also seen at p = 64; the cases with M = 4, 16 were numerically unstable, and they showed larger ∆(k) than the stable case (M = 64) for k ⪆ 10.In the unstable cases, ∆(k) was comparable to E(k) for k ⪆ 30, and was larger than E(k) in the dissipation regime (k ⪆ 100).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.

Conclusion
This paper investigated the ensemble data assimilation (DA) for the large eddy simulation (LES) by the lattice Boltzmann method coupled with the local ensemble transform Kalman filter (LBM-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 observing system simulation experiment (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 fullyassimilated 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 × 32 against the numerical resolution of 256 × 256.The LETKF, in contrast, showed good results up to the observation points of 8 × 8, 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 & Majda, 2010;Gottwald & 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 high-performance computing (HPC) techniques of the LBM-LETKF toward extreme-scale 3D LESs.

Figure 2 .
Figure 2. Typical energy spectrum and vorticity map.

Figure 3 .
Figure 3. Error growths with various magnitudes of the initial perturbations.(Left) time history of the errors.(Right) vorticity maps at the initial and final states with/without the perturbation (δ 0 = 10 −5 ).

Figure 4 .Figure 5 .Figure 6 .
Figure 4. Schematic illustration of the observing system simulation experiment (OSSE) for the 2D turbulence simulation.(Top left) Ground truth.(Bottom left) Forecast ensemble, which is the set of trial runs to reproduce the ground truth.Since the ground truth is invisible and the system is chaotic, the forecast ensemble may result in different states.(Top middle) Observation, which is created from the ground truth with artificial observation noises.(Bottom right) the dataassimilated state or the analysis ensemble, that is generated using the observation data to reproduce the ground truth.

Figure 7 .
Figure7.Energy spectra in the OSSEs with the nudging and with the LETKF at the observation interval p = 1.Also shown are the energy spectra of the ground truth and its observation data with artificial noises.
Figure 8. Dependencies of the RMSEs on the interval p of spatially-sparse observation points.|δ(∞)| refers to the time average of |δ(t)| over the converged state in 70, 000 ≤ t/∆t ≤ 80, 000.The black dotted line shows the RMS of the observation noise, ε o = 0.142.Missing points in the LETKF indicate the occurrence of catastrophic filter divergence.

Figure 11 .
Figure 11.Vorticity maps for the ensemble mean and the individual ensemble members at the ensemble size of M = 64 and the observation interval of p = 64.The ensemble mean shows the same result as Figure 10 (bottom, right).

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