Improved State Estimation for turbulent flows
combining Ensemble Data Assimilation and Machine
Learning

Miguel M Valero, Marcello Meldi

► To cite this version:

Miguel M Valero, Marcello Meldi. Improved State Estimation for turbulent flows combining En-semble Data Assimilation and Machine Learning. Computers and Fluids, 2025, 300, pp.106757. 10.1016/j.compfluid.2025.106757. hal-05488619

HAL Id: hal-05488619

https://hal.science/hal-05488619v1

Submitted on 2 Feb 2026

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Improved State Estimation for turbulent flows
combining Ensemble Data Assimilation and Machine
Learning

Miguel M. Valero<sup>a,*</sup>, Marcello Meldi<sup>a</sup>

<sup>a</sup>Univ. Lille CNRS ONERA Arts et Métiers ParisTech Centrale Lille UMR 9014-LMFL- Laboratoire de Mécanique des fluides de Lille - Kampé de Feriet F-59000 Lille France

**Abstract**

A novel strategy is proposed to improve the accuracy of state estimation and the reconstruction from low-fidelity models and sparse data from sensors. This strategy combines ensemble Data Assimilation (DA) and Machine Learning (ML) tools, exploiting their complementary features. ML techniques rely on the data produced by DA methods during analysis phases to train physics-informed corrective algorithms, which are then coupled with the low-fidelity models when data from sensors is unavailable. The methodology is validated via the analysis of the turbulent plane channel flow test case for Reτ ≈ 550. Here, the low-fidelity model consists of coarse-grained simulations coupled with the Immersed Boundary Method (IBM), while observation is sampled by a highly refined body-fitted calculation. The analysis demonstrates the capabilities of the algorithm based on DA and ML to accurately predict the flow features with significantly reduced computational costs. This approach exhibits potential for future synergistic applications of DA and ML, leveraging the robustness and efficiency of ML models alongside the physical interpretability ensured by DA algorithms.

Keywords: EnKF, ML, RFR, IBM, wall turbulence

Email address: miguel.martinez_valero@ensam.eu (Miguel M. Valero)

<footnote>*Corresponding author</footnote>
# 1. Introduction

Advances in the development of numerical techniques in fluid mechanics are rapidly removing several barriers to the analysis of realistic flow configurations of industrial interest. Current applications of data-driven techniques in Computational Fluid Dynamics (CFD), including Data Assimilation (DA) [1, 2] and Machine Learning (ML) [3, 4], have provided accurate results when combined with simulation strategies such as Reynolds-Averaged Navier-Stokes (RANS) approaches [5, 6, 7, 8, 9, 10, 11, 12] and Large Eddy Simulation (LES) [13, 14, 15, 16, 17, 18, 19, 20, 21]. With the progressive increase in computational resources available at computational centres, one can envision the blossoming of futuristic applications such as Digital Twins (DTs) [22, 23, 24, 25, 26] for realistic fluid mechanics cases in the coming decades.

Among the issues that must be challenged to reach such a level of technological maturity, one central aspect deals with the optimisation of the performance of data-driven tools when only sparse data (observation) is available. For ML tools, and in particular, for deep learning, a limited amount of data for training and validation can lead to a lack of convergence of the algorithms, severely impacting the accuracy and robustness of the prediction. The techniques developed to handle this problem are usually referred to as *frugal Artificial Intelligence* [27, 28, 29, 30]. Among these, DA is a very efficient complementary tool. In fact, approaches such as the Ensemble Kalman Filter (EnKF) [1, 31] can be used to produce a complete state estimation on the physical domain investigated, combining the prediction of a model and sparse observation. This state estimation can be very rich in terms of data produced, providing complete information to the ML tools for training. This strategy can, therefore, bypass the potential scarcity of observation. Rigorous mathematical frameworks have been proposed in the literature to combine DA and ML [32, 33, 34], and applications using CFD have been reported [35, 36, 37, 38]. However, the prediction of the DA analysis itself can also be severely degraded because of the scarcity of data. This problem is observed, for example, if observation is available only over a limited time window. During the aforementioned window, DA can succeed in obtaining an accurate prediction, thanks to the combination of the state estimation and the resulting parametric optimisation of the system. However, once the observation is no longer available, the DA state update cannot be performed. This can lead to progressive degradation of the results, particularly if the DA
optimisation of the numerical model is not entirely satisfying. Usually, this
problem happens because of intrinsic limitations in the structural features of
the model itself, such as relying on the eddy viscosity hypothesis for RANS
closures. Therefore, if the DA state augmentation provides poor results in
terms of state estimation, the ML training phase will also be sub-optimal. In
the present work, we explore the usage of ML tools to augment the predictive
capabilities of an EnKF. To this purpose, a Random Forest Regression (RFR)
algorithm [39] is trained to reproduce the model optimisation as well as the
state update obtained in the analysis phase of the EnKF. This algorithm is
then integrated as a black-box solver for the simulation of a permanent flow
in order to assess its accuracy. The validation of the algorithm is performed
via the analysis of the turbulent plane channel flow for $Re_τ \approx 550$. This
test case exhibits fully developed turbulent structures, and it is, therefore,
a complex case to handle for DA and ML algorithms despite the simplicity
of the geometry. Here, the RFR tools are used for two different purposes.
First, they are trained to mimic the effect of an Immersed Boundary Method
(IBM), whose free coefficients are optimised in the DA analysis phase. Sec-
ond, they are trained to learn the nonlinear state update process performed
by the EnKF. It is shown that the global ML tool matches the accuracy of
the complete DA strategy, but it is not affected by the time window of inves-
tigation and by the availability of further data. These findings highlight the
efficiency of the methodology for the study of cases for which data can be
limited in time, but for which analysis over significantly long time windows
may be required. This is, for instance, the case of several internal flows, such
as those in combustion engines or turbomachinery.

The article is structured as follows. In §2, we introduce the numerical
equations and the numerical background, with a particular emphasis on the
data-informed IBM developed using the DA algorithm, along with a presen-
tation of the ensemble-learning ML method employed. In §3, the test case
for analysis is presented, highlighting the limitations of the DA algorithm
and how ML techniques can address these challenges. Later, in §4, we ex-
plore the potential of various ML models as a complete alternative to the DA
tool, and we test them under different conditions from the training dataset
in §5. Finally, in §6, we provide conclusions and suggest directions for future
research.
## 2. Numerical ingredients

The numerical tools used in the present analysis are now briefly introduced. These tools include the strategies for CFD simulations as well as the DA and ML techniques, which are combined.

### 2.1. Numerical solver and Immersed Boundary Method

The numerical simulations are performed using a discretisation of the Navier-Stokes equations for incompressible flows and Newtonian fluids:

$$ \nabla \cdot \mathbf{u} = 0 \qquad (1) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot (\mathbf{u}\mathbf{u}) = -\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}_P \qquad (2) $$

Here, $\mathbf{u}$ is the velocity field, $p$ is the pressure normalised over the density $\rho$, $t$ is the time, and $\nu$ is the kinematic viscosity of the fluid. The coordinates $\mathbf{x} = (x, y, z)$ indicate the physical frame of reference. The term $\mathbf{f}_P$ represents a volume forcing affecting the momentum equation. In this analysis, volume forcing is modelled via the Immersed Boundary Method to consider the presence of the physical wall in the domain. The source term $\mathbf{f}_P$ included in Equation 2 is a classical penalisation scheme [40]:

$$ \mathbf{f}_P = \begin{cases} \mathbf{0} & \text{if } \mathbf{x} \in \Omega_f \\ -\nu \mathbf{D}(\mathbf{u} - \mathbf{u}_{ib}) & \text{if } \mathbf{x} \in (\Sigma_b \cup \Omega_b) \end{cases} \qquad (3) $$

$\Omega_f$, $\Omega_b$, and $\Sigma_b$ correspond to the fluid, solid, and interface regions, respectively. $\mathbf{u}_{ib}$ is here the target velocity, and it is used to obtain a Dirichlet velocity condition at the wall. The coefficients $D_{ij}$ of the tensor $\mathbf{D}$ control the intensity of the volume source term. Guidelines in the literature indicate that their value can be adjusted to optimise performance while ensuring the stability of the numerical algorithm [41, 42]. However, optimised values are often highly sensitive to the local features of the flow (in particular to attached or separated boundary layers) and therefore, global optimisation can provide unsatisfying results.

The numerical discretisation is based on a Finite Volume (FV) approach using a PISO scheme [43, 44, 45, 46]. This algorithm performs a recursive resolution of the discretised set of equations obtained from Equations 1-2 until convergence is obtained. The solver used is the PisoFoam algorithm included in the open-source platform OpenFOAM [47]. Space derivatives are
discretised using second-order centred schemes. A second-order backwards scheme is applied to advance the solution in time. Details about the grid resolution and the time steps are provided in the test case description in §3.1.

## 2.2. Data-informed Immersed Boundary Method using Data Assimilation

Discussion in §2.1 indicated that a global calibration of the IBM coefficients $D_{ij}$ may lead to inaccurate results, as optimised values are related to the local features of the flow at the level of the specific grid element. To this purpose, Valero & Meldi [48] developed a data-informed procedure for the calibration of the IBM coefficients $D_{ij}$. This procedure relies on an ensemble of $N_e$ model realisations and high-fidelity observation, which were combined using a Data Assimilation (DA) technique based on the Ensemble Kalman Filter (EnKF). The two main ingredients used in this data-driven approach are:

* The discrete dynamic model $\mathcal{M}$ corresponds to the unsteady IBM numerical solver presented in §2.1. This model offers a quasi-continuous representation of the flow field $\mathbf{u}$, though its accuracy is affected by the grid refinement selected.

* The observation $\mathbf{Y}$ is obtained by sampling instantaneous data from a high-fidelity, body-fitted Direct Numerical Simulation (DNS). The sensors would provide local values of the physical field (velocity at the wall) as well as global quantities (time evolution of the friction coefficient $C_f$).

Both sources of information are affected by uncertainties, which are here modelled as unbiased Gaussian distributions. The EnKF optimally adjusts the predicted flow field $\mathbf{u}$ to balance the dynamics of both sources of information according to their respective uncertainties. The key elements of the procedure are now discussed, while a complete description of the method is provided in Appendix A. The EnKF scheme is structured in two steps:

1. A forecast ($f$) step where the physical state $\mathbf{u}^f$ is advanced in time using the numerical model obtained by the discretisation of Equations 1–2. Using the DA formalism, this forecast step is expressed as a time advancement of the flow state $\mathbf{u}$ from the time step $k-1$ to $k$ using the model $\mathcal{M}$:

$$ \mathbf{u}_{i,k}^f = \mathcal{M}(\theta_{i,k:k-1}) \mathbf{u}_{i,k-1} \quad (4) $$
Here, the time advancement is considered for each ensemble member i. The array $\theta$ includes free parameters that can be modified to control the model's prediction $M$. An example can be numerical values in boundary conditions or the local values of the coefficients $D_{ij}$ governing the IBM.

2. If observation is sampled at the time step $k$, an analysis (a) phase is performed where the DA procedure uses the available information to improve the flow prediction and the model's accuracy. To do so, the model prediction at the forecast is combined with observation via the so-called Kalman gain $K$:

$$ \mathbf{u}_{k,i}^{a} = \mathbf{u}_{k,i}^{f} + \mathbf{K}_{k} \left( \mathbf{Y}_{k,i} - \mathcal{H}(\mathbf{u}_{k,i}^{f}) \right) \quad (5) $$

The projection operator $\mathcal{H}$ usually consists of manipulation (interpolation and/or integration) of the physical field predicted by the model to obtain a consistent comparison with the observation sampled by the sensors. Equation 5 is usually used for two distinct operations. The first one, usually referred to as state augmentation or state estimation, updates the physical state $\mathbf{u}^{f \to a}$. The second one, referred to as parametric optimisation, uses results from the previous operation to optimise some model's free parameters $\theta^{f \to a}$ to reduce the discrepancy between model and observation. In this work, the two operations are performed simultaneously thanks to the extended state approach [1].

