# An immersed boundary method using online sequential data assimilation

Miguel M Valero, Marcello Meldi

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

## To cite this version:

Miguel M Valero, Marcello Meldi. An immersed boundary method using online sequential data assimilation. Journal of Computational Physics, 2024, 524, pp.113697. 10.1016/j.jcp.2024.113697. hal-05488604

HAL Id: hal-05488604

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

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An Immersed Boundary Method using online sequential Data Assimilation

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

<sup>a</sup>Univ. Lille CNRS ONERA Arts et Metiers ParisTech Centrale Lille UMR 9014-LMFL- Laboratoire de Mecanique des fluides de Lille - Kampe de Feriet F-59000 Lille France

## Abstract

A data-driven strategy relying on the Ensemble Kalman Filter (EnKF) is here used to augment the accuracy of a continuous Immersed Boundary Method (IBM). The latter is a classical penalty method accounting for the presence of the immersed body via a volume source term, which is included in the Navier-Stokes equations. The model coefficients of the penalisation method, which are usually selected by the user, are optimised using the data-driven strategy. The parametric inference is governed by the physical knowledge of local and global features of the flow, such as the no-slip condition and the shear stress at the wall. The C++ library CONES (Coupling OpenFOAM with Numerical EnvironmentS) developed by the team is used to perform an online investigation, coupling on-the-fly data from synthetic sensors with results from an ensemble of coarse-grained numerical simulations. The analysis is performed for a classical test case, namely the turbulent plane channel flow with $Re_τ = 550$. The results, which are compared with high-fidelity Direct Numerical Simulation (DNS), show that the data-driven procedure exhibits remarkable accuracy despite the ensemble members' relatively low grid resolution. The findings present open perspectives of application in dynamic complex systems, such as digital twins.

*Keywords:* IBM, DA, EnKF, CONES, wall turbulence

<footnote>*Corresponding author</footnote>

<footnote>Email address: miguel.martinez_valero@ensam.eu (Miguel M. Valero)</footnote>
# 1. Introduction

The development of reliable, efficient, and cost-effective numerical tools for accurately predicting multi-physics problems is a timely key challenge in Computational Fluid Dynamics (CFD). In fact, the advancement of predictive tools to investigate flow configurations, including several complex aspects (turbulence, compressibility effects, fluid-structure interaction, transport of active and passive scalars), must take into account requirements exhibiting goal rivalry. Future paradigms in the numerical development of flow solvers strive for increased accuracy of the predictive computational strategies while reducing the computational resources required. This latter criterion is not only associated with a global reduction of the cost and the accessibility of methods on different computational machines and architectures but also with the decrease of the carbon footprint related to calculations in order to comply with progressively more strict regulations. Therefore, efficient numerical modelling in fluid mechanics is essential for progress in fundamental studies as well as industrial applications.

Among the flow problems previously introduced, the accurate numerical representation of near-wall flow features for bodies immersed in turbulent flows needs critical advancement. The prediction of numerous flow features of unstationary flows, such as aerodynamic forces, is driven by the precise representation of localised near-wall dynamics. This aspect is particularly relevant and, at the same time, challenging for the flow prediction around complex geometries. In this case, classical body-fitted approaches may have to deal with high deformation of the mesh elements, possibly leading to poor numerical prediction. Additionally, the simulation of moving bodies may require prohibitively expensive mesh updates. Several numerical strategies have been proposed in the last decades to handle these two problematic aspects. Among these proposals, the Immersed Boundary Method (IBM) [1, 2, 3, 4, 5] has emerged as one of the most popular approaches. The IBM includes a spectrum of tools that operate on non-body-fitted grids, i.e., the mesh is not tailored around the shape of the immersed body, which is used as an internal boundary condition for classical body-fitted approaches. For IBM, numerical strategies are developed instead to mimic the presence of the body. The techniques that rely, for example, on adding source terms in the dynamic equations or imposing discrete values at the centre of grid elements for some physical quantities investigated are usually distinguished by continuous and discrete methods. Despite the significant differences between
the numerous approaches in the literature, state-of-the-art methods can account for body movement and deformation with reduced computational cost. However, the complete resolution of near-wall features, particularly for high Reynolds regimes, is challenging in IBM. Owing to the regularity of the grid and the difficulties in applying arbitrary local directional stretching, a more significant number of mesh elements is usually required in IBM to obtain similar accuracy of body-fitted tools [5].

Another difficulty for the numerical representation of near-wall turbulence, for both IBM and body-fitted methods, is linked to the performance of turbulence/subgrid-scale modelling in the proximity of the wall [6, 7, 8]. Computational resources required for complete flow resolution at the wall may be prohibitive, and, at the same time, wall modelling can show a lack of accuracy in regions exhibiting flow recirculation or strong pressure gradients. In recent years, several approaches from the Estimation Theory [9] have been employed to complete turbulence closures and generalised wall models to obtain universal predictive tools. Among these techniques, data-driven approaches for Uncertainty Quantification and Propagation (UQ-UP) [10], Data Assimilation (DA) [11] and Machine Learning (ML) [12] have been extensively used to improve the accuracy of turbulence closures [13, 14, 15, 16, 17, 18, 19], subgrid-scale models [20, 21, 22, 23, 24, 25, 26] and, very recently, IBM tools [27, 28]. These works have identified theoretical and practical issues and proposed efficient parametric/structural corrections to such models. However, most of these analyses also indicate that a universal predictive model is arguably not achievable because of the high sensitivity to the physical features of the flows as well as to the numerous numerical details of the simulation process. In addition, data-driven methods require the availability and manipulation of data sets to perform their parametric optimisation. The information required from such databases can become prohibitively large when *data hungry* techniques such as deep learning are used. Considering that the optimised turbulence closures via data-driven techniques do not grant universal predictive features, one may wonder if running a classical high-accuracy simulation may be computationally less expensive and more accurate in many cases.

In the present work, a *physics infused* data-driven strategy based on DA is proposed to improve the accuracy of a classical IBM, namely the penalisation model [29]. The optimisation of the IBM is obtained via a parametric inference strategy relying on heterogeneous information, which can include actual measurements at sensors and data derived from physical knowledge of
the flow. The latter is mapped in a data set including local and global features, and it is made available over synthetic sensors, which are appositely created and placed. This information is used to enhance the performance of the model. Whenever possible, exploiting physical information is preferable to storing and using Big Data. First, the loading and manipulation of databases are alleviated, significantly reducing the computational costs. Second, we can quickly achieve a more considerable degree of flexibility about the number and positioning of the sensors. This last aspect can be beneficial in enhancing the convergence rate of the underlying optimisation processes, mitigating critical problems such as over-constraining.

The test case chosen for the present investigation is the turbulent plane channel flow for $Re_τ ≈ 550$. This test case, which has been extensively studied in the literature [6], is driven by the shear effects at the wall. In this framework, the physical knowledge of the near-wall dynamics observed for this test case is infused into the DA algorithm to optimise the IBM. In particular, the Ensemble Kalman Filter (EnKF) [11, 30] is used to infer the parametric behaviour of the penalisation method mimicking the near-wall flow dynamics. As three-dimensional, scale-resolving simulations are performed within the EnKF, the C++ platform CONES recently developed by the team [31] is used to perform the DA strategy on the fly. This means that the integration of the physical information is performed while an ensemble of numerical simulations is running and producing instantaneous flow fields without the need to stop and restart the calculations to perform the optimisation process.

The goal here is not to obtain a universally efficient IBM model, which is arguably a difficult, if not impossible, task. Instead, the main interest of the approach proposed is the capacity to optimise *on the fly* the performance of IBM when streaming data is provided. This scenario is relevant, for example, for digital twins for complex flows [32, 33, 34, 35], where data from a realistic physical application can come in streaming towards a numerical model. In this case, the digital counterpart, which is able to obtain a complete map of the flow, can be used, for example, to determine on-the-fly control laws for the physical twin, where sparse data can be measured from localised sensors. On the other hand, this limited information can be used to augment the accuracy of the prediction of the numerical model in a similar way to the investigation proposed in the present work. While the gap between the present investigation and future real applications of DA-assisted digital twins is large, developing efficient strategies for the assimilation of
sparse data to improve global numerical modelling is beneficial for a wide array of disciplines. One could envision applications in turbomachinery or aerodynamic systems with moving actuators, where the digital twin would enable real-time control based on the sampled flow conditions. This control would help mitigate instabilities and also anticipate or prevent occurrences of extreme events. Other applications include urban flows, where real-time monitoring of pollutant dispersion supports the design of more sustainable urban layouts, or biomedical engineering, where real-time sensor measurements in arteries or respiratory systems can be combined with modelling of physiological flows, contributing to more effective medical decision making. One last application which is strongly connected with IBM is atmospheric reentry, where the debris moves over a trajectory in the physical domain, and continuous surface ablation is observed. Within the perimeter of these applications, an online/in-streaming optimisation of available models such as IBM is essential.

The article is structured as follows. In §2, the numerical equations are presented, together with the formulation of the IBMs used and the DA algorithm. In §3, the test case is introduced and discussed. Preliminary results are presented, and the setup of the DA experiment is detailed. In §4, the main results are discussed. The sensitivity of the data-driven method to changes in its hyperparametric description, as well as the underlying CFD model, is investigated. In §5, the conclusions are drawn, and future perspectives are discussed.

## 2. Numerical Ingredients

### 2.1. *Dynamic equations and numerical solver*

An incompressible flow for a Newtonian fluid can be described by the classical form of the Navier-Stokes equations:

$$ \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) $$

$\mathbf{u}$, $p$, and $\nu$ are the velocity field, the pressure field (normalised over the density $\rho$), and the kinematic viscosity, respectively. The term $\mathbf{f}_P$ represents a normalised volume force, which could be, for example, the IBM source term used to account for the presence of the immersed body.
Numerical simulations are performed using the C++ finite-volume open-source software OpenFOAM [36]. CFD solvers available on this platform, which are based on Finite Volume discretisation, have been extensively used by the scientific community in recent years for both academic studies and industrial applications [37, 38, 39, 40]. The calculations are performed using second-order centred schemes for spatial discretisation, and a second-order backward scheme is selected for time discretisation. The numerical system is resolved using a Pressure-Implicit with Splitting of Operators (PISO) algorithm [41, 42, 43, 44], which is detailed in Appendix B adding the velocity-dependent forcing term in the momentum equations. The simulations performed are run on grids of different resolutions without including subgrid/turbulence closures [6] for the three-dimensional test case. Preliminary tests have shown that, for the meshes used, the inclusion of such additional terms in the equations is responsible for a global degradation of the accuracy.

## 2.2. Immersed Boundary Method

Since the original work by Peskin [1], the term IBM has been used to include a large spectrum of methods that simulate immersed (or embedded) boundaries in viscous flows on grids that do not conform to these boundaries. A qualitative representation is shown in figure 1. The usage of Cartesian grids permits formulating efficient spectral, finite differences, or finite volume approximations of the partial derivative equations used to represent the physical system. In the present work, we focus on an IBM strategy which relies on the usage of source terms included in the dynamic equations to account for the presence of the immersed body.