The analysis phase is performed only when observation $\mathbf{Y}$ is available at time step $k$ during unsteady simulations. The accuracy of the Kalman filter depends on the time window between successive analyses as well as on the size of the time step of the model. For the former, the DA procedure may diverge if the sample rate does not respect the Nyquist theorem with respect to the flow turnover time [15]. On the other hand, the size of the time step of the model must be small enough to ensure efficient physical propagation between analysis phases. In this analysis, the EnKF is performed online, which means that ensemble simulations are not stopped and restarted between analysis phases, but they are put on standby while the DA algorithm performs the state estimation and the parametric optimisation. This efficient algorithmic structure is obtained via a team-developed application, CONES (Coupling OpenFOAM with Numerical EnvironmentS) [12], available on https://github.com/MiguelMValero/CONES. This tool facilitates efficient coupling and field exchange between OpenFOAM modules and
the DA algorithm code, eliminating the need for simulation interruptions or
restarts, which could otherwise incur higher computational costs than the
DA procedure itself.

If the underlying model is unbiased and complete (i.e., its prediction
is exact once its free coefficients are optimised), state augmentation is ex-
pected to become negligible after the EnKF converges, as the physical state
predicted by the model aligns with the high-fidelity space. However, if the
model is biased and cannot exactly reproduce the true solution even with an
optimised parametric description (e.g., turbulence modelling [49, 50]), state
augmentation plays a crucial role in reducing the discrepancy between predic-
tions and observations, thereby enhancing accuracy. Hence, DA prediction
using state augmentation and parametric optimisation is, in principle, more
precise or equivalent to a flow solution obtained solely with a DA-optimised
model. However, several factors affect the choices made for the DA strategy.
One key criterion deals with the computational resources required. Simul-
taneously updating both the system's state and the model's coefficients is
significantly more computationally expensive than optimising the model's
free parameters alone. These differences can span several orders of magni-
tude. The primary reason for this is that the number of degrees of freedom
resolved by the model $n_{DF}$ is usually significantly larger than the number of
free parameters to be optimised $n_\theta$ ($n_{DF} \gg n_\theta$). In 3D simulations for incom-
pressible flows without dimensionality reduction techniques like localisation
[1], $n_{DF}$ is three times the number of cells in the domain, which corresponds
to the velocity vector in each grid element. On the other hand, $n_\theta$ typically
consists of only a few values if a global optimisation of the models is per-
formed. Additionally, the algebraic operations required by the DA algorithm
are computationally intensive, compounded by the need to send and receive
arrays containing a high number of data points and the intrusive modifi-
cations to the solver, both of which further increase the computation time.
Besides, here, the elapsed simulation time $t_{sim}$ and the computation cost C.C
can differ significantly since the former only accounts for a single realisation
(plus the DA algorithm occasionally), whereas the latter must include the
calculation of the $N_e$ realisations, so $C.C. \approx N_e t_{sim}$.

A second criterion affecting the accuracy of the DA algorithm is that the
synergistic combination of system augmentation and parametric optimisation
requires the presence of high-fidelity data within the time interval considered
during the simulation. If, for the sake of simplicity, one considers that the
first analysis phase is performed for $t = 0$, the complete DA scheme with state
<figure><img src="image_2.png" /><figcaption>Figure 1: Difference between the state augmentation and parametric optimisation via the EnKF and the pure parametric optimisation once the model’s free parameters $\theta_k^a$ are converged.</figcaption></figure>

estimation and parametric optimisation can be performed for the observation window $t_y \in [0, t_{max}]$, where $t_{max}$ corresponds to the last set of samples of the time series. For predictions beyond the observation window ($t > t_{max}$), the only way to enhance the flow is by replacing the original model’s parameters with the updated free coefficients from the EnKF ($\theta_k^a$) and running a single realisation, as shown in Figure 1. In the case of the present analysis, a CFD run with the IBM model is first augmented using state estimation and parametric optimisation for $t \le t_{max}$, and a CFD with optimised parameters for the local coefficients $D_{ij}$ is run for $t > t_{max}$. As the IBM model is not complete, one can expect a progressive degradation of the results for $t > t_{max}$, which is investigated in §3.2.

## 2.3. Machine Learning tools: Random Forest Regression

The ML tools used in the present work are now described. The Random Forest Regression (RFR) [39] is an ensemble learning method that builds multiple decision trees during training and outputs the average of their pre-
dictions to enhance accuracy and model robustness. It is widely applied in regression tasks, and it is known for its flexibility and ability to handle both numerical and categorical features. Combining bootstrapping and random feature selection, it aggregates the predictions of multiple decision trees, thus mitigating overfitting and enhancing model stability. Several critical steps are involved in generating decision trees within an RFR model. Initially, one must determine the number *R* of decision trees to be included in the forest. Increasing the number of trees can improve the accuracy and convergence of the algorithm, but also demands more considerable computational resources.

Figure 2 provides a schematic overview of the method. Each decision tree $r \in [1, ..., R]$ is constructed using a subset of the entire training dataset and selects splitting criteria or nodes to progress in tree depth. Commonly, the splitting criterion aims to minimise the variance of the target variable within resulting subsets, i.e., the split point is chosen to minimise the Mean Square Error (MSE) of the child nodes ($MSE_{child}$) and to maximise the reduction in MSE ($\Delta$MSE) from the parent to the child nodes, which is computed as the difference between the MSE in the parent node and the weighted sum of the MSEs in the resulting child nodes. These quantities are calculated using the following equations:

$$ \text{MSE} = \frac{1}{m} \sum_{i=1}^{m} (\mathcal{Y}_i - \hat{\mathcal{Y}}_i)^2 \quad (6) $$

$$ \Delta \text{MSE} = \text{MSE}_{\text{parent}} - \left( \frac{m_{\text{left}}}{m} \text{MSE}_{\text{left}} + \frac{m_{\text{right}}}{m} \text{MSE}_{\text{right}} \right) \quad (7) $$

Here, *m* represents the number of samples of the parent node, $\mathcal{Y}_i$ and $\hat{\mathcal{Y}}_i$ indicate the actual and predicted target value for the sample *i*, and $m_{left}$ and $m_{right}$ refer to the number of samples in the left and right child nodes. The MSEs for the parent, left child, and right child nodes are $MSE_{parent}$, $MSE_{left}$ and $MSE_{right}$, respectively.

At each parent node, the dataset is partitioned into two subsets based on a specific feature $\mathcal{X}_j$ and split point, yielding two child nodes. The left child node typically represents data meeting the split criteria (e.g., feature value $\mathcal{X}_{i,j} \le$ split point), while the right child node does not (e.g., feature value $\mathcal{X}_{i,j} >$ split point). Splitting nodes continues recursively for each child node until predefined stopping criteria are met, such as the maximum depth of each decision tree or the minimum number of samples required to split an
<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic representation of a parent node and the criteria to split it into two child nodes for a random decision tree r.</figcaption></figure>

internal node or form a leaf node. These hyperparameters play a crucial role
in preventing overfitting. By navigating these steps, RFR constructs a robust
ensemble model capable of capturing complex relationships within data while
maintaining generalisability. Additionally, it provides several advantages over
deep learning architectures, such as neural networks (NNs), particularly in
terms of interpretability and input preparation. Unlike NNs, it does not
require feature normalisation and still delivers strong reliability. Moreover,
NNs usually require larger datasets to attain comparable accuracy and are
more susceptible to overfitting issues [51, 52].

In the context of CFD, RFR has recently been employed in the prediction of 2D scenarios of turbulent flows via the blending of RANS models for a NACA profile [53] and other test cases like a turbulent flat plate, an axisymmetric subsonic jet, and a wall-mounted hump [54]. Also, enhanced physics-informed RFR have been applied in the prediction of the Reynolds-stress anisotropy tensor for improved turbulence models in RANS simulations [55, 56, 57].
### 3. Performance of the EnKF when using finite observation time windows

This section discusses the limitations of the EnKF when observation is available only over a limited time window. As previously mentioned, the availability of limited datasets over time can significantly impact the accuracy of flow predictions, particularly in scenarios involving complex modelling challenges, such as turbulence modelling and IBM. First, the test case of investigation is presented. The DA experiment is then detailed, and the sensitivity of the prediction of the EnKF to the availability of observations is assessed. In particular, the performance of the EnKF within the observation window $t \le t_{max}$ (where state augmentation/estimation and parametric optimisation of the model can be performed) is compared with the accuracy obtained after the end of such time window $t > t_{max}$, where state estimation is not available anymore and the model advances in time using exclusively the DA-optimised free parameters.

#### 3.1. Test case: turbulent plane channel flow, $Re_{τ} \approx 550$

The turbulent plane channel flow is a cornerstone benchmark test case in Computational Fluid Dynamics (CFD), offering readily accessible high-fidelity Direct Numerical Simulations (DNS) open databases online [58, 59, 60, 61, 62]. A fully developed turbulent flow oriented in the streamwise direction $x$ is investigated between two parallel static walls positioned at $y = 0$ and $y = 2h$. Here, $y$ denotes the wall-normal direction, $z$ is the spanwise direction, and $h$ is the half-height of the channel. A scheme of the flow is found in Figure 3(a).

The friction Reynolds number investigated here is $Re_{τ} = u_{τ}h/ν \approx 550$, where $u_{τ} = \sqrt{τ_w/ρ}$ represents the friction velocity derived from the wall shear stress $τ_w$. A uniform mass flow rate over time is obtained by integrating a source term into the dynamic equations already available in the libraries of OpenFOAM. Periodic boundary conditions are used for the lateral boundaries. A coarse-grained body-fitted DNS, named hereafter DNS-BF, is conducted with high wall resolution, achieving a wall unit value $y^+ \approx 1$. Additionally, a coarse-grained simulation with the penalty IBM, called DNS-IBM, is simulated with a grid very similar to the one employed for DNS-BF. Notable differences lie in the treatment of the near-wall region: both simulations incorporate the immersed/physical walls at $y = 0, 2h$, yet for the
<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Scheme of the turbulent plane channel flow, and (b) distribution of the cell layers in the wall-normal direction for the solid region $\Omega_b$, the fluid region $\Omega_f$, and the interface regions $\Sigma_b$ for the IBM simulations.</figcaption></figure>

DNS-IBM run, these positions correspond to the centre of the fourth Eulérian mesh element starting from the outer limit of the computational domain, where a no-slip boundary condition is imposed. As shown in Figure 3(b), this configuration results in two cell layers in the wall-normal direction representing the solid $\Omega_b$ region of the domain (grey colour), while three cell layers depict the interface $\Sigma_b$ between the fluid and solid regions (blue colour). These additional layers ensure numerical stability. As expressed in Equation 3, a penalisation term $f_P$ is included in the dynamic equations if the position of the considered grid cell is $x \in \Sigma_b \cup \Omega_b$. The penalisation coefficients $D_{ij}$ included in the tensor $D = DI$ are set with a uniform value of $D = 10^5$ in the mesh elements where $x \in \Sigma_b \cup \Omega_b$, and are set to zero elsewhere, following the work by Valero & Meldi [48]. To enforce the zero-velocity condition at the immersed walls, $u_{ib} = 0$.

In the aforementioned work, a strategy based on the EnKF was proposed to combine high-fidelity data sampled from a well-refined body-fitted DNS (referred to as R-DNS-BF) with the lower-fidelity DNS-IBM model previously presented. It was demonstrated that considering an observation window of 300 $t_A$, where $t_A$ is the advective time defined for this test case as $t_A = h/U_c$, with $U_c$ the time-averaged velocity at the centreline, the DA algorithm successfully achieves accurate flow prediction. This result is accomplished by combining the EnKF state estimation with the local optimisation of the IBM model coefficients $D_{ij}$. In this case, variations of the coefficients $D_{ij}$ in the wall-normal direction $y$ were investigated, which required the optimisation of five coefficients for each component of the volume forcing (see Figure 3(b)).
This last optimisation task reduces the bias of the model, which, therefore,
can follow the time evolution of the observed quantities with increased ac-
curacy. However, in cases where the structural limitations of the model can
preclude an exact representation of the test case investigated, even with an
optimised parametric description, the update obtained by state estimation
can also be beneficial in obtaining increased accuracy. The analysis in §3.2
aims to identify the contribution of these two DA concurring aspects.

3.2. Accuracy of the parametric optimisation using the EnKF

The differences in accuracy between the full DA algorithm (which com-
bines state estimation and parametric optimisation) and numerical simu-
lations with optimised DA models are now discussed for the test case of
interest. These two DA-augmented IBM runs are referred to as DNS-IBM-
DA<sub>s.e.</sub> and DNS-IBM-DA<sub>p.o.</sub>, respectively. In the first case, data from the
high fidelity simulation R-DNS-BF (in the form of instantaneous values of
the friction coefficient C<sub>f</sub> and the wall velocity **u** = **0** on several wall sensors)
are combined with IBM realisations during 2 500 analysis phases, which are
regularly spaced over an observation window of 300 t<sub>A</sub>. During these analy-
sis phases, the flow field is updated via state estimation, and the parameters
D<sub>ij</sub> of the IBM are optimised. The complete results of the investigation are
provided and discussed in the work by Valero & Meldi [48]. For the case
DNS-IBM-DA<sub>p.o.</sub>, a single IBM run is performed, using the optimised coeffi-
cients D<sub>ij</sub> determined by the case DNS-IBM-DA<sub>s.e.</sub>. One could see that the
case DNS-IBM-DA<sub>p.o.</sub> corresponds to a numerical simulation performed after
a DA investigation if, for any reason, the streaming data from the targeted
application is interrupted.

The results of this study are illustrated in Figure 4, where the data is
presented in non-dimensional form ⟨·⟩<sup>+</sup> by employing u<sub>τ</sub> and the viscous
wall unit δ<sub>ν</sub> = ν/u<sub>τ</sub> for normalisation. The + suffix indicates normalisations
using the u<sub>τ</sub> calculated by each simulation, while the ★ suffix indicates nor-
malisation using the u<sub>τ</sub> obtained by the simulation R-DNS-IBM. Statistical
homogeneity is observed across the x and z directions, so the statistical mo-
ments are averaged over these two directions. Further details on the test
cases can be found in Table 1. The mean velocity profile ⟨U⟩<sup>+</sup> = ⟨u<sub>x</sub>⟩/u<sub>τ</sub>
and the components of the Reynolds stress tensor ⟨u'<sub>i</sub> u'<sub>j</sub>⟩<sup>+</sup> = ⟨u'<sub>i</sub> u'<sub>j</sub>⟩/u<sub>τ</sub><sup>2</sup> of
the simulation R-DNS-BF, which is here used as a reference, were already
validated with datasets in the literature [59, 60] in our previous work [48].
One can see in Figure 4(a)–(b) that the mean flow predicted by the data-
informed (DNS-IBM-DA$_{s.e.}$) run is extremely close to the reference simula-
tion, unlike the DNS-BF and DNS-IBM simulations. The combined effect
of state estimation and parametric optimisation is, therefore, successful in
obtaining an accurate prediction. More precisely, DA not only ensures more
precise enforcement of the no-slip condition at the immersed wall compared
to the pure IBM simulation but also enhances the friction velocity $u_τ$, with
a discrepancy of only 1.46% when compared with the value obtained by the
R-DNS-BF run. These enhancements are achieved without incurring the
high computational costs C.C that were required for the R-DNS-BF run.
In Table 1, C.C is normalised over the computational costs of the coarse-
grained body-fitted simulation, indicated as $C.C^* = C.C/C.C_{DNS-BF}$. As
highlighted in §2.2, it is important to observe that state estimation demands
significantly higher computational costs compared to parametric optimisa-
tion alone—over an order of magnitude difference—due to the additional
overhead of the EnKF algorithm and the need to run an ensemble of simu-
lations. Furthermore, when comparing the curves between DNS-IBM-DA$_{s.e.}$
and DNS-IBM-DA$_{p.o.}$, a clear degradation is observed in the latter. While
its accuracy is still significantly better than the one achieved by the simula-
tions DNS-BF and DNS-IBM, the results show the effect of the cumulative
errors produced by the model bias discussed in previous sections, preventing
the achievement of the same level of accuracy that was observed for the DA
solution with state estimation.

Evaluating the second-order statistics, as illustrated in Figure 4(c)–(f),
presents a nuanced challenge. Overall, subtle improvements are observed in
the Reynolds stress tensor components with the DA procedure, with minimal
distinction between DNS-IBM-DA<sub>s.e.</sub> and DNS-IBM-DA<sub>p.o.</sub>. The limited dif-
ferences between the two DA strategies arise because the second-order fields
are passively updated when local velocity values serve as observations. In
both approaches, enhanced accuracy is noted in ⟨u'<sub>y</sub>u'<sub>y</sub>⟩<sup>+</sup>, while slight im-
provements are evident in ⟨u'<sub>x</sub>u'<sub>x</sub>⟩<sup>+</sup> and ⟨u'<sub>x</sub>u'<sub>y</sub>⟩<sup>+</sup>. Conversely, minor degrada-
tions are observed in ⟨u'<sub>z</sub>u'<sub>z</sub>⟩<sup>+</sup>. Upon comparison, the amalgamation of state
estimation with parametric optimisation generally yields more precise curves
than parametric optimisation alone, except for ⟨u'<sub>y</sub>u'<sub>y</sub>⟩<sup>+</sup>.

One can see that the effect of state estimation is, therefore, beneficial for
the global accuracy of the solution, in particular when structural deficiencies
of the model may preclude a completely satisfying optimisation. This is
actually the case for most applications in turbulent flows. In the following,
<figure><img src="image_5.png" /><figcaption>Figure 4: Comparison of the main statistical moments of the velocity field. Results are shown for simulations (···) DNS-IBM-DA<sub>s.e.</sub>, (—) DNS-IBM-DA<sub>p.o.</sub>, (—) DNS-BF, (···) DNS-IBM, and (—) R-DNS-BF.</figcaption></figure>

we develop ML techniques that target the emulation of the DA procedures.

4. Machine Learning for state estimation and parametric optimi-
sation

The usage of DA and ML tools in the framework of scale-resolving stream-
ing data is now presented. The resulting models can be trained and poten-
tially updated online when observation for DA is available, and they can be
used as a DA surrogate when observation is not available. Two strategies
are proposed here. The first strategy, which is conceptually the simplest,
focuses on the DA parametric optimisation phase. Specifically, a black-box
IBM model, hereafter referred to as IBM<sub>RFR</sub>, is designed to emulate the DA-
augmented IBM tool. The second ML strategy, termed SE<sub>RFR</sub>, captures the
nonlinear effects of state estimation, offering the capability to update the
flow state during time windows when observations are unavailable.
<table><caption>Table 1: Summary of the numerical runs considered for the study. †For DNS-IBM-DA<sub>s.e.</sub>, the computational costs are also normalised with respect to the ensemble members N<sub>e</sub> to represent the elapsed simulation time. This normalisation is expressed as C.C<sup>*</sup>/(N<sub>e</sub> + 1), where C.C<sup>*</sup> represents the computational costs, and an additional CPU core is assumed for performing the analysis phases with the EnKF.</caption><thead><tr><th></th><th>N<sub>x</sub> &times; N<sub>y</sub> &times; N<sub>z</sub></th><th>&Delta;x<sup>*</sup></th><th>&Delta;z<sup>*</sup></th><th>&Delta;y<sup>*</sup><sub>min</sub></th><th>&Delta;y<sup>*</sup><sub>max</sub></th><th>L<sub>x</sub>/h</th><th>L<sub>z</sub>/h</th><th>u<sub>&tau;</sub></th><th>C.C<sup>*</sup></th></tr></thead><tbody><tr><td>R-DNS-BF</td><td>1 024 &times; 256 &times; 512</td><td>9.9</td><td>6.6</td><td>1</td><td>11.2</td><td>6&pi;</td><td>2&pi;</td><td>0.0480</td><td>&asymp; 12k</td></tr><tr><td>DNS-BF</td><td>128 &times; 61 &times; 64</td><td>39.5</td><td>26.3</td><td>1.08</td><td>52</td><td>3&pi;</td><td>&pi;</td><td>0.0430</td><td>1</td></tr><tr><td>DNS-IBM</td><td>128 &times; 64 &times; 64</td><td>39.5</td><td>26.3</td><td>1.3</td><td>52</td><td>3&pi;</td><td>&pi;</td><td>0.0375</td><td>1.03</td></tr><tr><td>DNS-IBM-DA<sub>s.e.</sub></td><td>128 &times; 64 &times; 64</td><td>39.5</td><td>26.3</td><td>1.3</td><td>52</td><td>3&pi;</td><td>&pi;</td><td>0.0487</td><td>37.60<sup>&dagger;</sup></td></tr><tr><td>DNS-IBM-DA<sub>p.o.</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.0509</td><td>1.03</td></tr><tr><td>DNS-IBM-ML<sub>p.o.</sub></td><td>128 &times; 64 &times; 64</td><td>39.5</td><td>26.3</td><td>1.3</td><td>52</td><td>3&pi;</td><td>&pi;</td><td>0.0516</td><td>7.03</td></tr><tr><td>DNS-IBM-ML<sub>s.e.</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.0493</td><td>8.78</td></tr></tbody></table>

## 4.1. Strategy one: Machine Learning experiment to obtain a black-box IBM tool

The ML strategy proposed here mimics the parametric optimisation procedure performed by the EnKF in order to obtain a black-box tool providing the same performance. DA is an efficient complementary tool for ML, as the instantaneous flow data can be used to train the algorithm. Figure 5 illustrates the values of the coefficients D<sub>ij</sub> employed in the IBM simulations. One can see the differences produced by the DA optimisation process when comparing the prior state of simulation DNS-IBM and the final state used for the run DNS-IBM-DA<sub>p.o.</sub>. By exploiting the rich databases obtained by the DA algorithm, a comprehensive database of local penalty forcing values f<sub>P</sub> can be obtained over several time steps. This quantity serves as an output for the black-box ML tool, and therefore it is selected as labelled data or dependent variable Y. The features or independent variables X selected for this database include the local velocity field values u = (u<sub>x</sub>, u<sub>y</sub>, u<sub>z</sub>) and the position of the considered cell layer in the wall-normal direction y/h. An axial symmetry condition is enforced at half-height of the channel, where mesh elements near the top wall y<sub>top</sub>/h contribute to the database via the expression y/h = 2 - y<sub>top</sub>/h. Given the abundance of continuous data, leveraging the Random Forest Regression (RFR) ensemble supervised learning method appears to be a promising starting point. This method of-
<figure><img src="image_6.png" /><figcaption>Figure 5: Penalty coefficients $D_{ij}$ improved by means of (—) DNS-IBM-DA$_{s.e.}$, concerning the ones employed in the classical penalisation (····) DNS-IBM. The standard deviation $\sigma_{ij}$ of the prior for the initial Gaussian perturbation of the ensemble is represented in light grey ($\sigma_{ij} = 0.05 D_{ij}$).</figcaption></figure>

fers resilience against overfitting and outliers compared to other regression algorithms like linear regression, and can effectively handle nonlinear relationships in the data. Implementation is simplified through the C++ library `dlib` [63], making integration into the CFD solvers of OpenFOAM straightforward. Also, it employs all available CPU cores of the system to serialise the ML models, optimising costs during the learning procedure.

Preliminary tests have shown that, from a few million sampling data, the influence of hyperparameters on training quality is not particularly strong. Consequently, hyperparameter values are selected with a focus on preserving computational efficiency. For the test cases investigated here, the number of decision trees is set to $R = 100$, where each individual decision tree considers a third of the whole dataset for split criteria in the parent nodes. The minimum sampling size to form a leaf node is set to 5, mitigating the risk of overfitting.

The training data consists of 76 consecutive instantaneous snapshots of one of the ensemble members of the run DNS-IBM-DA$_{s.e.}$, which are obtained over a time span of $1.5 t_A$, with a time step of $\Delta t = 0.02 t_A$. The time window for sampling here is set in the range [$298.5 t_A$, $300 t_A$] in order to grant convergence of the optimisation of the coefficients $D_{ij}$. The choice of the specific ensemble member is not relevant, as each of them is a different Monte-Carlo realisation of the system uncertainty, and they all provide the same statistical information. The dataset used for training consists of $m_T = 6225920$ sampling data points. The training process is conducted online, generating three distinct models to predict each of the components of the forcing term $\hat{f}_P = (\hat{f}_{P_x}, \hat{f}_{P_y}, \hat{f}_{P_z})$.
<table><caption>Table 2: Metrics computed for each component of the forcing term $\hat{f}_P$.</caption><thead><tr><th></th><th colspan="3">IBM<sub>RFR</sub></th></tr><tr><th></th><th>f&#x0302;<sub>P<sub>x</sub></sub></th><th>f&#x0302;<sub>P<sub>y</sub></sub></th><th>f&#x0302;<sub>P<sub>z</sub></sub></th></tr></thead><tbody><tr><th>R<sup>2</sup></th><td>99.72%</td><td>67.02%</td><td>88.28%</td></tr><tr><th>NRMSE</th><td>25.38%</td><td>324.40%</td><td>55.71%</td></tr></tbody></table>