The approach considered is the classical penalisation method proposed by Angot et al. [29]. This continuous method employs Darcy’s law to obtain a closed expression for the source term $f_P$ in the body interface region $\Sigma_b$ and the solid region $\Omega_b$. To do so, a spatial dependent tensor $D(x)$ must be defined. The resulting penalty term is modelled as follows:

$$f_P = \begin{cases} \mathbf{0} & \text{if } x \in \Omega_f \\ -\nu D(u - u_{ib}) & \text{if } x \in (\Sigma_b \cup \Omega_b) \end{cases} \quad (3)$$

$u_{ib}(x, t)$ is a target velocity representing the immersed body’s physical behaviour. If the body surface is not moving, then $u_{ib} = \mathbf{0}$. The components $D_{ij}$ of the tensor $D$ are usually controlled to obtain the best compromise
<figure><img src="image_3.png" /><figcaption>Figure 1: Qualitative representation of a Cartesian grid used for IBM calculations. $\Omega_f$ and $\Omega_b$ represent the fluid and body domains, respectively. The blue area $\Sigma_b$ represents the body interface.</figcaption></figure>

between accuracy and stability of the numerical solver [5]. The choice performed is crucial in the proximity of the grid transition between the flow $\Omega_f$ and the interface $\Sigma_b$ regions.

This study considers a second IBM approach to quantify the limitations of the previously described penalty method, aiming to assess its efficacy and pinpoint potential constraints or drawbacks within the simulation framework. This second method was already validated in *OpenFOAM* by Constant [40] and is based on the work by Uhlmann [45] and Pinelli *et al.* [46], but it is improved by including a *Reproducing Kernel Particle Method* (RKPM) [47]. A complete description of this method is reported in Appendix C. This discrete IBM relies on two complementary physical spaces. The *Eulerian* domain is described by the grid used for the calculation, while the *Lagrangian* Markers are discrete points representing the immersed body's surface. The method determines the source term via a two-step procedure, which is usually referred to as *interpolation* and *spreading*. The main difference between the two algorithms is that the penalisation method is strictly local, i.e., the forcing depends exclusively on the flow field predicted in the correspondent mesh element. On the other hand, the discrete method relies on an interpolation stencil to communicate between the *Eulerian* grid and the *Lagrangian* Markers. The user can select the size of the grid elements for this structure. Larger interpolation stencils connect a more significant number of mesh elements with each *Lagrangian* Marker, which improves the numerical algorithm's stability and the solution's smoothness. However, larger compu-
tational stencils are also responsible for higher computational requirements. In addition, the size of the computational stencil (for the discrete method) and the width of the interface region $\Sigma_b$/value of the coefficients $D_{ij}$ (for the penalisation model) can be responsible for diffused interfaces, which can affect the precision of the numerical results.

## 2.3. Sequential Data Assimilation: The Ensemble Kalman Filter

Data Assimilation (DA) techniques [11, 48] have emerged in recent times as a powerful tool to enhance the reliability of numerical simulations in Fluid Dynamics. These approaches enable the integration of high-accuracy observation (DNS or experimental data, for example) into reduced-order numerical models to obtain improved predictions of the underlying physical system. DA approaches are usually grouped into two leading families: variational and sequential tools. The former, which includes well-known methods such as 3DVar and 4DVar, resolves the DA problem as an optimisation task in which initial conditions and dynamic models are determined to minimise a prescribed cost function. These models are extremely accurate but may not converge to analyse multi-scale time-evolving physical processes [49]. Applications in fluid mechanics mainly deal with stationary configurations [50, 51, 52]. On the other hand, sequential methods mostly rely on Bayesian techniques to resolve the DA problem. One of the most powerful methods in sequential DA is the Ensemble Kalman Filter (EnKF) [30, 53], which is an ensemble-based tool derived from the Kalman Filter [54]. This technique has proved to be efficient in reconstructing turbulent flows in recent years for both stationary and unstationary configurations [25, 26, 55, 56]. Additionally, it has previously been utilised with the immersed boundary projection method for two-dimensional laminar flow scenarios [57, 58].

The key elements representing the EnKF are now introduced. A complete description of this method is reported in Appendix A. Discussion here is mainly devoted to general DA procedures, while how the algorithms are integrated in the context of the present study is presented in §3.2. Let us consider the evolution of the physical variables of the flow from the instant $k$ to $k+1$. A forecast ($f$) step is performed, where an ensemble of state vectors $\mathbf{u}_i$ ($i \in [1, m]$), where $m$ is the number of ensemble members) is advanced in time using a model $\mathcal{M}$:

$$ \mathbf{u}_{i,k}^f = \mathcal{M}_{i,k:k-1}(\theta) \mathbf{u}_{i,k-1} \quad (4) $$
Here, $\theta$ represents a set of free parameters that can be used to calibrate the model. In the present work, $\mathcal{M}$ is, for example, the discretised set of the Navier-Stokes equations in the CFD solver, while $\theta$ is some parametric description affecting the model prediction such as the coefficients of the tensor $\mathbf{D}$ in equation 3. If observation (i.e., high-fidelity data) is available at time $k$, an analysis (a) phase is performed, where the DA procedure updates the system state. This observation usually consists of local values for physical quantities sampled in experiments or high-fidelity simulations. The latter is the case for the present work. The Kalman gain $\mathbf{K}_k$ governs the state update in the analysis phase. This matrix takes into account the flow variables predicted in the forecast, as well as the confidence in the accuracy of the model and of the observation:

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

The matrix $\mathbf{Y}_k$ includes the observations available, while $\mathcal{H}$ represents a mathematical operator used to project the model state in the observation space. The following state-of-the-art modifications are also included in the present EnKF formulation:

* The EnKF is naturally used to update the system's physical state, but it can also be used to update free coefficients $\theta$ governing the model $\mathcal{M}$. In the present work, this task is performed via an extended state approach [11]. In practice, the parameters are included in the state vector $\mathbf{u}'_i = [\mathbf{u}_i \theta_i]^T$ for each realisation, and they are updated with the physical field using equation 5.

* Some deterministic inflation [11] is included to increase the variability of the state predicted by the EnKF:

$$ \mathbf{u}_{i,k}^{a} = \bar{\mathbf{u}}_{k}^{a} + \lambda_{k} \left( \mathbf{u}_{i,k}^{a} - \bar{\mathbf{u}}_{k}^{a} \right) \quad (6) $$

where $\lambda_k \ge 1$. This procedure usually improves the rate of convergence and the accuracy of the updated system's state $\mathbf{u}_k^a$. However, the intensity of inflation has to be carefully calibrated. Very high values may produce unwanted spurious oscillations, which can cause divergence of the DA procedure. However, small values may preclude the proper identification of local/global optimised configuration, particularly if the prior estimate is not accurate.
* Covariance localisation of the Kalman gain $K_k$ is used to gradually set to zero the correlation between the observations and the system's state with increasing distance. This choice controls the emergence of some undesired spurious fluctuations due to sampling errors. Covariance localisation is usually performed premultiplying $K_k$ with a used-defined matrix $L_k$. This matrix, for which $[L]_{ij} \in [0, 1]$, is defined using an exponential decay function in space so that the update of the model state tends to zero when variables calculated far from the sensors are considered. The elements $i, j$ of the matrix are obtained as:

$$[L]_{i,j} = e^{-\left(\frac{\Delta x_{ij}^2}{\eta_x} + \frac{\Delta y_{ij}^2}{\eta_y} + \frac{\Delta z_{ij}^2}{\eta_z}\right)} \quad (7)$$

where the terms $\Delta x_{ij}$, $\Delta y_{ij}$, $\Delta z_{ij} \ge 0$ represent the distance in the $x$, $y$ and $z$ direction between the point of evaluation of the model and the location of the sensors. The hyperparameters $\eta_x, \eta_y, \eta_z$ are used to control the behaviour of the matrix $L_k$. If the computational grid used by the model and the location of the sensors do not change in time, $L$ is also time-independent. This is the case for the present investigations.

Looking at equations 3 and 5 for the IBM and for the EnKF, respectively, one can see a similarity in the way the two terms may affect the prediction of the velocity field. However, their effect is complementary. The penalisation method provides a local and linear effect with respect to the velocity in the mesh element. The intensity is controlled by the coefficients of the stiffness tensor. On the other hand, the state update in equation 5 is i) strongly non-local, as information at every sensor contributes to the state update in the whole localised domain and ii) non-linear with the velocity, as the Kalman gain is obtained via numerous matrix operations involving products of the velocity field. Therefore, the state estimation and the parametric optimisation performed by the DA tools look for a balance between the two elements, providing the best accuracy and reducing the discrepancy of the numerical model with the observation.

### 3. Test case: turbulent plane channel flow, $Re_\tau = 550$

The turbulent plane channel flow is a fundamental benchmark test case in CFD, and open databases are available online, providing high-accuracy DNS data for reference. On the one hand, the simple geometry of this test
case allows the exclusion of a number of physically complex aspects in the
optimisation process, such as, for example, the separation of the boundary
layer. On the other hand, because of the velocity fluctuations observed in
the wall's proximity, traditional IBMs usually fail to provide an accurate flow
prediction. For all of these reasons, this test case is a suitable candidate for
an ambiguous assessment for IBM augmentation via DA.

The friction Reynolds number considered is $Re_τ = u_τh/ν = 550$. The friction velocity is defined as $u_τ = \sqrt{τ_w/ρ}$ where $τ_w$ is the wall shear stress, while $h$ is the half-height of the channel. The coordinate system is set so that $x$ is the streamwise direction, $y$ is the wall-normal direction, and $z$ is the spanwise direction. The two walls are positioned at $y = 0, 2h$, respectively. Periodic boundary conditions are applied to the side walls, and the mass flow rate is conserved in time using a source term for the dynamic equations already implemented in *OpenFOAM*. Results are usually presented in non-dimensional form (suffix +) using the friction velocity $u_τ$ and the viscous wall unit $δ_ν = ν/u_τ$ for normalisation. A second non-dimensional form used in the present analysis (suffix ⋆) relies on $u_τ$ calculated by a reference simulation (R-DNS-BF) introduced in §3.1. In addition, statistical moments are also averaged over the two directions $x$ and $z$, for which statistical homogeneity is observed.

3.1. *Reference simulation and preliminary results*

A database of classical numerical simulations (i.e., without data assimila-
tion) has been performed to obtain preliminary results. Details are reported
in the first four rows of table 1. The physical domain is discretised using a
uniform distribution of the elements in the *x* and *z* directions. A geometric
distribution is used in the normal direction, increasing the grid size while
moving away from the wall. Therefore, $\Delta y_{min}^*$ represents the grid size at
the wall, and $\Delta y_{max}^*$ is the size at the centre of the channel. This analysis
has been performed to obtain a suitable initial *prior* state to be optimised
with DA. The database includes several simulations obtained with different
techniques:

* One body-fitted DNS (R-DNS-BF), which represents the reference sim-
ulation that will be used to validate the results. Statistical moments for
this simulation have been obtained for 300 advective times t<sub>A</sub> = h/U<sub>c</sub>,
where U<sub>c</sub> is the averaged velocity at the centreline for y = h. The
averaged velocity profile ⟨U⟩<sup>+</sup> = ⟨u<sub>x</sub>⟩/u<sub>τ</sub> and the components of the
Reynolds stress tensor $\langle u_i' u_j' \rangle^+ = \langle u_i' u_j' \rangle / u_\tau^2$ compare well with results by Del Álamo and Jiménez [59], and Hoyas and Jiménez [60] for similar $Re_\tau$ as shown in figure 2. The minor discrepancies observed could be due to a small difference in $u_\tau$ (around 2.55%), to the rate of convergence of the statistical moments, to the different discretisation strategies or to the different distribution of the mesh elements in the wall-normal direction.

* One body-fitted DNS (DNS-BF), which is run on a smaller domain and uses a lower grid resolution when compared with the simulation R-DNS-BF. Details of the grid are provided in table 1. One can see that a factor of 4 coarsens the resolution in the x and z directions. In the vertical direction, the number of mesh elements is also reduced by a factor of 4, and the ratio coefficient of the geometric distribution of the elements is different. This last choice has been performed to obtain a resolution similar to the wall region compared to the refined grid used for the simulation R-DNS-BF.

* An IBM simulation with the continuous forcing (penalisation) method, which is performed on a grid very similar to the one used for DNS-BF shown in figure 3. The differences emerge in the proximity of the wall. For the IBM, the wall mesh element has its centre in $y = 0$ (or $y = 2h$ at the top), and three additional layers of mesh elements of the same size are included in the opposite direction of the flow domain. These layers are placed to ensure the algorithm's numerical stability. We observed that including more solid layers would not change the results of the calculations while increasing the computational costs. The calculation is referred to as DNS-IBM-CF. The bottom and top walls are not moving, so $u_{ib} = 0$. The volume force $f_P$ is non-zero at the interface region $\Sigma_b$ (light blue region in figure 3(b)), which consists of the three closest mesh elements in the y direction to the immersed walls ($y = 0, 2h$), and the solid domain $\Omega_b$, which is in this case represented by a layer of two cells in the y direction (grey in figure 3(b)). The tensor $D$ determining the behaviour of the IBM is set so that $D = D_I$, where $D = 10^5$ for the mesh elements in $\Sigma_b \cup \Omega_b$ and $D = 0$ elsewhere. These values are chosen according to studies reported in the literature [5, 29].

* An additional IBM run performed with the discrete forcing method on
<table><caption>Table 1: Grid resolution, domain size and time step $\Delta t$ used for the generation of the database of simulations. The calculated friction velocity $u_{\tau}$ is also reported for every simulation.</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<sub>min</sub><sup>*</sup></th><th>Δy<sub>max</sub><sup>*</sup></th><th>L<sub>x</sub>/h</th><th>L<sub>z</sub>/h</th><th>u<sub>τ</sub></th><th>Δt (t<sub>A</sub>)</th></tr></thead><tbody><tr><td>R-DNS-BF</td><td>1024 × 256 × 512</td><td>9.9</td><td>6.6</td><td>1</td><td>11.2</td><td>6π</td><td>2π</td><td>0.0480</td><td>0.004</td></tr><tr><td>DNS-BF</td><td>128 × 58 × 64</td><td>39.5</td><td>26.3</td><td>1.09</td><td>56.6</td><td>3π</td><td>π</td><td>0.0430</td><td>0.02</td></tr><tr><td>DNS-IBM-CF</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.0375</td><td>0.02</td></tr><tr><td>DNS-IBM-DF</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.0397</td><td>0.02</td></tr><tr><td>DNS-IBM-DA</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.0487</td><td>0.02</td></tr></tbody></table>

the same grid as DNS-IBM-CF and referred to as DNS-IBM-DF. The interface region $\Sigma_b$ is chosen to be the same for the two methods in order to provide a rigorous comparison. The Lagrangian Markers for the discrete method are positioned at the centre of each mesh element for $y = 0, 2h$, and the computational stencil is made by $3 \times 3 \times 3$ elements.

Results obtained from the simulations of the database are analysed to evaluate the accuracy of the calculations performed with the coarse grids and to assess the efficacy of IBM in this case. Comparisons are performed against the simulation R-DNS-BF, which is considered to be the true physical state. First, the predicted friction velocity $u_{\tau}$ is studied. The values obtained by the different simulations, which are reported in table 1, show that the coarse-grained simulations significantly under-predict $u_{\tau}$. The discrepancy in predicting this quantity is substantially more prominent for the two IBM simulations when compared with the run DNS-BF. This observation confirms IBM's difficulties in accurately representing the physical features of the flow for this test case.

Similar considerations can be drawn via the analysis of the mean streamwise velocity as well as the components of the Reynolds stress tensor, which are shown in figure 4. The investigation of the velocity field, which is depicted in figure 4(a) and (b), reveals that the prediction obtained using the two IBM strategies is similar. However, one can see in figure 4(b) that, using the values prescribed for $D$, the IBM based on penalisation is not able to respect the no-slip condition $u_{ib} = 0$ for the velocity at the wall. The simulation DNS-BF obtained a marginally more accurate prediction of the velocity
<figure><img src="image_4.png" /><figcaption>Figure 2: Comparison of the main statistical moments of the velocity field. Results are shown for the simulation (—) R-DNS-BF and (---) the database by Del Álamo and Jiménez [59].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: (a) General view of the grid employed for the plane channel with IBM. (b) Cartesian grid employed when using the IBM to analyse the turbulent plane channel. $\Omega_f$ and $\Omega_b$ represent the fluid and body domains, respectively (the latter displayed in grey). The body interface is the blue area $\Sigma_b$ (three cell layers in the y direction).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Comparison of the main statistical moments of the velocity field. Results are shown for simulations (—) DNS-BF, (···) DNS-IBM-CF, (---) DNS-IBM-DF, and (—) R-DNS-BF.</figcaption></figure>

field, mainly because of the more precise prediction of $u_{\tau}$. However, the discrepancy with the reference R-DNS-BF simulation is significant for all the other database runs. The analysis of the components of the Reynolds stress tensor, which are shown in figure 4(c)–(f), also indicates a low accuracy for the coarse-grained simulations. One could deduce that the simulation DNS-IBM-CF is globally the least accurate, even if significant differences can be observed for every component of the Reynolds stress tensor.

In summary, the present analysis shows a global lack of accuracy in the coarse-grained simulations. In particular, the IBM run DNS-IBM-CF, which uses a penalisation term, provides low accuracy and is not able to obtain the essential no-slip condition $u_{ib} = 0$ at the wall.

## 3.2. Data Assimilation experiment

A DA strategy is proposed here to infuse physical knowledge within an IBM numerical solver. The online analysis aims i) to improve the accuracy of such a solver and ii) to do so with a limited increase in computational resources required. Discussion in §2.3 indicated how the DA methods rely on a *model*, providing a quasi-continuous representation of the physical variables
in the domain of investigation, as well as some local observation. These
elements are now discussed in detail.

The *model* here considered is the numerical solver where the test case used for the preliminary simulation is DNS-IBM-CF. As seen in §3.1, such a model was not able to provide accurate estimations of the statistical moments. In addition, a no-slip condition at the wall was not obtained. The number of ensemble realisations, which permits exploring the spaces associated with the parametric description and the physical solution, is set to *m* = 40. The value chosen for *m* is based on recommendations provided in the literature from analyses combining the EnKF with CFD solvers [25, 61]. The data-driven strategy is used to dynamically enhance the velocity field of the whole domain **u** at each analysis phase as well as to infer the local value of the diagonal components of the tensor **D**. Therefore, the entire parametric space to be optimised consists of the three components *D*<sub>xx</sub>, *D*<sub>yy</sub> and *D*<sub>zz</sub> for 128 × 64 × 64 = 524 288 grid elements, for a total of 1 572 864 degrees of freedom to be determined in the parametric space investigated. However, for IBMs, coefficients for mesh elements ∈ Ω<sub>f</sub> are always set to zero. This choice is conserved here, and therefore, the optimisation procedure is performed for the mesh elements in the solid Ω<sub>b</sub> and the interface Σ<sub>b</sub> regions. This implies that, out of the 64 layers in the *y* direction, only 5 at the bottom and 5 at the top are considered. In addition, the following strategies are applied to reduce the complexity of the problem:

1. Statistical symmetry around the plane $y = h$ is used to consider only 5 layers. While symmetry is not observed for the instantaneous flow, viscous phenomena in the wall's proximity reduce the intensity of velocity fluctuations, mitigating the effect of this approximation.

2. Similarly, statistical invariance due to homogeneity in the *x* and *z* is used to neglect the dependency of the coefficients to those directions, i.e., $D_{ii}(\mathbf{x}, t) = D_{ii}(y, t)$ for $i = x, y, z$.

Thus, the optimisation task is reduced to a space of $\theta = 15$ degrees of freedom, which is the determination of 3 constants for 5 grid layers in the $y$ direction. A second essential aspect is the *prior*, which is used as the initial condition for the model. The prior, which includes the physical condition imposed at the initial time as well as the free parameters prescribed, plays an essential role in the rate of convergence towards the optimised solution by DA [57]. For this work, the initial physical state imposed for each ensemble member is obtained from interpolating the simulation R-DNS-BF for
$t=0$. In addition, the velocity field for the mesh elements in the subviscous layer has been set to $\mathbf{u} = [u_{\tau}y/\delta_{\nu}, 0, 0]$. The initial state for the parametric description $\theta$ is now discussed. Due to the absence of data from the reference simulation, the values $D_{ii}$ utilized in the ensemble runs are obtained using a Gaussian distribution $D_{ii} = \mathcal{N}(\bar{D}_0, \sigma_0^2) = \mathcal{N}(10^5, 2.5 \cdot 10^7)$ for the mesh elements $\in \Omega_b \cup \Sigma_b$. This implies that each ensemble member is obtained imposing a perturbation $\sigma_0 = 0.05 \bar{D}_0$ to the coefficients used for the preliminary simulation DNS-IBM-CF, for which $D_{ii} = \bar{D}_0 = 10^5$ for every mesh element in the solid and in the interface regions. This variance observed over the ensemble realisations enables the exploration of optimised configurations that may exhibit important differences with the prior conditions while ensuring they remain within the bounds of physically plausible solutions. These conservative choices for physical state and parametric description have been performed to unambiguously identify the sensitivity of the solution to the parametric variation, speeding up the optimization procedure by EnKF. Preliminary tests have shown how the choice for the prior physical state does not affect the convergence of the algorithm towards the optimum if a physically acceptable solution is imposed. On the other hand, as for every DA tool [11], the initial choice for the problem's parametric description can affect the optimum rate of convergence and its position in the parametric space investigated. First tests spanning several initial conditions for the parameters ($D_{ii} \in [10^4, 5 \cdot 10^5]$) have shown maximum discrepancies of the order of 5% in the physical states predicted with the optimised coefficients, indicating a satisfying level of robustness of the algorithm for this test case. However, additional initial conditions, such as $D_{ii} \ge 10^6$, may lead to stabilisation issues, while setting $D_{ii} = 0$ results in degradation of the flow physics.