The model performance is assessed via two metrics utilising $p = 20\%$ of the entire training dataset $m_T$: the coefficient of determination $R^2$ to evaluate how well the model captures the variance in data and the Normalised Root Mean Square Error (NRMSE) to provide insight into the average of the prediction errors. For normalisation, selecting the interquartile range—$Q_1$ and $Q_3$ representing the 25th and 75th percentiles of the actual penalty force $f_P$—is advantageous due to its robustness against outliers. Their expressions are respectively:

$$R^2 = 1 - \frac{\sum_{i=1}^{p} (\mathcal{Y}_i - \hat{\mathcal{Y}}_i)^2}{\sum_{i=1}^{p} (\mathcal{Y}_i - \bar{\mathcal{Y}})^2} \quad (8)$$

$$NRMSE = \frac{1}{Q_3 - Q_1} \sqrt{\frac{\sum_{i=1}^{p} (\mathcal{Y}_i - \hat{\mathcal{Y}}_i)^2}{p}} \quad (9)$$

Preliminary investigations and tests dealing with the ML training show that the convergence of the obtained black-tool IBM depends on the volume forcing component. More precisely, the forcing component working in the streamwise $x$ direction, which is the largest in magnitude, exhibits significantly better diagnostics than the other two components. This suggests a significant presence of noise in the samples in the wall-normal $y$ and spanwise $z$ directions. This observation is corroborated by Figure 6, where several lower orders of magnitude are evident for the source term in these directions compared to the streamwise value. Therefore, the decision is made to solely predict the streamwise component of the forcing, yielding $\hat{f}_P = (\hat{f}_{P_x}, 0, 0)$, defining this particular simulation as DNS-IBM-ML<sub>p.o.</sub>.

This model incorporates the normalisation of the features and the labelled data, ensuring its applicability across diverse conditions that are not limited only to the training dataset. This is performed, first, by scaling the
<figure><img src="image_7.png" /><figcaption>Figure 6: Averaged penalty source term $\langle \boldsymbol{f}_P \rangle$ of the database employed for training. In light grey, we represent the standard deviation.</figcaption></figure>

features $(\cdot)^*$ relative to the global friction velocity determined through DA with parametric optimisation, i.e., $u_{\tau_{\text{DA/p.o.}}} = 0.0509$ (refer to Table 1). Consequently, for each local sample indexed by $i$, the feature vector is $\mathcal{X}_i = (u_{x_i}^*, u_{y_i}^*, u_{z_i}^*, y_i^*) = (u_{x_i}/u_{\tau_{\text{DA/p.o.}}}, u_{y_i}/u_{\tau_{\text{DA/p.o.}}}, u_{z_i}/u_{\tau_{\text{DA/p.o.}}}, (y_i/h) u_{\tau_{\text{DA/p.o.}}}/\nu)$. Second, the output of the model consists of the forcing term normalised by considering $\boldsymbol{f}_{P_i}^* = \boldsymbol{f}_{P_i} \delta_{\nu_{\text{DA/p.o.}}}/u_{\tau_{\text{DA/p.o.}}}^2$.

The generation of the black-box IBM via the RFR model (IBM<sub>RFR</sub>) is shown in Figure 7 for the simulation DNS-IBM-ML<sub>p.o.</sub>. In this representation, each local sample $i$ of the training dataset includes the instantaneous velocity values and the wall-normal distance at the cell centre of the respective mesh element under consideration. The predicted force $\hat{\boldsymbol{f}}_P^*$ is estimated and converted to $\hat{\boldsymbol{f}}_P$ on the fly within the CFD formalism (see Appendix B for more details about the implementation of RFR inside the PISO algorithm).

## 4.2. Validation of Random Forest Regression in the training dataset

The black-box IBM<sub>RFR</sub> is now assessed using conditions consistent with the training data. The analysis of first-order and second-order statistics, which are shown in Figure 8, indicates that the performance achieved by the ML-generated model closely matches that of DA with parametric optimisation, albeit with a slight degradation observed for the component $\langle u'_y u'_y \rangle^+$. In any case, the black-box IBM is systematically more accurate than the coarse-grained body-fitted (DNS-BF), and IBM runs for most of the statistical moments. This highlights the excellent convergence of ML training thanks to the richness of data provided by the DA algorithm. In addition, as shown in Table 1, the simulation elapsed time is approximately five times
<figure><img src="image_8.png" /><figcaption>Figure 7: Overview of RFR algorithm for the prediction of the streamwise forcing term $\hat{f}_{P_x}$.</figcaption></figure>

smaller in DNS-IBM-ML<sub>p.o.</sub> compared to DNS-IBM-DA<sub>s.e.</sub>. This discrepancy increases to roughly 5 ($N_e + 1$) when factoring in total CPU time/cost $C.C^*$. Also, as shown in Figure 8(b), the zero-velocity condition at the wall is properly constrained and minimal discrepancies are observed for the prediction of $u_\tau$ (around 1.38%) when compared with the run DNS-IBM-DA<sub>s.e.</sub>.

The investigation of the instantaneous flow patterns shown in Figure 9, highlights a close similarity between the Q-criterion near-wall isocontours of DNS-IBM-ML<sub>p.o.</sub> and those of the DA strategy. The training performed using instantaneous flow data allows the IBM<sub>RFR</sub> model to obtain a qualitatively correct representation of the main flow structures.

However, the analysis of the time spectra of the fluctuating velocity $u'$ at different localised sensors, which are obtained via Fast Fourier Transform [FFT, 64], shows some limitations of IBM<sub>RFR</sub>. Here, the dimensionless wavenumber $\kappa^+ = 2\pi f\nu/(U_c u_\tau)$ is estimated from the FFT frequency $f$ by assuming Taylor's frozen turbulence hypothesis [65], consistent with the approach in [48]. The results, which are presented in Figure 10, show that the predicted spectrum is very similar to the one obtained by the runs DNS-IBM and DNS-BF. The prediction for this quantity is, therefore, significantly less accurate than the run IBM-DNS-DA<sub>p.o.</sub> run, which systematically gets closer to the reference results of simulation R-DNS-BF. Several factors may contribute to the ML model's degradation relative to the simulation used for training. First, the $u_\tau$ value used for normalisation during the training task is the one obtained in the DA run, but, as shown in Table 1, a small yet non-negligible discrepancy is observed for the prediction obtained with the
<figure><img src="image_9.png" /><figcaption>Figure 8: Comparison of the main statistical moments of the velocity field. Results are shown for simulations (---) DNS-IBM-ML<sub>p.o.</sub>, (—) DNS-IBM-DA<sub>p.o.</sub>, and (—) R-DNS-BF.</figcaption></figure>

IBM<sub>RFR</sub>. In addition, the choice to simplify the black-box volume forcing term so that $\dot{f}_P = (\hat{f}_{P_x}, 0, 0)$ may provide some sources of inaccuracy in developing instantaneous physical strategies. However, one can see from the comprehensive analysis of the results that the IBM<sub>RFR</sub> globally behaves like the DA prediction used for training, and its behaviour is never worse than the simulations DNS-BF and DNS-IBM performed on the same grid.

This investigation validates the possibility of using instantaneous data to train a black-box IBM model. This task is of moderate difficulty because, even if the coefficients $D_{ij}$ are varying in the $y$ direction for the DA algorithm, the volume forcing is still a linear function of the velocity field. In the next section, a more complex ML tool is developed in order to capture the nonlinear dynamics associated with the field update of state estimation.

### 4.3. Strategy two: Machine Learning experiment to obtain a DA state estimation surrogate tool

The analysis of the DA performance, when applied to the turbulent plane channel flow, highlighted the importance of the flow update due to state estimation in order to obtain accurate results. However, it was also shown that this task can be performed only when observation is available, therefore
<figure><img src="image_10.png" /><figcaption>Figure 9: Isocontours of Q-criterion calculated for $y/h = 0.18$ ($\Delta y^* \approx 100$). The outcomes are presented for simulations: (a) DNS-IBM-ML$_{p.o.}$, (b) DNS-IBM-DA$_{p.o.}$, (c) DNS-IBM, and (d) R-DNS-BF.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: FFT calculated by sampling the fluctuating velocity field $u'$ located at (first row) $\Delta y^* \approx 30$ and (second row) $\Delta y^* \approx 56$. Results are shown for the simulations (---) DNS-IBM-ML<sub>p,o.</sub>, (—) DNS-IBM-DA<sub>p,o.</sub>, (···) DNS-IBM-CF, (—) DNS-BF, and (—) R-DNS-BF. With a dashed black line, we represent the -5/3 slope, which corresponds to the trend in the inertial range.</figcaption></figure>
precluding systematic application, in particular if data is not available for
a long time window. Here, RFR techniques are used to obtain a model
mimicking the flow correction of state estimation. The training of such a
model is performed by sampling flow information on the fly while applying the
DA algorithm. This case is significantly more complex than the generation
of the black-box IBM. First, from a physical perspective, the model itself is
more intricate since it involves the representation of the Kalman gain matrix
$K_k$ from the DA algorithm, which captures the optimal nonlinear correlation
between the state and high-fidelity observations. This task is essential for
an accurate representation of nonlinear dynamics. However, an important
simplification is that the test case investigated here is statistically steady.
This implies that once the EnKF dissipates initial discrepancies due to the
selected prior state, the action of the DA tool is statistically the same over
the observation window. In particular, the optimised coefficients for $D_{ij}$ are
stationary in time, which simplifies the task of setting the boundaries for the
observation window used for the ML training.

Figure 11 shows the scheme of the tasks performed. CONES [12] is employed to establish a connection between CFD, DA, and ML tools for the case of investigation, which is the turbulent plane channel flow previously investigated. The black-box IBM<sub>RFR</sub> $\hat{f}_P$ derived in §4.3, is included within the Navier–Stokes equations through the PISO algorithm (see Appendix B). The second ML algorithm aims to create a generative model mimicking the correction via state estimation. This tool is referred to as RFR for state estimation (SE<sub>RFR</sub>).

The input and output features used for the training of SE<sub>RFR</sub> are now
discussed. The feature matrix and the dependent values are defined, re-
spectively, as $\mathcal{X} = (u_{x,k}^{f+}, u_{y,k}^{f+}, u_{z,k}^{f+}) = (u_{x,k}^{f}/u_{\tau_k}, u_{y,k}^{f}/u_{\tau_k}, u_{z,k}^{f}/u_{\tau_k})$, and $\mathcal{Y} = (\Delta u_{x,k}^{+}, \Delta u_{y,k}^{+}, \Delta u_{z,k}^{+}) = ((u_{x,k}^{a} - u_{x,k}^{f})/u_{\tau_k}, (u_{y,k}^{a} - u_{y,k}^{f})/u_{\tau_k}, (u_{z,k}^{f} - u_{z,k}^{f})/u_{\tau_k})$,
where the instantaneous value of $u_{\tau_k}$ is employed for normalisation in the fea-
tures. This step clearly shows the contribution of the black-box model as an
additive correction term, taking the forecast as input and providing as out-
put the difference with the field produced by the analysis phase. $k$ stands for
every analysis phase where data is sampled in the observation window. In
this setup, data from one of the EnKF members $n_e$ is collected in the form of
dimensionless forecast velocities $u_{k,n_e}^{f+}$ and dimensionless analysis velocities
$u_{k,n_e}^{a+}$ at the 10 near-wall layers located either in the solid $\Omega_b$ or interface $\Sigma_b$
regions (5 at the top wall and 5 and the bottom wall) over a time period of $1.8 t_A$, which corresponds to $k_{max} = 15$ consecutive analysis phases of the EnKF. The total amount of data is $m_T = 1228800$ sampling data points, selected to provide significant representation while optimising computational efficiency and minimising RAM usage, requiring only 7.5 GB of memory. The same approach is applied for the training of the IBM$_{RFR}$ model, where the features are also normalised with respect to the instantaneous $u_{\tau_k}$. Data is collected at the same time instances, resulting in the same number of sampling data points as in the SE$_{RFR}$ model while reducing RAM usage to just 2.1 GB of memory. In both models, the hyperparameters during training are consistent with those discussed in §4.3. One can see in Table 3 that the metrics $R^2$ and NRMSE for all predictions show an excellent rate of convergence.

<figure><img src="image_12.png" /><figcaption>Figure 11: Overview of the DA code CONES employed for state estimation and parameter optimisation, highlighting the components that are modelled by either IBM$_{RFR}$ or SE$_{RFR}$.</figcaption></figure>

Figure 12 shows the results of the simulation DNS-IBM-ML$_{s.e.}$, which integrates both RFR models (IBM$_{RFR}$ and SE$_{RFR}$) under the same conditions as the training dataset. The first-order statistics now perfectly align with the homologous DNS-IBM-DA$_{s.e.}$ and the R-DNS-BF runs, with a discrepancy of 2.71% for the value of $u_{\tau}$ compared to the reference fine-grained simulation
<table><caption>Table 3: Metrics computed for all dependent variables estimated via either IBM<sub>RFR</sub> or SE<sub>RFR</sub>.</caption><thead><tr><th></th><th>f&#x0302;<sub>P<sub>x</sub></sub></th><th>&Delta;u&#x0302;<sub>x</sub></th><th>&Delta;u&#x0302;<sub>y</sub></th><th>&Delta;u&#x0302;<sub>z</sub></th></tr></thead><tbody><tr><td>R<sup>2</sup></td><td>99.91%</td><td>99.98%</td><td>99.97%</td><td>99.98%</td></tr><tr><td>NRMSE</td><td>15.59%</td><td>7.73%</td><td>3.16%</td><td>1.49%</td></tr></tbody></table>

and computational costs $C.C^* = 8.78$, matching the order of magnitude of the coarse-grained body-fitted simulation (DNS-BF) as detailed in Table 1. This demonstrates the critical role of state estimation in unsteady simulations in obtaining accurate field predictions. Additionally, all Reynolds stress tensor components, except for $\langle u'_y u'_y \rangle^+$, show strong consistency with the curves obtained using DA.

Instantaneous features are shown in Figure 13 for the spatial structure of the streaks in a plane parallel to the wall. The R-DNS-BF simulation exhibits the shortest streaks and the smallest velocity fluctuations, whereas the DNS-IBM simulation produces the longest streaks with the highest velocity fluctuations. The DNS-IBM-DA<sub>s.e.</sub> and DNS-IBM-ML<sub>s.e.</sub> simulations fall between these extremes, even if the mesh resolution is identical to the one for the run DNS-IBM. For brevity, the spectral analysis is omitted, as it reveals no significant deviations from the curves obtained for the DNS-IBM-ML<sub>p.o.</sub> simulation.

These results highlight the potential for simultaneously applying DA and ML algorithms on the fly, facilitating advanced applications in complex flow scenarios. For instance, in DTs, real-world physical systems could provide sparse spatial and temporal measurements obtained experimentally and integrate them into a coarse-grained/low-fidelity simulation by means of a sequential DA algorithm to produce an optimised continuous state. Subsequently, an ML tool could leverage this extensive, high-fidelity dataset to accurately refine operational parameters in real time, enhancing system performance and preventing premature failure, especially under extreme conditions.

However, one could argue about how the model behaves under circumstances that deviate from those used for sampling the training dataset. This important point is discussed and investigated in §5.
<figure><img src="image_13.png" /><figcaption>Figure 12: Comparison of the main statistical moments of the velocity field. Results are shown for simulations (---) DNS-IBM-ML<sub>s.e.</sub>, (···) DNS-IBM-DA<sub>s.e.</sub>, and (—) R-DNS-BF.</figcaption></figure>

5. Beyond training conditions: different grid refinement and $Re_τ$ numbers

In this section, two open questions about the application of DA and ML tools outside the training conditions are studied. Firstly, we investigate the impact of altering mesh refinement to achieve a wall resolution of $y^+ \approx 1$, corresponding to spatial interpolation scenarios. Secondly, we explore the effect of varying $Re_τ$ and, consequently, modifying the near-wall physics.

For the test cases with different grid refinements, we investigate two configurations with higher resolution near the wall. In the first arrangement, labelled as IBM-ML-grid 1, the grid presents fewer mesh elements than the one with the training conditions (IBM-DNS-ML<sub>s.e.</sub>) by coarsening the cells in the streamwise and spanwise directions by a factor of two, while simultaneously refining the wall-normal direction by the same factor. In contrast, the second configuration, denoted as IBM-ML-grid 2, contains more mesh elements by refining only the wall-normal direction by a factor of two. Table 4 provides a summary of these two complementary simulations, where it is important to note that, in these two runs, all quantities $\Delta(\cdot)^*$ are normalised with respect to $u_τ$ from R-DNS-BF.
<figure><img src="image_14.png" /><figcaption>Figure 13: Visualisation of the near-wall streaks for $y/h = 0.03$ ($\Delta y^* \approx 15$). Results are shown for simulations (a) DNS-IBM-ML$_{s.e.}$, (b) DNS-IBM-DA$_{s.e.}$, (c) DNS-IBM, and (d) R-DNS-BF.</figcaption></figure>
<table><caption>Table 4: Overview of supplementary simulations conducted to validate the ML model. $\Delta u_{\tau} = \frac{|u_{\tau} - u_{\tau_{ref}}|}{u_{\tau_{ref}}}$ is normalised by the reference friction velocity $u_{\tau_{ref}}$ at each Reynolds number. Reference values are taken from R-DNS-BF for test cases with $Re_{\tau} \approx 550$, from Moser et al. [66] for $Re_{\tau} \approx 395$, and from Hoyas et al. [60] for $Re_{\tau} \approx 950$.</caption><thead><tr><th></th><th>N<sub>x</sub> × N<sub>y</sub> × N<sub>z</sub></th><th>Δx<sup>*</sup></th><th>Δz<sup>*</sup></th><th>Δy<sup>*</sup><sub>min</sub></th><th>Δy<sup>*</sup><sub>max</sub></th><th>L<sub>x</sub>/h</th><th>L<sub>z</sub>/h</th><th>u<sub>τ</sub></th><th>Δu<sub>τ</sub></th></tr></thead><tbody><tr><td>DNS-IBM-ML<sub>s.e.</sub></td><td>128 × 64 × 64</td><td>39.5</td><td>26.3</td><td>1.3</td><td>67.6</td><td>3π</td><td>π</td><td>0.0493</td><td>0.0271</td></tr><tr><td>IBM-ML-grid 1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.0464</td><td>0.0333</td></tr><tr><td>IBM-DA-grid 1</td><td>64 × 128 × 32</td><td>79.0</td><td>52.7</td><td>1</td><td>30</td><td>3π</td><td>π</td><td>0.0454</td><td>0.0542</td></tr><tr><td>IBM-grid 1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.0312</td><td>0.3500</td></tr><tr><td>BF-grid 1</td><td>64 × 122 × 32</td><td>79.0</td><td>52.7</td><td>0.82</td><td>24.7</td><td>3π</td><td>π</td><td>0.0356</td><td>0.2583</td></tr><tr><td>IBM-ML-grid 2</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.0483</td><td>0.0063</td></tr><tr><td>IBM-DA-grid 2</td><td>128 × 128 × 64</td><td>39.5</td><td>26.3</td><td>1</td><td>30</td><td>3π</td><td>π</td><td>0.0458</td><td>0.0458</td></tr><tr><td>IBM-grid 2</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.0329</td><td>0.3146</td></tr><tr><td>BF-grid 2</td><td>128 × 122 × 64</td><td>39.5</td><td>26.3</td><td>0.82</td><td>24.7</td><td>3π</td><td>π</td><td>0.0361</td><td>0.2479</td></tr><tr><td>ML-ReTau395</td><td>256 × 64 × 128</td><td>30.6</td><td>20.4</td><td>1.01</td><td>52.6</td><td>6π</td><td>2π</td><td>0.0496</td><td>0.0236</td></tr><tr><td>IBM-ReTau395</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.0392</td><td>0.2283</td></tr><tr><td>BF-ReTau395</td><td>256 × 58 × 128</td><td>30.6</td><td>20.4</td><td>0.84</td><td>43.9</td><td>6π</td><td>2π</td><td>0.0440</td><td>0.1339</td></tr><tr><td>Moser et al. [66]</td><td>384 × 257 × 384</td><td>9.7</td><td>4.8</td><td>0.03</td><td>7.2</td><td>2π</td><td>π</td><td>0.0508</td><td>—</td></tr><tr><td>ML-ReTau950</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.0454</td><td>0.0005</td></tr><tr><td>IBM-ReTau950</td><td>64 × 64 × 32</td><td>83.7</td><td>55.8</td><td>2.8</td><td>143.8</td><td>3π/2</td><td>π/2</td><td>0.0312</td><td>0.2048</td></tr><tr><td>BF-ReTau950</td><td>64 × 58 × 32</td><td>83.7</td><td>55.8</td><td>2.3</td><td>120.1</td><td>3π/2</td><td>π/2</td><td>0.0361</td><td>0.3128</td></tr><tr><td>Hoyas et al. [60]</td><td>3072 × 385 × 2304</td><td>7.6</td><td>3.8</td><td>0.03</td><td>7.6</td><td>8π</td><td>3π</td><td>0.0454</td><td>—</td></tr></tbody></table>
These simulations are compared with some runs featuring the same grid refinement for a classical penalty IBM (IBM-grid 1 and IBM-grid 2, respectively), body-fitted simulations (BF-grid 1 and BF-grid 2, respectively), and the combination of parametric optimisation and state estimation by DA using the methodology already explained in §2.2 (IBM-DA-grid 1 and IBM-DA-grid 2, respectively). The prior for the latter simulations is based on the field at $300 t_A$ and the optimised parameters $D_{ij}$ from the DNS-IBM-DA$_{s.e.}$ simulation. A normal distribution with a standard deviation of 1% is applied to these parameters $D_{ij}$ in each of the 40 ensemble members to introduce variability into the system. Furthermore, multiplicative inflation $\lambda = 1.01$ is set to the model parameters after each analysis phase to further increase the variability of the ensemble, and clipping is performed to the mesh elements near the wall located either in the body $\Omega_b$, fluid-solid interface $\Sigma_b$ regions or the viscous sublayer of the fluid, i.e., $y^+ < 5$ (see Appendix A for the definitions and [48] for the practical implementations). It is insightful to compare the RFR predictions with those from DA, considering a priori the latter as the best-case scenario for the ML tool. This is because the data used to train the RFR model is derived from the enhanced predictions generated by DA.

The R-DNS-BF run is again used as a reference for comparison. In Figure 14, the main statistical moments of the first investigation are shown, revealing that the simulations employing DA and ML, i.e., IBM-DA-grid 1 and IBM-ML-grid 1, respectively, achieve an extremely close overall agreement with the reference fine simulation, especially for the first-order statistics. The RFR model built with data using the original grid refinement demonstrates a good agreement with the DA results for all Reynolds stress tensor components, in particular for $\langle u_x' u_x' \rangle^+$ and $\langle u_x' u_y' \rangle^+$. While slight degradation is observed in $\langle u_y' u_y' \rangle^+$ compared to DA, it is still significantly more accurate than the classical body-fitted simulation (BF-grid 1) and even more so than the classical penalisation IBM (IBM-grid 1). Furthermore, $\langle u_z' u_z' \rangle^+$ is closer to the reference simulation curve than the one obtained via DA. The accuracy of the prediction is also observed for the friction velocity $u_\tau$. The ML model can predict its value with a discrepancy of 3.33% against the value estimated by the reference body-fitted simulation.