The observation used in the present analysis is now described. The information employed deals with the flow's physical behaviour in the wall's proximity. First, synthetic sensors are created at the wall to obtain the no-slip condition for the velocity field, which the classical penalisation IBM tool did not correctly predict. A qualitative representation of the distribution of the synthetic sensors for $y=0$ is shown in figure 5(a). A total of 4096 sensors are used for which the condition $u_x = 0$ is imposed as surrogate observation. No constraint is imposed for $u_y$ and $u_z$. This choice complies with the intrinsic limitations of the EnKF, for which a matrix inversion of the size of the observed data [$n_o, n_o$] must be performed [11]. Also, thanks to the solenoidal features of the resolved numerical schemes, the inferred field for $u_x$
also affects the other velocity components without the need to over-constrain the system to be optimized. We could have also included synthetic sensors in the solid region for $y < 0$, but this would have incremented the computational resources required for the EnKF without significantly improving the accuracy. A directly measured physical information, which is the instantaneous friction coefficient $C_f$, is also used as observation. Its relation with the shear stress at the wall $\tau_w$ and $u_\tau$ is the following:

$$ \tau_w = \rho\nu \left. \frac{\partial\langle u_x \rangle}{\partial y} \right|_{y=0} = \rho u_\tau^2 = \frac{1}{2} \rho C_f U_c^2 \quad (8) $$

where $U_c$ is the mean velocity at the centreline for $y = h$. The friction coefficient $C_f(t)$ sampled during the simulation R-DNS-BF is thus used to augment the IBM via the application of the EnKF. The usage of this last quantity is actually a direct emulation of what would happen in a digital twin for a realistic application. In fact, one can envision an experiment in a wind tunnel in which forces affecting an immersed body (car model, building model) are measured by a six-component balance in real time. These measurements can be manipulated to obtain the correspondent drag and lift coefficients, which can be forwarded in-streaming to the EnKF. Therefore, $C_f(t)$ obtained from simulation R-DNS-BF can be thought of as a replica of a typical output from an experiment. The confidence in the surrogate and measured observation is now discussed. First of all, the uncertainty affecting each sensor is supposed to be statistically independent so that the matrix $R$ in equation A.5 can be reasonably approximated to a diagonal matrix. Gaussian noise is artificially used to perturb the observations. For the synthetic sensors associated with the no-slip condition at the wall, the standard deviation of such noise has been set to $0.05 U_c$. For the sensor measuring the friction coefficient, the standard deviation is set to a uniform value of $0.05 \langle C_f \rangle$. This value is very close to the variance of the time evolution of $C_f(t)$ observed for the simulation R-DNS-BF (see figure 5(b)).

The DA procedure combining the *model* and the *observation* is now detailed. The sequential features of the EnKF are fully exploited, which means that the DA state estimation and optimisation are performed *on the fly* during the run. For simulations of this size, non-intrusive approaches are too expensive in terms of the computational resources required; in fact, simulations of the ensemble must be interrupted before every analysis phase and restarted afterwards. This operation time, which sums up the writing and reading of files for the physical solution, may increase the computational
<figure><img src="image_7.png" /><figcaption>Figure 5: (a) Position at the bottom wall of the sensors providing constraints for the streamwise direction ($u_x = 0$) when applying the EnKF. (b) Temporal evolution of $C_f(t)$ of R-DNS-BF, used as observation for the Data Assimilation experiment.</figcaption></figure>

costs by several orders of magnitude when compared with the forecast of the solution, and it is not compatible with dynamic flow estimation and control such as the one performed in digital twins. For this reason, the numerical simulations of the ensemble are coupled online with the DA code using the application *Coupling OpenFOAM with Numerical EnvironmentS* (CONES), recently developed by the research group [31] and uploaded on https://github.com/MiguelMValero/CONES. CONES parallelises the problem using an open-source coupler called CWIPI [62], which identifies two modules (simulations in OpenFOAM and the library with the EnKF algorithm) and establishes HPC communications among them according to the user's specifications. CPU cores assigned for the module of OpenFOAM are occupied when carrying out the forecast step and become inactive during the analysis phases, whereas those destined for the correction step work oppositely. Thus, this algorithm completely excludes the costly operations of interrupting/restarting the simulations of the ensemble. In this work, a DA analysis phase is performed every six time steps of the numerical forecast. Considering that the time step $\Delta t = 0.02t_A$, this implies that data are assimilated every $0.12t_A$. This relatively high frequency of assimilation with respect to the characteristic physical scales of the flow naturally excludes the risk of lack of convergence of the DA procedure [22]. As detailed in Appendix B, the pressure is updated from the system's state estimated by the EnKF by means of a Poisson equation. This additional step grants conservativity of the discretized Navier-Stokes equations. A visualization of the EnKF procedure is shown in figure 6.
<figure><img src="image_8.png" /><figcaption>Figure 6: Scheme of CONES application for sequential DA techniques.</figcaption></figure>
Finally, covariance inflation and localisation are discussed. These state-of-the-art procedures aim to improve the accuracy and robustness of the EnKF, as well as reduce the computational costs associated [11]. For the former, deterministic inflation is applied by using equation 6 with a constant $\lambda = 1.05$ during initial evolution phases for $t \in [0, 5 t_A]$. This choice increases the variability of the ensemble and prevents the collapse of the parametric description of the system during the first analysis phases. Localisation is applied by taking into account that observation/physical constraints are located in proximity to the wall. A matrix $\mathbf{L}$ premultiplying the Kalman Gain in equation A.5 is generated where all the coefficients are zero with the exception of those referring to the mesh elements located in the subviscous layer ($y^+ < 5$) and outside the physical space. These coefficients are determined using the decay exponential function described in equation 7. The parameters $\eta_x, \eta_y$, and $\eta_z$ are tuned to mitigate discontinuities during the update of the DA state, especially while transitioning from the subviscous layer to the outer wall regions. This discontinuity, limited to 0.3%, necessitates a $\eta_y$ value of $1.37 \cdot 10^{-5}$. Similar strategies are implemented in the remaining two directions, aligning with the periodic boundary conditions (specifically, $L_x/h = 3\pi$ and $L_z/h = \pi$), which results in values of $\eta_x = 15.29$ and $\eta_z = 1.70$, respectively. Physical localisation of the domain is also jointly performed to exclude the regions where the updates are zero for the EnKF calculation. This clearly reduces the computational costs demanded. The region conserved corresponds to the first seven cell layers in the wall-normal direction, that is, a total number of 57 344 mesh elements against the 524 288 elements of the complete mesh. The state vector is composed of $n = 3 \cdot 57 344 = 172 032$ degrees of freedom for the velocity field plus the fifteen coefficients describing the tensor $\mathbf{D}$. The choice of the parameters for localisation excludes from the outer-wall cells the EnKF state update, which is, however, affected by the near wall state correction during the forecast of the model, thanks to the non-local, non-linear features of the Navier-Stokes equations [63].

## 4. Results

The physics-infused procedure developed relies on on-the-fly state estimation and parametric optimisation using an ensemble of numerical simulations. These runs predict the instantaneous turbulent physical field. Results obtained by this procedure, whose associated run is referred to as DNS-IBM-DA, are first compared with the reference database presented in §3.1. A
sensitivity analysis of the hyperparameters governing the numerical model's performance and the DA strategy is also investigated.

### 4.1. *Comparison of DA results with reference simulations*

It is here reminded that the DA procedure based on the EnKF is used to optimize the parametric behaviour of an IBM penalization model, as well as to perform an update of the flow field. These operations are performed sequentially, operating directly over the predicted instantaneous flow field.

First of all, data in table 1 indicates that the infusion of the local no-slip constraint and the global knowledge of the friction coefficient $C_f$ is able to significantly improve the prediction of the friction velocity $u_τ$. In fact, the DNS-IBM-CF run obtains a significant underprediction of this latter quantity with a discrepancy of 21.88%. The physics-infused procedure obtains an overprediction of 1.46%, reducing the discrepancy by almost 15 times when compared with the classical penalisation IBM. Therefore, one can expect that the global accuracy of the flow is significantly increased, considering the essential role of $u_τ$ in the establishment of the physical features observed for this test case.

A first qualitative comparison of the flow features is shown in figure 7 for the isocontours calculated via the Q-criterion. One can see that, despite the lack of grid resolution, the DA simulation (DNS-IBM-DA) is able to capture fine structures that exhibit a better agreement with the reference DNS (R-DNS-BF). This result implies that the instantaneous flow update performed by the EnKF provides a significant improvement in accuracy when compared with the classical IBM approaches, which are run on grids of identical resolution. Similar conclusions can be drawn comparing the statistical moments obtained by the simulations, which are shown in figure 8. The normalized mean velocity profile $\langle U \rangle^+$ shown in figure 8(a) and (b) is significantly closer to the DNS reference both in the inner and the outer layers. In addition, one can see that the no-slip condition at the wall is well obtained by the DA simulation, within the confidence level prescribed for the observation. Improvements using the DA method can also be observed for the Reynolds stress tensor components, in particular when compared with the continuous IBM model. These improvements, however, vary depending on the component considered. High accuracy is obtained for the prediction of $\langle u'_x u'_y \rangle^+$, which is shown in figure 8(f). This term is strictly connected with the friction coefficient $C_f$/friction velocity $u_τ$, which is an observed variable. Therefore, the DA optimization acts more directly on predicting $\langle u'_x u'_y \rangle^+$. In addition, the
<figure><img src="image_9.png" /><figcaption>Figure 7: Isocontours of Q-criterion calculated for $y/h = 0.18$ ($y^* \approx 100$). Results are shown for simulations (a) DNS-IBM-CF, (b) DNS-BF, (c) DNS-IBM-DA, and (d) R-DNS-BF.</figcaption></figure>
estimation of this quantity via the penalisation IBM is already adequately accurate. Hence, the prior state is not far from the optimised state. On the other hand, a small degradation is observed for $\langle u'_z u'_z \rangle^+$ in figure 8(e).

The analysis is completed by comparing the time spectra shown in figure 9. The spectra are calculated using the 1-D discrete Fast Fourier Transform (FFT) [64] to the time-series of the fluctuating velocity $\mathbf{u}'(t)$ sampled at the locations $y^* = 30, 56$:

$$FFT(k) = \sum_{k=0}^{N-1} e^{-2\pi j \frac{kt}{N}} \mathbf{u}'(t) \qquad (9)$$