For the second configuration, the main statistical moments represented in Figure 15 demonstrate that the grid resolution in the streamwise and spanwise directions is crucial for accurately predicting flow statistics. The DA and ML runs (IBM-DA-grid 2 and IBM-ML-grid 2, respectively) yield results closely aligned with the R-DNS-BF benchmark, with the ML sim-
<figure><img src="image_15.png" /><figcaption>Figure 14: Comparison of the main statistical moments of the velocity field. Results are shown for simulations (---) IBM-ML-grid 1, (—) IBM-DA-grid 1, (—) BF-grid 1, (···) IBM-grid 1, and (—) R-DNS-BF.</figcaption></figure>

ulation nearly overlapping the reference curves when predicting the first-order statistics and $\langle u_x' u_y' \rangle^+$. Also, one can see that, compared with the DA methodology, the ML simulation performs similarly or slightly better across all the Reynolds stress tensor components. As observed, the ML run outperforms the homologous BF-grid 2 simulation except for $\langle u_y' u_y' \rangle^+$. Significant improvements are also evident when compared to the IBM-grid 2 run, although the latter predicts the Reynolds stress tensor components $\langle u_y' u_y' \rangle^+$ and $\langle u_z' u_z' \rangle^+$ more accurately. The prediction of the wall shear stress is also excellent, with only a discrepancy in the $u_\tau$ value of 0.63% against the R-DNS-BF run.

The second analysis using the ML model comes at its extension and validation across different $Re_\tau$ values. In particular, two conditions for the Reynolds number, namely $Re_\tau = 395$ and $Re_\tau = 950$, are investigated with the same grid resolution as in DNS-IBM-ML<sub>s.e.</sub>. These configurations are selected, first of all, because high-fidelity databases are available for these Reynolds numbers, allowing for accurate normalisation of the quantities $\Delta(\cdot)^*$ in Table 4 against $u_\tau$ from the reference data. Additionally, the different near-wall physics at these $Re_\tau$ provide valuable insights. In comparison
<figure><img src="image_16.png" /><figcaption>Figure 15: Comparison of the main statistical moments of the velocity field. Results are shown for simulations (---) ML-grid 2, (—) IBM-DA-grid 2, (—) BF-grid 2, (···) IBM-grid 2, and (—) R-DNS-BF.</figcaption></figure>

with the R-DNS-BF run with $Re_τ \approx 550$, the configuration with $Re_τ \approx 395$ exhibits milder velocity gradients in the proximity of the wall. On the contrary, a configuration with $Re_τ \approx 950$ presents a significantly thinner region with high-velocity gradients, and the influence of the outer flow on the near-wall region becomes stronger, increasing the complexity of the flow. Also, to ensure a comparable number of flow structures in these additional simulations, the size of the computational domain is halved in the streamwise and spanwise directions for the configurations with $Re_τ \approx 395$, and doubled for the simulations with $Re_τ \approx 950$. The bulk velocity and kinematic viscosity $ν$ are recalculated to ensure that the centreline velocity $U_c$ remains unchanged.

Given the differences in the physics for these cases compared to the original ML model and the lack of instantaneous high-fidelity data to be observed, DA runs were not performed for this analysis. Even if data were available, it is highly possible that the new DA-optimised parametric configuration would exhibit differences from the results obtained for the run DNS-IBM-DA$_{s.e.}$, establishing difficulties for the comparisons between the tools. The DNS database by Moser et al. [66] serves as the reference for the coarse-grained simulations with $Re_τ \approx 395$. One can see in Figure 16
<figure><img src="image_17.png" /><figcaption>Figure 16: Comparison of the main statistical moments of the velocity field for Reτ = 395. Results are shown for simulations (·-·) ML-ReTau395, (—) BF-ReTau395, (···) IBM-ReTau395, and (—) Moser et al. [66].</figcaption></figure>

that the ML model, alluded to as ML-ReTau395, matches the reference data when predicting the first-order statistics and shows a relatively close match in predicting $\langle u'_z u'_z \rangle^+$ and $\langle u'_x u'_y \rangle^+$. Overall, the ML-ReTau395 run predicts the flow features more accurately than their homologous body-fitted (BF-ReTau395) and classical penalisation IBM (IBM-ReTau395) runs with the same grid resolution. Besides, the prediction of the $u_τ$ is notably accurate, with only a slight discrepancy of 2.36% with respect to the DNS data. This suggests that the ML model can effectively predict wall-bounded turbulent boundary layers when the wall resolution is better than that used for the training. This finding is expected considering that the grids used for the case at Reτ ≈ 550 were already able to resolve the viscous sublayer of the flow. Therefore, no modification in the physical behaviour is observed.

At last, the case for Reτ ≈ 950 is considered. The DNS database by Hoyas et al. [60] is used for comparison, with the corresponding curves represented in Figure 17. This case highlights the limitations of the ML model (ML-ReTau950) when applied to flow conditions that differ significantly from the training environment. Although the predicted first-order statistics and $u_τ$ show a reasonable agreement with the reference data, second-order statis-
tics exhibit poorer resolution compared to the original body-fitted and IBM
approaches, particularly for $\langle u'_y u'_y \rangle^+$ and $\langle u'_z u'_z \rangle^+$. The main problem here
is associated with previous considerations for the case $Re_τ ≈ 395$. The ML
model is trained using data from a wall-resolved simulation, while, in this
case, the first grid point is already outside the training conditions spatially.
Therefore, the physics learnt during the training does not match what the
algorithm needs to replicate in this case. To accurately capture the flow fea-
tures under these conditions, the ML model would either need to incorporate
wall modelling physics as already performed with physics-informed neural
networks (PINNs) for LES models [67, 68], quantify the uncertainty in the
coarse-grained DNS simulations to bound all physically realisable states, like
in physics-informed RFRs for RANS simulations [55, 56, 57]), or be trained
on a broader database that includes higher $y^+$ values. If the ML algorithm
could be trained over time with DA fields obtained with different $Re_τ$, such
as $Re_τ = 395$ and $Re_τ = 1200$, that would expose the model to different
conditions in terms of $u_τ$. That would also imply that the ML tool would
observe different values of $y^+$, for example $y^+ ∈ [0.5, 3]$, therefore enriching
its predictive capabilities. Such variations in $Re_τ$ could be associated with
a change of regime during the DA observation window, if the latter is long
enough. In this case, if one wanted to use the trained ML tool for interme-
diate $Re_τ$ such as $Re_τ = 550$ or $Re_τ = 950$, the model would probably be
able to provide a satisfying correction of the numerical field obtained. In
summary, the present analysis indicates that DA can provide more complete
snapshots for ML training, therefore improving the robustness of the tools
thus obtained. However, existing limitations in ML extrapolation capabilities
cannot be completely bypassed, unless a relatively small quantity of sensor
data is available to be augmented via DA approaches.

# 6. Conclusions

The research activities in this work propose using Machine Learning tools
to learn and mimic the state estimation corrections performed in the anal-
ysis phase of DA strategies. This research advancement opens perspectives
of more systematic application of such state corrections, which are limited
to the data availability in classical DA procedures. The application here is
performed by combining the Ensemble Kalman Filter, a sequential Data As-
similation algorithm, and Random Forest Regression, an ensemble-learning
Machine Learning method. The generated black-box models consist of a cor-
<figure><img src="image_18.png" /><figcaption>Figure 17: Comparison of the main statistical moments of the velocity field. Results are shown for simulations (—) ML-ReTau950, (—) BF-ReTau950, (···) IBM-ReTau950, and (—) Hoyas *et al.* [60].</figcaption></figure>

rection term accounting for DA's state estimation as well as a surrogate IBM. For the latter, the training data is taken from a classical penalisation IBM, whose free coefficients are optimised by the EnKF. The strategy is assessed by investigating the turbulent plane channel flow test case. The analysis of the flow statistics, including the friction velocity $u_τ$, the mean streamwise velocity $\langle U \rangle^+$ and the Reynolds stress tensor components $\langle u_i' u_j' \rangle^+$, indicates that the black-box models are able to obtain a predictive accuracy of the same level as the DA approaches. However, the analysis of the spectra of the fluctuating component of the velocity field shows some degradation of the results, whose accuracy is in line with the prior models used in the investigation. When applied to configurations outside of the training range, such as those due to variations in wall grid resolution and $Re_τ$, the ML black-box tools show good features of robustness, computational efficiency, and adaptability to similar physical configurations. However, limitations are clearly observed when the models are used to extrapolate physical information which is not available in the training set.

The hybrid approach combining DA with ML presented in this work shows significant potential for applications involving parametric optimisa-
tion and state estimation via on-the-fly training. The team's current research explores this DA-ML combination in statistically unsteady cases, such as an oscillatory cylinder with horizontal sinusoidal motion. This strategy holds promise for delivering fast, accurate, and interpretable results in various applications, suggesting a valuable pathway for future work in the field where both methodologies could be employed with streaming data in real industrial configurations.

## Acknowledgments

The authors would like to thank Prof. Nicolas Hascoët for his initial insights into the Machine Learning models used in this study. This project was provided with computer and storage resources by GENCI at TGCC, thanks to the grant 2023-A0142A01741 on the supercomputer Joliot Curie's SKL partition. This research was conducted in the framework of the ANR-JCJC-2021 IWP-IBM-DA project.

## Appendix A. Ensemble Kalman Filter algorithm

The Ensemble Kalman Filter (EnKF) is a Monte-Carlo approach of the Kalman Filter (KF) [69], where both the physical system state **u** (typically the velocity field discretised over the computational domain in CFD applications) and the covariance error matrix **P**—which quantifies the discrepancy **e** between the system state predicted by a dynamic model *M*(*θ*) defined by some parameters *θ*, and its true state (i.e., **P** = E[**e**(*e*)<sup>T</sup>])—are propagated in time. This propagation is achieved by defining an ensemble of realisations called ensemble members. Each ensemble member has distinct state vectors **u**<sub>*i*</sub> (*i* ∈ [1, *N*<sub>*e*</sub>], where *N*<sub>*e*</sub> is the number of ensemble members) to represent the system state. The average of these vectors $\bar{\mathbf{u}}$ is assumed to correspond to the true state, while their variability reflects the error covariance matrix **P**. To enhance performance, the *extended-state* EnKF approach [1] can be employed. In this method, certain model parameters *θ* are also propagated alongside the system state to reduce model bias. Hence, the system state **u}' = [**u** *θ*]<sup>T</sup>.

The EnKF operates in a two-step procedure. The first step corresponds to a time advancement in the numerical model $\mathcal{M}(\theta)$, known as forecast, represented as $\bullet^f$. By assuming that this progression is effectuated from time
k − 1 to time k , the system state predicted by the model can be estimated as:

$$
\boldsymbol{u}_{i,k}^{f} = \mathcal{M}(\theta_{i,k:k-1}) \boldsymbol{u}_{i,k-1} \qquad (A.1)
$$

If observation is available in time k , an analysis phase, denoted by · $^{a}$, is
performed, where the EnKF algorithm updates the system state and the
model to optimally combine the dynamics of both the model and observa-
tions, effectively minimising the covariance error matrix of the updated state
P k . However, rather than explicitly computing this matrix—which would
be computationally prohibitive due to its size of [ n DF , n DF ], where n DF is
the system’s degrees of freedom and equal to 3 times the number of mesh
elements for 3D simulations in CFD applications—the EnKF uses an al-
ternative approach. Instead, the algorithm defines an anomaly matrix X ,
which measures the deviation of the ensemble members from the mean u . If
the ensemble is representative enough, and there is statistical independence
among the ensemble members, the forecast covariance matrix P k can be
approximated by:

$$
\mathbf{P}_k^f \approx \mathbf{X}_k^f (\mathbf{X}_k^f)^T \qquad (A.2)
$$

where the anomaly matrix **X**<sub>k</sub><sup>f</sup> has size [n<sub>DF</sub>, N<sub>e</sub>] and is estimated from the
state vectors **u**<sub>i</sub> and the mean **u**&#773; as:

$$
\bar{\mathbf{u}}_k^f = \frac{\sum_{i}^{N_e} \mathbf{u}_{i,k}^f}{N_e} \qquad (A.3)
$$

$$
\left[ \mathbf{X}_k^f \right]_i = \frac{\mathbf{u}_{i,k}^f - \bar{\mathbf{u}}_k^f}{\sqrt{N_e - 1}} \qquad (A.4)
$$

We require an ensemble of n<sub>o</sub> observations to ensure dimensional consistency during algebraic operations. The resulting observation matrix **Y**<sub>k</sub>, which has dimensions [n<sub>o</sub>, N<sub>e</sub>], is derived from perturbing the observation vector **y**<sub>k</sub> of size [n<sub>o</sub>, 1]. This perturbation is modelled using Gaussian noise ε, where the mean corresponds to the actual observation values in **y**<sub>k</sub>. The standard deviation reflects the observation uncertainty, encapsulated in a diagonal matrix denoted as **R**<sub>k</sub> with dimensions [n<sub>o</sub>, n<sub>o</sub>] and called the observation covariance matrix. The diagonal structure of this matrix assumes that uncertainties are Gaussian, unbiased, and independent. For an infinite
ensemble size, $\lim_{N_e \to +\infty} \mathbb{E}[(\epsilon - \mathbb{E}[\epsilon])(\epsilon - \mathbb{E}[\epsilon])^T] = \mathbf{R}$, which becomes time-independent. Thus, each entry of the observation matrix can be expressed as $[\mathbf{Y}_k]_i = \mathbf{y}_k + \epsilon_i$, where $\epsilon_i \sim \mathcal{N}(0, \mathbf{R})$.

The EnKF also calls for the use of a sampling matrix $\mathcal{H}(\mathbf{u}_k^f)$, where $\mathcal{H}$ acts as an interpolator that projects the forecast system state into the observation framework. Similar to the system's state $\mathbf{u}_k^f$ in Equation A.4, this projection is composed by a matrix $\mathbf{S}_k^f$ where each column represents a normalised anomaly estimated from the mean $\overline{\mathcal{H}(\mathbf{u}_k^f)}$ as:

$$ \overline{\mathcal{H}(\mathbf{u}_k^f)} = \frac{\sum_{i}^{N_e} [\mathcal{H}(\mathbf{u}_k^f)]_i}{N_e} \qquad (A.5) $$

$$ [\mathbf{S}_k^f]_i = \frac{[\mathcal{H}(\mathbf{u}_k^f)]_i - \overline{\mathcal{H}(\mathbf{u}_k^f)}}{\sqrt{N_e - 1}} \qquad (A.6) $$

The Kalman gain matrix $\mathbf{K}_k$ describing the most optimal correlation between the state and the observations can be defined as expressed in Equation A.7 [70, 71]. In our context, we apply an element-wise multiplication with a covariance localisation matrix $\mathbf{L}$. This matrix is uniform for static grids or observation sensors, with coefficients $L_{ij} \in [0, 1]$. This approach enhances the significance of coefficients associated with closer distances $\Delta\mathbf{r}_{ij} = (\Delta x_{ij}, \Delta y_{ij}, \Delta z_{ij})$ between the centre of each mesh element and the observation location, compared to those further away. The parameters $(\eta_x, \eta_y, \eta_z)$ act as tuning factors, balancing the spatial scale of the system's physical representation with the smoothing needed to mitigate discontinuities at the boundaries of the updated region.

$$ \mathbf{K}_k = \mathbf{L} \odot \mathbf{X}_k^f (\mathbf{S}_k^f)^T \left[ \mathbf{S}_k^f (\mathbf{S}_k^f)^T + \mathbf{R} \right]^{-1} \qquad (A.7) $$

$$ L_{ij} = e^{-\left[ \left( \frac{\Delta x_{ij}}{\eta_x} \right)^2 + \left( \frac{\Delta y_{ij}}{\eta_y} \right)^2 + \left( \frac{\Delta z_{ij}}{\eta_z} \right)^2 \right]} \qquad (A.8) $$

All the elements described are crucial for calculating the updated system's state $\mathbf{u}_k^a$, which is subsequently used to advance the numerical model $\mathcal{M}(\theta_{i,k+1:k})$ in the following iteration $k + 1$. Typically, multiplicative inflation $\lambda_k$ is applied to perturb the updated state, as the filter tends to converge on local
optima that are closer to the forecast state rather than the desired global optimum.

$$ \mathbf{u}_k^a = \mathbf{u}_k^f + \mathbf{K}_k (\mathbf{Y}_k - \mathcal{H}(\mathbf{u}_k^f)) \quad (A.9) $$

$$ \mathbf{u}_k^a = \bar{\mathbf{u}}_k^a + \lambda_k (\mathbf{u}_k^a - \bar{\mathbf{u}}_k^a) \quad (A.10) $$

The step-by-step implementation of the EnKF is detailed in the Algorithm 1.

<table><caption>**Algorithm 1:** Scheme of the EnKF used in the present study.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Advancement in time of the state vectors:<br>u<sub>i,k</sub><sup>f</sup> = M<sub>k:k-1</sub> u<sub>i,k-1</sub></td></tr><tr><td>2</td><td>Generation of an observation matrix from the confidence level given to the observation data:<br>Y<sub>i,k</sub> = y<sub>k</sub> + ε<sub>i</sub>, with ε<sub>i</sub> ~ N(0, R)</td></tr><tr><td>3</td><td>Estimation of the ensemble means (system state and projection matrix):<br>u&#773;<sub>k</sub><sup>f</sup> = 1/N<sub>e</sub> &sum;<sub>i=1</sub><sup>N<sub>e</sub></sup> u<sub>i,k</sub><sup>f</sup>, H(u&#773;<sub>k</sub><sup>f</sup>) = 1/N<sub>e</sub> &sum;<sub>i=1</sub><sup>N<sub>e</sub></sup> H(u<sub>i,k</sub><sup>f</sup>)</td></tr><tr><td>4</td><td>Computation of the anomaly matrices (system state and projection matrix):<br>X<sub>k</sub> = (u<sub>i,k</sub> - u&#773;<sub>k</sub>) / &radic;(N<sub>e</sub>-1), S<sub>k</sub> = (H(u<sub>i,k</sub><sup>f</sup>) - H(u&#773;<sub>k</sub><sup>f</sup>)) / &radic;(N<sub>e</sub>-1)</td></tr><tr><td>5</td><td>Calculation of the Kalman gain (with covariance localisation):<br>K<sub>k</sub> = L &odot; X<sub>k</sub><sup>f</sup>(S<sub>k</sub>)<sup>T</sup> [S<sub>k</sub>(S<sub>k</sub>)<sup>T</sup> + R]<sup>-1</sup></td></tr><tr><td>6</td><td>Update of the state matrix:<br>u<sub>i,k</sub><sup>a</sup> = u<sub>i,k</sub><sup>f</sup> + K<sub>k</sub> (Y<sub>i,k</sub> - H(u<sub>i,k</sub><sup>f</sup>))</td></tr><tr><td>7</td><td>Application of (deterministic) covariance inflation:<br>u<sub>i,k</sub><sup>a</sup> = u&#773;<sub>k</sub><sup>a</sup> + λ<sub>k</sub> (u<sub>i,k</sub><sup>a</sup> - u&#773;<sub>k</sub><sup>a</sup>)</td></tr></tbody></table>
## Appendix B. Random Forest Regression in PISO algorithm

The integration of the RFR model into the CFD solver involves modifying the PISO algorithm [45] for parametric optimisation and introducing an additional step for state estimation. During this process, the velocity and pressure fields are iteratively updated to satisfy the momentum and Poisson equations. The RFR requires an explicit estimation of the penalisation volume force $\hat{f}_P$ inside Equation 2. At time $t$, with a time step advancement of $\Delta t$, and considering $j \in [1, J]$ as a single iteration within the PISO loop, the following steps are executed:

1. Resolution of the momentum equation, from where $u_{t,0}$ is obtained.

$$A \mathbf{u}_{t,0} - \mathbf{b} - \hat{\mathbf{f}}_{P_{t-\Delta t}} = -\nabla p_{t-\Delta t} \quad (B.1)$$