$N = 7500$ is the total number of samples. The sampling time is $\Delta t = 0.04 t_A$ over a period $T = 100 t_A$. Using Taylor's hypothesis for frozen turbulence [65], results in figure 9 are shown for wavenumbers $\kappa = 2\pi f/U_c$, where $f$ is the temporal frequency of the FFT transform. In addition, $\kappa$ is adimensionalized with respect to the kinematic viscosity $\nu$ and the friction velocity $u_\tau$ so that $\kappa^+ = \kappa\nu/u_\tau$. To smooth out the curves, a first-order Butterworth low-pass filter [66] is used to eliminate the undesired noise. One can see that the coarse-grid simulations IBM-DNS-CF and DNS-BF produce similar spectra, and the fluctuation energy starts to decay for relatively low wavenumbers. This observation is related to two concurring phenomena: the lack of grid resolution and the poor accuracy in estimating the wall shear stress. On the other hand, DA results get significantly closer to the reference DNS results, indicating that the accurate prediction of $u_\tau$ provides a beneficial effect on the representation of the velocity fluctuating field. One could also observe that the gap between the reference DNS and the DA run is around half-octave, i.e., proportional to the difference in resolution of the grids used for calculations. This confirms that the usage of the EnKF was able to compensate for the modelling error efficiently and that the discrepancies observed are intrinsically tied to the grid resolution. The accurate performance of the EnKF for this case is tied with the usage of multiple sources of information (no-slip condition and $C_f$) and its capabilities to handle them efficiently, accounting for their confidence level. In fact, tests using only the synthetic velocity sensors did not provide the same level of accuracy as the simulation presented, and they showed a higher sensitivity to DA hyperparameters and the prior state. This favourable feature is not necessarily granted by other DA approaches or Machine Learning tools, for which handling different types of data may lead to issues for the convergence of the trained models
<figure><img src="image_10.png" /><figcaption>Figure 8: Comparison of the main statistical moments of the velocity field. Results are shown for simulations (—) DNS-IBM-DA, (···) DNS-IBM-CF, (—) DNS-BF, and (f) R-DNS-BF.</figcaption></figure>

[67]. Also, it is worth discussing that for some IBM proposals in the literature, the no-slip condition and the wall shear stress act as dual variables [68]. This deterministic relation appears to be approximated in a stochastic sense by the EnKF approach used here. The precision of such approximation is connected with the level of confidence prescribed for the model and for the observation.

## 4.2. Sensitivity of the physics-infused strategy to mesh refinement and placement of sensors

The analysis of the results performed in §4.1 highlighted how the infusion of physical information in the numerical process can improve the quality of the results even when simulations are performed using coarse grids. In the present section, the sensitivity of the DA strategy is tested against variations in the key elements that constitute the global methodology. More precisely, two aspects are considered. Variations in the grid resolution for the model represent the first one. The second one consists of a different distribution and density of the sensors providing observation. Four additional DA runs are performed for this purpose, varying the indicated parameters. Features of such simulations are reported in table 2. The simulation DNS-IBM-DA
<figure><img src="image_11.png" /><figcaption>Figure 9: FFT calculated by sampling the fluctuating velocity field $\mathbf{u}'$ located at (first row) $y^* \approx 30$ and (second row) $y^* \approx 56$. Results are shown for the simulations (—) IBM-DNS-DA, (···) IBM-DNS-CF, (—) DNS-BF, and (—) R-DNS-BF.</figcaption></figure>

from the previous section is now taken as the reference for comparison with
the rest of the data-driven simulations. As one can see, case 1 uses the same
mesh as the reference DA run, while cases 2, 3, and 4 are performed using a
modified grid. For these cases, the resolution in the streamwise direction $x$
and in the spanwise direction $z$ is two times coarser, while a higher resolution
is employed in the normal direction $y$. For the latter, a smaller expansion
ratio between consecutive mesh elements $r_y = \Delta y_i/\Delta y_{i-1}$ is used. This
parameter, which is equal to $r_y = 1.171$ for the grid employed in the DNS-
IBM-DA run, is reduced to $r_y = 1.060$ in these studies. This new grid is
composed of $N = 262144$ elements. The second aspect that is considered
is the number of synthetic sensors employed to locally inform the model to
respect the constraint $u_x = 0$ for $y = 0$. One can see in figure 10 that
different densities of the sensors, as well as different distributions, have been
analysed. In particular, one can see that probes have been positioned both
at the centre of the mesh element and at their interface. In addition, in case
3, the problematic aspect of multiple sensors within one mesh element has
been investigated.

The prior state for the new DA runs has been chosen using the same
criteria previously presented for the DNS-IBM-DA case, and the optimisation
targets the values of the diagonal elements of the tensor $D$. The prior state for the parameters is, for these cases, the solution obtained by the DNS-IBM-DA run. This choice, which has been performed to obtain a faster and more robust convergence of the optimisation procedure, also allowed for a lower deterministic inflation level to be applied. In this case, a constant value of $\lambda = 1.01$ is chosen for all parameters throughout the time range $t \in [0, 300 t_A]$. The level has been selected empirically to obtain the best trade-off between fast convergence and stability of the optimisation procedure.

Figure 11 illustrates the main statistical moments of the flow. One can see that the choice of a different grid affects the prediction of the physical quantities investigated, while a weak sensitivity to the position and number of sensors is observed. In particular, results for case 1 are very similar to those obtained by the run DNS-IBM-DA despite the smaller number of sensors located at the cell edges. A quantification of the differences between these two runs is provided by the error in the prediction of the friction coefficient, calculated as:

$$ \Delta_{C_f} = \frac{|\langle C_{f,\text{DA}} \rangle - \langle C_{f,\text{R-DNS-BF}} \rangle|}{\langle C_{f,\text{R-DNS-BF}} \rangle} \quad (10) $$

The results, which are reported in table 2, indicate a discrepancy of around 3% for the reference DA run and 1% for case 1. The analysis is now extended to the three DA runs performed over the modified grid. Results obtained for the statistical moment of the velocity field in figure 11 indicate that degradation of the performance of the DA algorithm is observed. This observation is expected because the IBM is, in this case, less accurate owing to the lower grid resolution. However, it has been verified that the DA runs obtained using this grid are sensibly better than classical simulations on the same mesh, confirming the improvement in accuracy of the DA strategy, whatever numerical model is used for the forecast. One can see that results from cases 2 and 4 are almost identical, also in terms of measured discrepancy $\Delta_{C_f}$. On the other hand, results for case 3 are less precise, with an increase of $\Delta_{C_f}$ when compared to the other two DA realisations on the same grid. A direct conclusion that can be drawn from this analysis is that while features of the numerical model used for DA are important, the number and position of sensors can also play a role in the global accuracy of the results. More precisely, when sensors are positioned too close to one another, the uncertainties affecting the measurements can be responsible for convergence issues in the optimisation process. In addition, the hypothesis of uncorrelated uncertainty of the observation (i.e., $R$ being a diagonal matrix) becomes questionable when sensors
<table><caption>Table 2: Summary of the DA runs performed for the sensitivity analysis.</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<sub>min</sub><sup>*</sup></th><th>&Delta;y<sub>max</sub><sup>*</sup></th><th>L<sub>x</sub>/h</th><th>L<sub>z</sub>/h</th><th>u<sub>x</sub> = 0</th><th>&Delta;C<sub>f</sub></th></tr></thead><tbody><tr><td>R-DNS-BF</td><td>1024 &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>&mdash;</td><td>&mdash;</td></tr><tr><td>DNS-IBM-DA</td><td>128 &times; 64 &times; 64</td><td>39.5</td><td>26.3</td><td>1.3</td><td>67.6</td><td>3&pi;</td><td>&pi;</td><td>4096</td><td>0.033</td></tr><tr><td>case 1</td><td>128 &times; 64 &times; 64</td><td>39.5</td><td>26.3</td><td>1.3</td><td>67.6</td><td>3&pi;</td><td>&pi;</td><td>1024</td><td>0.012</td></tr><tr><td>case 2</td><td>64 &times; 128 &times; 32</td><td>79.0</td><td>52.7</td><td>1</td><td>30</td><td>3&pi;</td><td>&pi;</td><td>1024</td><td>0.149</td></tr><tr><td>case 3</td><td>64 &times; 128 &times; 32</td><td>79.0</td><td>52.7</td><td>1</td><td>30</td><td>3&pi;</td><td>&pi;</td><td>4096</td><td>0.195</td></tr><tr><td>case 4</td><td>64 &times; 128 &times; 32</td><td>79.0</td><td>52.7</td><td>1</td><td>30</td><td>3&pi;</td><td>&pi;</td><td>256</td><td>0.154</td></tr></tbody></table>

are close. On the other hand, if sensors are far from the other in terms of characteristic scales of the flow, the global lack of information may preclude a satisfying state estimation. It must also be stressed that in this case, only one piece of information, namely the streamwise velocity, has been observed for each sensor. The global picture is arguably going to be more complex when multiple physical quantities are provided at each location.

The analysis of second-order statistical moments in figure 11(c)–(f) shows that no differences are observed among the simulations with the same grid. This may imply that the behaviour of such statistics is less sensitive to variations of the friction coefficient/friction velocity in the range investigated here. One can also see that, even if every DA run improves the performance of the underlying IBM, only the realisations with the first grid refinement (DNS-IBM-DA and case 1) are actually able to comply with the indicated level of confidence (0.95 $\langle C_{f,R-DNS-BF} \rangle$).

At last, the optimised values for tensor $D$ and the IBM volume forcing $f_P$ are shown in figures 12 and 13. For the latter, which has been normalised over $U_{c,DNS}/t_{A,DNS}$, the representation includes the average value as well as its standard deviation. Here, results are selected for the most accurate realisation for each grid, i.e., cases 1 and 2. In addition, the volume forcing for the penalisation IBM without DA, using the grid for case 1, is included in figure 13(d)–(f). First of all, the parameter distributions obtained for the two meshes in figure 12 are very similar. This observation confirms the robustness of the optimisation procedure. In addition, for each component, the coefficients for the mesh elements of the interface $\Sigma_b$ play a predominant role when compared to the ones in the solid region $\Omega_b$. This result is in agreement with the properties of continuity and conservation that one can find in Roma's function (see equation C.4) used in the discrete IBM. However,
<figure><img src="image_12.png" /><figcaption>Figure 10: Location of the sensors where $u_x = 0$ is infused during the DA procedure. Snapshots are shown for $y = 0$.</figcaption></figure>
two important points need discussion:

* One would expect the order of magnitude of the parameters linked to $D_{yy}$ to be similar to those characterising $D_{zz}$ and, at the same time, to be significantly smaller than the coefficients defining $D_{xx}$. The results of the DA optimisation indicate that the values for $D_{yy}$ are, on the other hand, similar to those for $D_{xx}$, while the coefficients for $D_{zz}$ are four orders of magnitude lower. This apparently strange result is due to the low sensitivity of the global solution to the intensity of the forcing in the normal and spanwise directions in the layer below the wall of the interface region (i.e., the third one out of the five considered). On the other hand, the optimised value by DA for the layer at $y = 0$ is essential to obtain an accurate prediction of the flow. In particular, the EnKF provides small values of the coefficients for $y = 0$ ($D_{yy} \approx 10^4$, $D_{zz} \approx 10^2$). It has been tested that if higher values of the coefficients $D_{yy}$ and $D_{zz}$ are imposed at the wall, the accuracy is diminished. When optimised values for $D_{yy}$ and $D_{zz}$ for $y = 0$ are imposed, the distribution of the coefficient $D_{xx}$ is the most important IBM ingredient here to obtain an improved prediction of the wall shear stress and to enforce the no-slip condition at the wall. This result is of general interest for the analytic formulation of volume-forcing-based IBM. In fact, one can see that while a shape similar to Roma's delta function is obtained for $D_{xx}$, this is not necessarily true for the other two components. Considering that the DA optimisation showed that the important element is the value at $y = 0$ for $D_{yy}$ and $D_{zz}$, one could arguably use a delta function similar to Roma's proposal for these parameters, however restricting the value of the peak for $y = 0$ to the value optimised via DA.

* The optimised values for the coefficients $D_{xx}$, $D_{yy}$ and $D_{zz}$ produce important differences in the volume forcing when compared with the penalisation IBM without DA. For the streamwise direction, the comparison of figure 13(a) and (d) shows that the DA-augmented IBM exhibits a significantly stronger volume forcing for $y = 0$, around 50 times larger. An opposite trend is obtained for the most external layer of the interface region, where the intensity of $\langle f_{P_x} \rangle$ for the data-driven IBM goes to zero. The comparison of the forcing in the normal direction and the spanwise direction indicates that differences are here less pronounced. However, one can see that the standard deviation of the
<figure><img src="image_13.png" /><figcaption>Figure 11: Comparison of the main statistical moments of the velocity field for the data-driven cases. Results are shown for simulations (—) DNS-IBM-DA, (···) case 1, (—) case 2, (•—) case 3, (···) case 4, and (—) R-DNS-BF.</figcaption></figure>

forcing obtained via DA is occasionally lower, in particular for $y = 0$ and for $\langle f_{P_z} \rangle$.

The generalizability of the present DA approach is now discussed. For the sake of simplicity, the EnKF has been combined with a penalisation method. This choice reduces the complexity of the algorithms, and it improves the clarity of the presentation of all the numerical tasks performed during the optimisation process. However, the EnKF can be combined with every IBM driven by free coefficients. For example, if one considers the discrete IBM by Pinelli et al. [46], one could consider the local optimisation of the free coefficient $a_s$ in equation C.1 or free coefficients or constraints used in the delta functions.

### 4.3. Comparison of computational resources required

Results presented in §4 have shown how the physics-infused procedure relying on the EnKF provides a significant improvement in the predictive capabilities of a classical penalisation IBM model. In this section, the computational resources, i.e., scalar hours CC required to perform each of the simulations of the database in table 1 are discussed. While the simulations
<figure><img src="image_14.png" /><figcaption>Figure 12: Coefficients of tensor $D$ for the simulations (---) case 1 and (-) case 2. Similar trends are observed in cases with identical grid refinement.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: IBM volume forcing term $\langle f_P \rangle$ in the proximity of the bottom wall for (top row) the simulations (---) case 1, and (-) case 2, and (bottom row) the simulation (---) DNS-IBM-CF. The components of the volume forcing are shown for: (left column) the streamwise direction $x$, (central column) the normal direction $y$, and (right column) the spanwise direction $z$. The lines correspond to the mean values, while the shaded areas indicate the interval of the standard deviation.</figcaption></figure>
have not been run on the same machine—therefore, an exact comparison
of the computational resources used cannot be performed—the elements of
discussion provided aim to obtain insights about the efficiency of the data-
infused procedure, i.e., gain in accuracy versus increase in costs. To this pur-
pose, the resources required to perform the simulations are normalised over
$CC_{\text{DNS-BF}}$, which is the total costs employed to run the simulation DNS-BF
($CC^* = CC/CC_{\text{DNS-BF}}$). The numerical simulations are run with the same
numerical schemes and grid (except for the simulation R-DNS-BF and the
data-driven cases 2, 3 and 4). Therefore, differences observed in computa-
tional terms are not associated with the numerical solver.

The simulation DNS-IBM-CF exhibits very similar requirements when
compared with the coarse-grained DNS-BF: $CC^* = 1.027$ since the calcu-
lation of the penalization term is straightforward in this case and does not
require a significant amount of supplementary resources. On the other hand,
the value of the parameter $CC^* = 1.172$ for the simulation DNS-IBM-DF is
relatively higher. The reason is due to the communication and execution of
the interpolation and spreading steps, which are performed by two different
libraries. While this algorithmic cost can be improved with efficient coding,
it still represents a non-negligible increase in the computational costs since,
additionally, 0.74 scalar hours are required when launching the simulation to
load the stencils of the 16 384 Lagrangian markers (one Lagrangian marker
for each Eulerian cell located at $y = 0, 2h$).

Now, the simulation R-DNS-BF is considered. For this simulation, one
needs to take into account that the physical domain is 2 times larger in
the streamwise $x$ and spanwise $z$ directions, and the grid resolution is 4
times more refined in every direction, which implies that the total number of
mesh elements is 256 larger when compared with the grid used for the other
calculations. In addition, the time step is also five times smaller, which means
that the increase of computational resources is around $O(10^4)$ larger than the
simulation DNS-BF. Practical estimation provides a result of $CC^* \approx 12\,000$.

At last, the data-driven runs are considered. In this case, the compu-
tational costs are normalised by the ensemble number $m + 1$ (one addi-
tional CPU core is employed to perform the analysis phases) to represent
the physical/real-world time. It is noticed from the Algorithm 1 that the
cost associated with the matrix operations depends on the number of ob-
servations employed $n_o$, the number of ensemble members $m$, the number
of parameters $\theta$, and the number of degrees of freedom in our system $n$.
In this work, $m$ and $\theta$ are constant for every DA realisation. Therefore,
the physical time related to the data-driven procedures can be expressed as $CC^* = f(n, n_o)$. In figure 14(a), the physical times of the five DA runs analysed in §4.2 are plotted against the two variables $(n^*, n_o^*)$ adimensionalised over their maximum value $(n^{max} = 172.032$ and $n_o^{max} = 4.096)$. One can see that the five results for $CC^*$ from the available DA realisations are well-fitted by a log-level regression model. This is proved by estimating the coefficient of determination $R^2$. For a dependent variable $z$ whose predicted values are $\hat{z}_i$, its definition is $R^2 = \sum_i (\hat{z}_i - \bar{z})^2 / \sum_i (z_i - \bar{z})^2$, which represents the proportion of the total variation in the $z_i$ values from the average $\bar{z}$ explained by the regression equation. For the predicted $\widehat{CC}^*$, the following equation is derived with a $R^2 = 99.902\%$:

$$ \ln \left( \frac{\widehat{CC}^*}{m+1} \right) = 1.651 + 1.448 n^* + 0.511 n_o^* \quad (11) $$

It is observed that, in the empirical relation expressed in equation 11, $n^*$ is more important than $n_o^*$ in the range of investigation of the present study. The data-driven runs are now sorted by increasing the computational time. Case 4 provides a $CC^*/(m+1) = 8.797$, case 2 is 9.284, case 3 is 13.708, case 1 is 24.797, and DNS-IBM-DA is 37.598. This result clearly indicates that the hyperparameters driving the DA technique are critical both in terms of accuracy as well as for computational costs. Further investigations should be focused on determining these coefficients to optimise the algorithm, possibly by combining Machine Learning (ML) techniques altogether. Besides, as already expressed in §4.2, one can see that more observations do not necessarily imply better accuracy. To quantify such a statement, the parameter $\Delta_{\widehat{C_f}}$ from table 2 is also approximated by a log-level regression model and shown in figure 14(b). In this case, $R^2 = 99.269\%$ for the range covered in these studies. The resulting empirical formula is:

$$ \widehat{\Delta_{C_f}} = 0.219 - 0.223 n^* + 0.041 n_o^* \quad (12) $$

The positive coefficient of $n_o^*$ in equation 12 is an indication of the potential degradation of the accuracy of the DA technique due to overconstraints, such as for the DA run case 3 in the present work.

## 5. Conclusions

A physics-infused strategy has been developed to improve the accuracy of a classical Immersed Boundary Method (IBM) tool, namely the penali-
<figure><img src="image_16.png" /><figcaption>Figure 14: Real (red dots) and predicted (grey plane) values for (a) computational time and (b) the relative error of the friction coefficient $\Delta \overline{C_f}$ for the range covered in these studies.</figcaption></figure>

sation method. Physical information about the flow state, in the form of no-slip condition and stress at the wall, has been introduced using an online data-driven technique based on the Ensemble Kalman Filter (EnKF). This sequential tool has been used to update the physical state of the flow, as well as to optimise the parametric description of the penalisation IBM integrated into the dynamic equations. The analysis of the statistical quantities of the flow, which includes $u_\tau$, the mean streamwise velocity $\langle u_x \rangle$ and the components of the Reynolds stress tensor, indicates that the DA procedure is able to improve the accuracy of the prediction of the flow field, significantly reducing the discrepancy with a high-resolution DNS. In addition, the analysis of the instantaneous features, in the form of isocontours of the Q-criterion, shows that a similar organisation of the structures is obtained despite the different grid resolution. The sensitivity of the physics-infused procedures to variations in the model and the observation has then been analysed. For the former, a grid with a different resolution in each space direction has been selected. For the latter, the position and number of sensors have been modified. The results indicate a mild sensitivity to these key elements, lower than the one exhibited for classical simulations with the same mesh changes. Therefore, the procedure shows good robustness characteristics. In addition, the current version of the online DA procedure is highly competitive in terms of accuracy vs computational costs required, and it can be further improved using advanced EnKF strategies such as multiple local EnKFs with extended localisation, as performed in environmental studies [11]. The critical point
here is that the proposed technique partially bypasses the need to produce,
store, and manipulate extensive databases, reducing the computational re-
sources needed for these means.

While ongoing results are restricted to the analysis of an attached bound-
ary layer with zero external pressure gradient, the present analysis indicates
that the EnKF can be efficiently used to optimise online the performance of
an IBM when streaming data is forwarded to the numerical model. This find-
ing, which is important for future applications of IBM within digital twins,
needs further testing using more sophisticated algorithms for complex test
cases. For example, the team is currently investigating the flow in a centrifu-
gal pump [69], where the rotating impeller is modelled using IBM. The DA
algorithm combines such a CFD model with experimental PIV data from the
velocity field obtained using local experimental facilities. Here, significant
challenges arise in accurately combining turbulence modelling and IBM in
near-wall regions.