2. Estimation of the dimensional forcing term $\hat{\mathbf{f}}_{P_{t,0}}$ from $\hat{\mathbf{f}}_{P_{t,0}}^{+/*}$ with the predicted velocity $\mathbf{u}_{t,0}$ by means of RFR. By calling the operator $\mathcal{F}_{\text{p.o.}}^r$ to define the prediction of a single decision tree $r$ and normalising with respect to the instantaneous $u_{\tau_{t,0}}$ or the statistical $u_{\tau_{\text{DA/p.o.}}}$ friction velocity:

$$\hat{\mathbf{f}}_{P_{t,0}}^{+/*} = \frac{1}{R} \sum_{r=1}^{R} \mathcal{F}_{\text{p.o.}}^{r} \left( \mathbf{u}_{t,0}^{+/*}, y_{t,0}^{+/*} \right) \quad (B.2)$$

3. Estimation of the pressure $p_{t,j}$ through the Poisson equation. $A$ is a scalar field calculated from $A$, whereas $T'(\mathbf{u}_{t,j-1})$ is a tensor field containing the discretised form of all the terms on the left side of Equation B.1. To account for the update of the forcing term $\hat{\mathbf{f}}_{P_{t,j-1}}$, an additional term is included.

$$\nabla \cdot \frac{1}{A} \nabla p_{t,j} = \nabla \cdot \left( \frac{T'(\mathbf{u}_{t,j-1})}{A} \right) + \nabla \cdot \left( \frac{\hat{\mathbf{f}}_{P_{t,j-1}}}{A} \right) \quad (B.3)$$

4. Update of the velocity field $\mathbf{u}_{t,j}$ to satisfy the zero-divergence condition.

$$\mathbf{u}_{t,j} = \frac{T'(\mathbf{u}_{t,j-1})}{A} - \frac{1}{A} \nabla p_{t,j} + \frac{1}{A} \hat{\mathbf{f}}_{P_{t,j-1}} \quad (B.4)$$

5. Re-estimation of the forcing term with the new velocity field $\mathbf{u}_{t,j}$.

$$\hat{\mathbf{f}}_{P_{t,j}}^{+/*} = \frac{1}{R} \sum_{r=1}^{R} \mathcal{F}_{\text{p.o.}}^{r} \left( \mathbf{u}_{t,j}^{+/*}, y_{t,j}^{+/*} \right) \quad (B.5)$$
Equations B.3 to B.5 are solved iteratively until reaching convergence. For state estimation at time t, when j = J, the velocity $\mathbf{u}_{t,j}$ corresponds to $\mathbf{u}_k^f$, which is equivalent to the forecast velocity field at the $k^{th}$ analysis phase for the EnKF algorithm. To refine its estimation $\mathbf{u}_k^a$ as done with the EnKF in the analysis step, it is necessary to incorporate another RFR model, denoted by the operator $\mathcal{F}_{\text{s.e}}^r$ for a single decision tree r, which is able to predict the correction given to the velocity $\Delta \hat{\mathbf{u}}_k^{+/*}$, so that $\hat{\mathbf{u}}_k^{a+/*} = \mathbf{u}_k^{f+/*} + \Delta \hat{\mathbf{u}}_k^{+/*}$. Then, the dimensional $\hat{\mathbf{u}}_k^a$ is estimated.

$$ \Delta \hat{\mathbf{u}}_k^{+/*} = \frac{1}{R} \sum_{r=1}^{R} \mathcal{F}_{\text{s.e}}^r (\mathbf{u}_k^{f+/*}) \quad (B.6) $$

## References

[1] M. Asch, M. Bocquet, M. Nodet, Data Assimilation: Methods and Applications, SIAM, 2016.

[2] G. Evensen, F. Vossepoel, P. J. van Leeuwen, Data Assimilation Fundamentals, Springer Textbooks in Earth Sciences, Geography and Environment, 2022.

[3] A. Burkov, The Hundred-Page Machine Learning Book, Andriy Burkov, 2019.

[4] S. Brunton, B. Noack, P. Koumoutsakos, Machine learning for fluid mechanics, Annual Review for Fluid Mechanics 52 (2020) 477–508.

[5] C. Gorlé, G. Iaccarino, A framework for epistemic uncertainty of turbulent scalar flux models for Reynolds-averaged Navier-Stokes simulations, Physics of Fluids 25 (2013) 055105.

[6] W. N. Edeling, P. Cinnella, R. P. Dwight, H. Bijl, Bayesian estimates of parameter variability in the $\mathcal{K}-\varepsilon$ turbulence model, Journal of Computational Physics 258 (2014) 73–94.

[7] L. Margheri, M. Meldi, M. V. Salvetti, P. Sagaut, Epistemic uncertainties in RANS model free coefficients, Computers & Fluids 102 (2014) 315–335.
[8] B. Tracey, K. Duraisamy, J. Alonso, A paradigm for data-driven predictive modeling using field inversion and machine learning, 53rd AIAA Aerospace Sciences Meeting Session: Turbulence Modeling I (2015) 2015-1287.

[9] J. L. Wu, H. Xiao, E. Paterson, Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework, Physical Review Fluids 7 (2018) 074602.

[10] P. Srinivasan, L. Guastoni, H. Azizpour, P. Schlatter, R. Vinuesa, Predictions of turbulent shear flows using deep neural networks, Physical Review Fluids 4 (2019) 054603.

[11] P. S. Volpiani, M. Meyer, L. Franceschini, J. Dandois, F. Renac, E. Martin, O. Marquet, D. Sipp, Machine learning-augmented turbulence modeling for rans simulations of massively separated flows, Physical Review Fluids 6 (2021) 064607.

[12] L. Villanueva, M. M. Valero, A. S. Glumac, M. Meldi, Augmented state estimation of urban settings using on-the-fly sequential data assimilation, Comp. Fluids 269 (2024) 106118.

[13] M. Meldi, D. Lucor, P. Sagaut, Is the Smagorinsky coefficient sensitive to uncertainty in the form of the energy spectrum?, Physics of Fluids 23 (2011) 125109.

[14] A. Vollant, G. Balarac, C. Corre, Subgrid-scale scalar flux modelling based on optimal estimation theory and machine-learning procedures, Journal of Turbulence 18 (9) (2017) 854-878.

[15] M. Meldi, Augmented Prediction of Turbulent Flows via Sequential Estimators: Sensitivity of State Estimation to Density of Time Sampling for Available Observation, Flow, Turbulence and Combustion 101 (2018) 389-412.

[16] P. Chandramouli, E. Memin, D. Heitz, 4D large scale variational data assimilation of a turbulent flow with a dynamics error model, Journal of Computational Physics 412 (2020) 109446.

[17] A. Lozano-Durán, J. Bae, Self-critical machine-learning wall-modeled les for external aerodynamics, Annual Research Briefs 2020 - (2020) 1-15.
[18] V. Mons, Y. Du, T. Zaki, Ensemble-variational assimilation of statistical data in large-eddy simulation, Physical Review Fluids 6 (2021) 104607.

[19] G. Moldovan, A. Mariotti, G. Lehnasch, L. Cordier, M. Salvetti, M. Meldi, Multigrid sequential data assimilation for the large-eddy simulation of a massively separated bluff-body flow, Computers and Fluids 281 (2024) 106385.

[20] L. Villanueva, K. Truffin, M. Meldi, Synchronization and optimization of large eddy simulation using an online ensemble kalman filter, International Journal of Heat and Fluid Flow 110 (2024) 109597. doi:10.1016/j.ijheatfluidflow.2024.109597. URL https://linkinghub.elsevier.com/retrieve/pii/S0142727X24003229

[21] J. Plogmann, O. Brenner, P. Jenny, Adjoint-based assimilation of sparse time-averaged velocity data in large-eddy simulations (2024). arXiv:eprint=2407.11746.

[22] R. Molinaro, J.-S. Singh, S. Catsoulis, C. Narayanan, D. Lakehal, Embedding data analytics and CFD into the digital twin concept, Comp. Fluids 214 (2021) 104759.

[23] D. Wagg, K. Worden, R. Barthorpe, P. Gardner, Digital twins: State-of-the-art and future directions for modelling and simulation in engineering dynamics applications, Amer. Soc. Mech. Eng. 6 (3) (2020) 030901.

[24] B. Yang, S. Yang, Z. Lv, F. Wang, T. Olofsson, Application of digital twins and metaverse in the field of fluid machinery pumps and fans: A review, Sensors 22 (2022) 9294.

[25] A. Sharma, E. Kohasih, J. Zhang, A. Brintrup, A. Calinescu, Digital twins: State of the art theory and practice, challenges, and open research questions, J. Ind. Inf. Int. 30 (2022) 100383.

[26] L. Schena, P. Marques, R. Poletti, S. Ahizi, J. V. der Berghe, M. Mendez, Reinforcement twinning: From digital twins to model-based reinforcement learning, Journal of Computational Physics 82 (2024) 102421.
[27] S. S. H. Boosari, Predicting the dynamic parameters of multiphase flow in cfd (dam-break simulation) using artificial intelligence-(cascading deployment), Fluids 4 (2019) 44.

[28] C. White, D. Ushizima, C. Farhat, Fast neural network predictions from constrained aerodynamics datasets, AIAA Scitech 2020 Forum (2020).

[29] D. Kochkov, J. Smith, A. Alieva, Q. Wang, M. Brenner, S. Hoyer, Machine learning-accelerated computational fluid dynamics, Applied Mathematics 118 (21) (2021).

[30] X. Fan, D. Akhare, J.-X. Wang, Neural differentiable modeling with diffusion-based super-resolution for two-dimensional spatiotemporal turbulence (2024). arXiv:2406.20047.

[31] G. Evensen, The ensemble kalman filter for combined state and parameter estimation, IEEE Control Systems Magazine (2009) 83-104.

[32] J. Brajard, A. Carrassi, M. Bocquet, L. Bertino, Combining data assimilation and machine learning to infer unresolved scale parametrization, Phil. Trans. R. Soc. A 379 (2021) 20200086.

[33] A. Farchi, P. Laloyaux, M. Bonavita, M. Bocquet, Using machine learning to correct model error in data assimilation and forecast applications, Quaterly Journal of the Royal Meteorological Society 147 (2021) 3067-3084.

[34] R. Arcucci, J. Zhu, S. Hu, Y.-K. Guo, Deep data assimilation: Integrating deep learning with data assimilation, Applied Sciences 11 (1114) (2021).

[35] R. Villiers, V. Mons, D. Sipp, E. Lamballais, M. Meldi, Enhancing unsteady reynolds-averaged navier-stokes modelling from sparse data through sequential data assimilation and machine learning, Flow, Turbulence and Combustion (2025). doi:10.1007/s10494-024-00623-6.

[36] Z. Y. Wang, W. W. Zhang, A unified method of data assimilation and turbulence modeling for separated flows at high reynolds numbers (2022). arXiv:2211.00601.
[37] M. Quattromini, M. A. Bucci, S. Cherubini, O. Semeraro, Active learning of data-assimilation closures using graph neural networks (2023). arXiv:2303.03806.

[38] Y. Ling, A. Lozano-Durán, Numerically consistent data-driven subgrid-scale model via data assimilation and machine learning, AIAA SciTech Forum (2025).

[39] L. Breiman, Random forests, Machine Learning 45 (2001) 5–32.

[40] P. Angot, C.-H. Bruneau, P. Fabric, A penalization method to take into account obstacles in incompressible viscous flows, Numer. Math. 81 (1999) 497–520.

[41] R. Verzicco, Immersed boundary methods: Historical perspective and future outlook, Annual Review of Fluid Mechanics 55 (2022) 129 – 155.

[42] J. Kou, A. H. de Mendoza, S. Joshi, S. L. Clainche, E. Ferrer, Eigen solution analysis of immersed boundary method based on volume penalization: applications to high-order schemes, Journal of Computational Physics 449 (2022) 10817.

[43] J. Ferziger, M. Peric, Computational Methods in Fluid Dynamics, New York : Springer-Verlag, 1996.

[44] C. Greenshields, H. Weller, Notes of Computational Fluid Dynamics: General Principles, CFD Direct Ltd, 2022.

[45] R. Issa, Solution of the implicitly discretised fluid flow equations by operator-splitting, Journal of Computational Physics 62 (1) (1986) 40–65. doi:https://doi.org/10.1016/0021-9991(86)90099-9.

[46] H. K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Second Edition, Pearson Prentice Hall, 2007.

[47] OpenFOAM - Field Operation And Manipulation, https://www.openfoam.com.

[48] M. M. Valero, M. Meldi, An immersed boundary method using online sequential data assimilation, Journal of Computational Physics 524 (3 2025). doi:10.1016/j.jcp.2024.113697.
[49] D. C. Wilcox, Turbulence Modeling for CFD, 3rd Edition, DCW Industries Inc., La Canada CA, 2006.

[50] H. Xiao, P. Cinnella, Quantification of Model Uncertainty in RANS Simulations: A Review, Progress in Aerospace Sciences 108 (2019) 1–31.

[51] P. Roßbach, Neural networks vs. random forests - does it always have to be deep learning?, frankfurt School Blog (2018).

[52] S. Wang, C. Aggarwal, H. Liu, Using a random forest to inspire a neural network and improving on it, in: Proceedings of the 2017 SIAM International Conference on Data Mining, 2017.

[53] M. de Zordo-Banliat, G. Dergham, X. Merle, P. Cinnella, Space-dependent turbulence model aggregation using machine learning, Journal of Computational Physics 497 (2024) 112628.

[54] S. Cherroud, X. Merle, P. Cinnella, X. Gloerfelt, Space-dependent aggregation of data-driven turbulence models, Journal of Computational Physics 527 (2025) 113793.

[55] M. L. Kaandorp, R. P. Dwight, Data-driven modelling of the reynolds stress tensor using random forests with invariance, Comp. Fluid. 202 (2020) 104497.

[56] J.-X. Wang, J.-L. Wu, H. Xiao, Physics-informed machine learning approach for reconstructing reynolds stress modeling discrepancies based on dns data, Physical Review Fluids 2 (2017) 034603.

[57] J. Heyse, A. Mishra, G. Iaccarino, Data driven physics constrained perturbations for turbulence model uncertainty estimation, in: AAAI Spring Symposium: MLPS, 2021.

[58] J. Kim, P. Moin, R. Moser, Turbulence statistics in fully developed channel flow at low reynolds number, Journal of Fluid Mechanics 177 (1987) 133–166.

[59] J. C. del Alamo, J. Jiménez, Spectra of the very large anisotropic scales in turbulent channels, Phys. Fluids 15 (6) (2003).
[60] S. Hoyas, J. Jiménez, Reynolds number effects on the reynolds-stress budgets in turbulent channels, Phys. Fluids 20 (2008) 101511.

[61] M. Bernardini, S. Pirozzoli, P. Orlandi, Velocity statistics in turbulent channel flow up to $re_τ = 4000$, Journal of Fluid Mechanics 742 (2014) 171–191.

[62] A. Cimarelli, E. D. Angelis, P. Schlatter, G. Brethouwer, A. Talamelli, C. Casciola, Sources and fluxes of scale energy in the overlap layer of wall turbulence, Journal of Fluid Mechanics 771 (2015) 407–423.

[63] D. E. King, Dlib-ml: A machine learning toolkit, Journal of Machine Learning Research 10 (2009) 1755–1758.

[64] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical recipes: The Art of Scientific Computing. Third Edition, Cambridge University Press, 2017.

[65] G. I. Taylor, The spectrum of turbulence, in: Proc. R. Soc. Lond. A, Vol. 164, 1938, pp. 476–490.

[66] R. D. Moser, J. Kim, N. N. Mansour, Direct numerical simulation of turbulent channel flow up to $re_τ = 590$, Physics of Fluids 11 (4) (1999) 943–945.

[67] S. Maejima, K. Tanino, S. Kawai, Physics-informed machine-learning solution to log-layer mismatch in wall-modeled large-eddy simulation, Physical Review Fluids 9 (2024) 084609.

[68] X. Yang, S. Zafar, J.-X. Wang, H. Xiao, Predictive large-eddy-simulation wall modeling via physics-informed neural networks, Physical Review Fluids 4 (2019) 034602.

[69] R. E. Kalman, A new approach to linear filtering and prediction problems, Journal Basic Eng. 82 (1960) 35–45.

[70] A. Carrassi, M. Bocquet, L. Bertino, G. Evensen, Data assimilation in the geosciences: An overview of methods, issues, and perspectives, WIREs Climate Change 9 (2018).
[71] I. Hoteit, D.-T. Pham, M. E. Gharamti, X. Luo, Mitigating observation perturbation sampling errors in the stochastic enkf, Monthly Weather Review 143 (2015) 2918–2936.