However, the research work here presented clearly indicates how building
accurate digital twins relying on numerical simulation and streaming data
from experiments is an open challenge which demands advancement in many
different scientific topics. Computational costs are a central key issue, as nu-
merical simulation's elapsed time needed to run the code and provide results
must be shorter than the evolution time of the physical flow. In addition,
DA algorithms also have important computational demands, which can how-
ever been reduced via advanced localisation and clipping techniques [70]. A
third critical area of investigation deals with the number and positioning of
sensors as well as the quality of the sampled data needed to obtain accurate
and converged DA predictions. These three elements are not controllable
for every application due, for example, to difficulties in sensor positioning
in critical regions or structural limitations in frequency acquisition. All of
these elements discussed here highlight how disruptive advancement in re-
search needs to be achieved in the coming decades to achieve such ambitious
technical targets for digital twins in fluid mechanics applications.

Acknowledgements

This project was provided with computer and storage resources by GENCI
at TGCC thanks to the grant 2022-A0122A01741 on the supercomputer Jo-
liot Curie's SKL partition. The present research work has been developed in
the framework of the project ANR-JCJC-2021 IWP-IBM-DA.
Appendix A.

The Ensemble Kalman Filter (EnKF) is presented here. This method relies on an ensemble of Monte-Carlo realisations to account for the uncertain behaviour of the model. To this purpose, an ensemble of state vectors $\mathbf{u}_i$ ($i \in [1, m]$, where $m$ is the number of ensemble members) is used to propagate the uncertainty of the model $\mathcal{M}$ during the forecast ($f$) step from $k-1$ to $k$:

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

When observation is available in time, an analysis (a) phase is performed, where the DA procedure updates the system state and the model. The model realisations are combined to obtain an anomaly matrix $\mathbf{X}$, which measures the deviation of the ensemble with respect to the mean $\bar{\mathbf{u}}$. By using the hypothesis of statistical independence of the realisations, an ensemble approximation $\mathbf{X}$ of the error covariance matrix can be obtained. Assuming that the analysis phase is performed at the time step $k$ as before, the determination of such matrix is obtained via the following equations:

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

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

The states $\mathbf{u}$ and $\bar{\mathbf{u}}$ are manipulated as vectors of matrix size $[n, 1]$, where $n$ is the number of degrees of freedom characterising the random process. For a CFD application for incompressible flows, the most classical choice is to consider the full velocity field discretised over the $n_{GE}$ grid elements. Therefore, $n = 3 n_{GE}$. Similarly, the size of the matrix $\mathbf{X}$ is $[n, m]$. In order to obtain a well-posed numerical system, an ensemble of $n_o$ observations is also obtained via a perturbation of the observation vector $\mathbf{y}_k$ whose size is $[n_o, 1]$. This procedure allows obtaining an observation matrix $\mathbf{Y}_k$ of size $[n_o, m]$, whose $m$ columns are defined as $\mathbf{Y}_{i,k} = \mathbf{y}_k + \mathbf{e}_i$, $i = 1, 2, ..., m$. The random noise $\mathbf{e}_i$ is described by a Gaussian probability function $\mathbf{e}_i \sim \mathcal{N}(0, \mathbf{R}_k)$. As for the KF, $\mathbf{R}_k$ represents the observation covariance matrix.

The EnKF requires the use of a sampling matrix $\mathcal{H}(\mathbf{u}_k^f)$ to project the model into the position of the observations. Analogously to the system's state $\mathbf{u}_k^f$ in equation A.3, it can be decomposed into a matrix $\mathbf{S}_k^f$ where
each column represents a normalised anomaly, and the mean is defined as $\overline{\mathcal{H}(\mathbf{u}_k^f)} = \sum_{i}^{m} \mathcal{H}(\mathbf{u}_{i,k}^f)/m$:

$$ \left[ \mathbf{S}_k^f \right]_i = \frac{\mathcal{H}(\mathbf{u}_{i,k}^f) - \overline{\mathcal{H}(\mathbf{u}_k^f)}}{\sqrt{m-1}} \quad (A.4) $$

For an infinite ensemble size, $\lim_{m \to +\infty} \mathbb{E}[(e - \mathbb{E}[e])(e - \mathbb{E}[e])^T] = \mathbf{R}$, which becomes time-independent, and the Kalman gain matrix $\mathbf{K}_k$ describing the most optimal correlation between the state and the observations is simplified to [71, 72]:

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

All the elements presented constitute the essential ingredients to calculate the updated system's state $\mathbf{u}_k^a$, which is then used for the time advancement of the numerical model $\mathcal{M}$ in the following iteration $k + 1$:

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

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

## Appendix B.

The data-driven algorithm used to augment the IBM relies on numerical solvers based on the PISO algorithm [41]. The velocity and pressure field are iteratively updated to comply with the momentum equation and the Poisson equation. When the penalisation volume force is velocity-dependent $\mathbf{f}_P(\mathbf{u}_{t,j})$, it may be solved implicitly in equation 2. Let us consider the time $t$ and the time step advancement $\Delta t$. Considering $j \in [1, J]$ as a single iteration from the PISO loop, the following steps are performed:

1. Resolution of the momentum equation, from where $\mathbf{u}_{t,j}$ is obtained.

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

2. Estimation of the pressure $p_{t,j+1}$ through the Poisson equation. $A'$ is a scalar field calculated from $\mathbf{A}$ and $\mathbf{f}_P$, whereas $\mathbf{T}'(\mathbf{u}_{t,j})$ is a tensor field containing the discretised form of all the terms on the left side of equation B.1.
**Algorithm 1:** Scheme of the EnKF used in the present study.

**Input:** $\mathcal{M}$, $\mathbf{R}$, $\mathbf{y}_k$, and a prior/initial state system $\mathbf{u}_{i,0}$, where
usually $\mathbf{u}_{i,0} \sim \mathcal{N}(\bar{\mathbf{u}}_0, \sigma_0^2)$

for $k = 1, 2, ..., K$ do

for $i = 1, 2, ..., m$ do

1

Advancement in time of the state vectors:

$$ \mathbf{u}_{i,k}^{f} = \mathcal{M}_{k:k-1} \mathbf{u}_{i,k-1} $$

2

Generation of an observation matrix from the confidence level
given to the observation data:

$$ \mathbf{Y}_{i,k} = \mathbf{y}_k + \mathbf{e}_i, \text{ with } \mathbf{e}_i \sim \mathcal{N}(0, \mathbf{R}) $$

3

Estimation of the ensemble means (system state and
projection matrix):

$$ \bar{\mathbf{u}}_{k}^{f} = \frac{1}{m} \sum_{i=1}^{m} \mathbf{u}_{i,k}^{f}, \quad \overline{\mathcal{H}(\mathbf{u}_{k}^{f})} = \frac{1}{m} \sum_{i=1}^{m} \mathcal{H}(\mathbf{u}_{i,k}^{f}) $$

4

Computation of the anomaly matrices (system state and
projection matrix):

$$ \mathbf{X}_k = \frac{\mathbf{u}_{i,k} - \bar{\mathbf{u}}_k}{\sqrt{m-1}}, \quad \mathbf{S}_k = \frac{\mathcal{H}(\mathbf{u}_{i,k}^f) - \overline{\mathcal{H}(\mathbf{u}_k^f)}}{\sqrt{m-1}} $$

5

Calculation of the Kalman gain (with covariance localisation):

$$ \mathbf{K}_k = \mathbf{L} \odot \mathbf{X}_k^f (\mathbf{S}_k)^T [\mathbf{S}_k (\mathbf{S}_k)^T + \mathbf{R}]^{-1} $$

6

Update of the state matrix:

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

7

Application of (deterministic) covariance inflation:

$$ \mathbf{u}_{i,k}^{a} = \bar{\mathbf{u}}_{k}^{a} + \lambda_{k} (\mathbf{u}_{i,k}^{a} - \bar{\mathbf{u}}_{k}^{a}) $$
$$ \nabla \cdot \frac{1}{A'} \nabla p_{t,j+1} = \nabla \cdot \left( \frac{T'(u_{t,j})}{A'} \right) \qquad (B.2) $$

3. Update of the velocity field $u_{t,j+1}$ to satisfy the zero-divergence condition.

$$ u_{t,j+1} = \frac{T'(u_{t,j})}{A'} - \frac{1}{A'} \nabla p_{t,j+1} \qquad (B.3) $$

Equations B.2 and B.3 are solved iteratively until reaching convergence. In the case of the presence of high-fidelity data at time $t$, when $j = J$, $u_{t,j} = u_k^f$, which is the forecast velocity field at the $k^{th}$ analysis phase for the EnKF. Together with the coefficients $\theta_k^f$ of the tensor $D$, it constitutes the forecast system's state. The EnKF updates them to estimate $(u_k^a, \theta_k^a)$. This state, however, does not necessarily respect the Navier-Stokes equations, as the velocity field is updated while the pressure field is the one obtained in the forecast. In order to obtain a consistent final solution, an additional update of the pressure is performed via another resolution of the Poisson equation.

1. The volume force term $(f_P)_k^a$ is updated by using the new coefficients $\theta_k^a$ from the tensor $D$:

$$ (f_P)_k^a = -\nu D(\theta_k^a) u_k^a \qquad (B.4) $$

2. A Poisson equation is solved to update the pressure $p_k^a$. The complete flow field and the updated parameters $\theta_k^a$ are then used to solve the equation B.1 in $t + \Delta t$.

$$ \nabla^2 p_k^a = -\nabla \cdot (u_k^a \nabla u_k^a) + \nabla \cdot (f_P)_k^a \qquad (B.5) $$

## Appendix C.

The complete algorithm for the reference discrete IBM [40] is described below. The method relies on two steps, which are now detailed:

1. In the interpolation step, the physical variables describing the flow, which are calculated on the Eulerian mesh, are interpolated on the Lagrangian Markers. The source term is then calculated using the Lagrangian space. Its structural form is very similar to the one seen for the penalisation method in equation 3:
$$F_P = -a_s (U^* - U_{ib}) \qquad (C.1)$$

Capital letters are used to indicate the variables in the Lagrangian space. $U^*$ is the velocity field interpolated from the Eulerian grid while $a_s$ is a coefficient resulting from the discretization procedure. $U_{ib}$ is homologous to $u_{ib}$ in the Lagrangian framework.

2. After the source term is calculated on the Lagrangian space, it is projected on the Eulerian grid during the spreading step. This procedure allows obtaining a closed expression for the source term to be integrated within the computational solver.

For each time $t$ of the numerical simulation, the steps previously presented are integrated following this algorithm:

1. Estimation of an initial velocity $u_{t,j^*}$ following the expression (C.2) for the momentum predictor without a source term.

$$A u_{t,j^*} - b = -\nabla p_{t-\Delta t} \qquad (C.2)$$

2. Interpolation (I) of the velocity $u_{t,j^*}$ from the element $i$ of a subspace $D_s$ (Eulerian mesh) to the Lagrangian marker $s$, with $s \in [1, N_s]$.

$$I[u_{t,j^*}]_s = [U^*](s) = \sum_{i \in D_s} (u_{t,j^*})_i \delta_d(x_i - \mathcal{X}_s) \Delta x \qquad (C.3)$$

$x_i$ and $\mathcal{X}_s$ refer to the position in the Eulerian and Lagrangian frameworks, respectively. $\Delta x$ describes the Eulerian quadrature ($\Delta x = \Delta x \Delta y \Delta z$ in the Cartesian frame of coordinates) and the interpolation kernel $\delta_d$ is a discretised delta function based on the Euclidean distance $r = (x_j - \mathcal{X}_s)/d$ proposed by Roma et al. [73]:

$$\delta_d = \begin{cases} \frac{1}{3}(1 + \sqrt{-3r^2 + 1}) & 0 \le r \le 0.5 \\ \frac{1}{6}(5 - 3r - \sqrt{-3(1-r)^2 + 1}) & 0.5 \le r \le 1.5 \\ 0 & \text{otherwise} \end{cases} \qquad (C.4)$$

The kernel is multiplied with a prescribed function $\mathcal{G}$ to account for the directional stretching of the mesh elements.
3. Calculation of $F_P$ on the $N_s$ Lagrangian markers by employing the equation C.1 already presented in §2.2. The target values over the $N_s$ Lagrangian markers are $U_{ib} = 0$ if the body surface is not moving.

4. Spreading step to obtain the Eulerian counterpart $f_P$:

$$f_P(\mathbf{x}_i) = \sum_{i \in D_s} (F_P)_k \delta_d(\mathbf{x}_i - \mathcal{X}_k) \epsilon_k \quad (C.5)$$

The k-index denotes a loop over the Lagrangian markers whose support contains the Eulerian node j, and $\epsilon_k$ is the Lagrangian quadrature.

$$C\epsilon = 1 \quad (C.6)$$

where the vectors $\epsilon = (\epsilon_1, ..., \epsilon_{N_s})^T$ and $1 = (1, ..., 1)^T$ with size $N_s$. C is the matrix defined by the product between the $k^{th}$ and $l^{th}$ interpolation kernels:

$$C_{kl} = \sum_{i \in D_s} \delta_d(\mathbf{x}_i - \mathcal{X}_k) \delta_d(\mathbf{x}_i - \mathcal{X}_l) \quad (C.7)$$

5. Repetition of the momentum predictor, but with the inclusion of $f_P$, as shown in equation B.1, to obtain $\mathbf{u}_{t,j}$. However, different from the continuous IBM, in this case, $f_P \neq f_P(\mathbf{u}_{t,j})$:

$$A \mathbf{u}_{t,j} - \mathbf{b} - f_P = -\nabla p_{t-\Delta t} \quad (C.8)$$

6. Computation of the Poisson equation and the momentum corrector (equations C.9 - C.10) iteratively until achieving convergence. A is a scalar field obtained from $A$, and $T(\mathbf{u}_{t,j})$ contains the information of all the terms on the left side of equation C.8.

$$\nabla \cdot \frac{1}{A} \nabla p_{t,j+1} = \nabla \cdot \left( \frac{T(\mathbf{u}_{t,j})}{A} \right) \quad (C.9)$$

$$\mathbf{u}_{t,j+1} = \frac{T(\mathbf{u}_{t,j})}{A} - \frac{1}{A} \nabla p_{t,j+1} \quad (C.10)$$

## References

[1] C. S. Peskin, Flow patterns around heart valves, J. Comp Phys. 10 (1972) 252–271.
[2] C. S. Peskin, The immersed boundary method, Acta Numerica 11 (2002) 1–39.

[3] R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261.

[4] W. Kim, H. Choi, Immersed boundary methods for fluid-structure interaction: a review, International Journal of Heat and Fluid Flow 75 (2019) 301–309.

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

[6] S. B. Pope, Turbulent flows, Cambridge University Press, 2000.

[7] P. Sagaut, Large-eddy simulation for incompressible flows. An introduction, Springer-Verlag, 2005.

[8] D. C. Wilcox, Turbulence Modeling for CFD, 3rd Edition, DCW Industries Inc., La Canada CA, 2006.

[9] D. Simon, Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches, Wiley, 2006.

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

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

[12] K. Duraisamy, H. Xiao, G. Iaccarino, Turbulence modeling in the age of data, Annual Review of Fluid Mechanics 51 (2019) 1–23.

[13] 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.

[14] W. N. Edeling, P. Cinnella, R. P. Dwight, H. Bijl, Bayesian estimates of parameter variability in the $K-ε$ turbulence model, Journal of Computational Physics 258 (2014) 73–94.
[15] L. Margheri, M. Meldi, M. V. Salvetti, P. Sagaut, Epistemic uncertainties in RANS model free coefficients, Computers & Fluids 102 (2014) 315-335.

[16] 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.

[17] 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.

[18] 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.

[19] 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, Phys. Rev. Fluids 6 (2021) 064607. doi:10.1103/PhysRevFluids.6.064607.

[20] 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.

[21] 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.

[22] 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.

[23] 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.

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

[26] 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, submitted, available on ArXiv 2212.13831 (2022) 1-25.

[27] X. Fan, J. Wang, Differentiable hybrid neural modeling for fluid-structure interaction, Journal of Comp. Physics 496 (2024) 112584.

[28] R. Sundar, D. Majumdar, D. Lucor, S. Sarkar, Physics-informed neural networks modelling for systems with moving immersed boundaries: Application to an unsteady flow past a plunging foil, Journal of Fluids and Structures 125 (2024) 104066. doi:https://doi.org/10.1016/j.jfluidstructs.2024.104066.

[29] P. Angot, C. Bruneau, P. Fabrie, A penalization method to take into account obstacles in incompressible viscous flows, Numerische Mathematik 81 (1999) 497-520.

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

[31] 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.

[32] 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.

[33] 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.

[34] 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.
[35] 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.

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

[37] G. Tabor, M. Baba-Ahmadi, Inlet conditions for large eddy simulation: A review, Computers & Fluids 39 (4) (2010) 553–567.

[38] M. Meldi, M. V. Salvetti, P. Sagaut, Quantification of errors in large-eddy simulations of a spatially evolving mixing layer using polynomial chaos, Physics of Fluids 24 (2012) 035101.

[39] B. Selma, M. Désilets, P. Proulx, Optimization of an industrial heat exchanger using an open-source CFD code, Applied Thermal Engineering 69 (1–2) (2014) 241–250.

[40] E. Constant, F. Javier, M. Meldi, P. Meliga, E. Serre, An immersed boundary method in OpenFOAM: Verification and validation, Comp. Fluid. 157 (2017) 55–72.

[41] 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.

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

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

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

[45] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comp. Phys. 209 (2) (2005) 448–476.

[46] A. Pinelli, I. Z. Naqavi, U. Piomelli, J. Favier, Immersed-boundary methods for general finite-difference and finite-volume navier-stokes solvers, J. Comp Phys. 229 (24) (2010) 9073–9091.
[47] W. K. Liu, S. Jun, M. Y. F. Zhang, Reproducing kernel particle methods, Int. J. Numer. Meth. Fluid. 20 (8) (1995) 1081-1106.

[48] S. B. Daley, Atmospheric Data Analysis, Cambridge University Press, 1991.

[49] Z. Sirkes, E. Tziperman, Finite Difference of Adjoint or Adjoint of Finite Difference?, Monthly Weather Review 125 (1997) 3373-3378.

[50] G. Artana, A. Camilleri, J. Charlier, E. Memin, Strong and weak constraint variational assimilations for reduced order fluid flow modeling, Journal of Computational Physics 231 (2014) 3264-3288.

[51] D. P. G. Foures, N. Dovetta, D. Sipp, P. J. Schmid, A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean flow reconstruction, Journal of Fluid Mechanics 759 (2014) 404-431.

[52] V. Mons, O. Marquet, Linear and nonlinear sensor placement strategies for mean flow reconstruction via data assimilation, Journal of Fluid Mechanics 923 (2021) A1.

[53] M. Katzfuss, J. R. Stroud, C. K. Wickle, Understanding the ensemble kalman filter, The Amer. Stat. 70 (4) (2016) 350-357.

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

[55] J. W. Labahn, H. Wu, B. Coriton, J. H. Frank, M. Ihme, Data assimilation using high-speed measurements and LES to examine local extinction events in turbulent flames, Proceedings of the Combustion Institute 37 (2019) 2259-2266.

[56] X. Zhang, C. Michelin-Ströfer, H. Xiao, Regularized ensemble Kalman methods for inverse problems, Journal of Computational Physics 416 (2020) 109517.

[57] A. F. C. da Silva, T. Colonius, Flow state estimation in the presence of discretization errors, Journal of Fluid Mechanics 890 (2020) A10.

[58] M. L. Provost, J. Eldredge, Ensemble kalman filter for vortex models of disturbed aerodynamics flows, Phys. Rev. Fluid. 6 (2021) 050506.
[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] G. Moldovan, G. Lehnasch, L. Cordier, M. Meldi, A multigrid/ensemble Kalman filter strategy for assimilation of unsteady flows, Journal of Computational Physics 443 (2021) 110481.

[62] A. Reflox, B. Courbet, A. Murrone, P. Villedieu, C. Laurent, Cedre software, Aero. Lab 2 (2011) 1-10.

[63] M. Meldi, A. Poux, A reduced order model based on Kalman Filtering for sequential Data Assimilation of turbulent flows, Journal of Computational Physics 347 (2017) 207-234.

[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] S. Butterworth, On the theory of filter amplifiers, Experimental Wireless and the Wireless Engineer 7 (1930) 536-541.

[67] S. Kamm, S. S. Veekati, T. Müller, N. Jazdi, M. Weyrich, A survey on machine learning based analysis of heterogeneous data in industrial automation, Computers in Industry 149 (2023) 103930. doi:https://doi.org/10.1016/j.compind.2023.103930.

[68] K. Taira, T. Colonius, The immersed boundary method: A projection approach, Journal of Computational Physics 225 (2) (2007) 2118-2137. doi:https://doi.org/10.1016/j.jcp.2007.03.005.

[69] A. Dazin, G. Cavazzini, G. Pavesi, P. Dupont, S. Coudert, G. Ardizzon, G. Caignaert, G. Bois, High-speed stereoscopic piv study of rotating instabilities in a radial vaneless diffuser, Experiments in Fluids 51 (2011) 83-93.
[70] L. Villanueva, Développement d'outils d'assimilation de données pour l'estimation augmentée d'écoulements internes, Ph.D. thesis, ENSMA, Poitiers (2024).

[71] 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).

[72] 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.

[73] A. M. Roma, C. S. Peskin, M. J. Berger, An adaptive version of the immersed boundary method, J. Comp. Phys. 153 (1999) 509-534.