# Data Assimilation and Uncertainty Quantification for the RANS prediction of separated flows

Maxime Casanova, Vincent Mons, Pedro Stefanin Volpiani, O. Marquet, Denis Sipp

► To cite this version:

Maxime Casanova, Vincent Mons, Pedro Stefanin Volpiani, O. Marquet, Denis Sipp. Data Assimilation and Uncertainty Quantification for the RANS prediction of separated flows. ETMM15, Sep 2025, Dubrovnik, Croatia. pp.262-267, 10.5281/zenodo.17280519. hal-05384351

HAL Id: hal-05384351

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

Submitted on 26 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# DATA ASSIMILATION AND UNCERTAINTY QUANTIFICATION FOR THE RANS PREDICTION OF A SEPARATED FLOW

M. Casanova¹, V. Mons¹, P.S. Volpiani¹, O. Marquet¹ and D. Sipp¹

¹ DAAA, ONERA, Institut Polytechnique de Paris, 92190, Meudon, France

maxime.casanova@onera.fr

## Abstract

This work presents a formalism for enhancing RANS models by data assimilation, that incorporates an assessment of the quality of the assimilated results through uncertainty quantification. The framework is applied to the two-dimensional separated flow in a converging-diverging channel with a Reynolds number of 12,600. The Spalart-Allmaras (SA) model is selected for its relative simplicity compared to other RANS models. The methodology is evaluated using sparse mean velocity measurements that are used to infer among two types of corrective terms in the RANS-SA equations: one acting on the momentum equations, and the other on the modified eddy-viscosity equation, via the turbulent production term. Results show that, in the case of sparse data, correcting the turbulent production term yields more accurate posterior mean flow reconstructions. In this case, the UQ strategy, based on a low-rank approximation of the Hessian operator, successfully predicts posterior standard deviations that align with the true residual errors, highlighting regions of inaccurate flow reconstruction.

## 1 Introduction

Turbulent flows, characterized by their chaotic and non-linear dynamics, present significant challenges in engineering design. While Computational Fluid Dynamics (CFD) solvers are essential instruments for accurate flow analysis and prediction, high-fidelity simulations come at an often prohibitive computational cost for industrial applications. Reynolds-average Navier-Stokes (RANS) models address this challenge by modeling turbulence scales through additional transport equations, thus reducing the cost of the resulting simulations. Despite the widespread use of RANS models, their applicability is limited by additional assumptions and simplifications, which can introduce inaccuracies. In particular, RANS models have been shown to yield inaccurate predictions for highly separated and transitional flows (G. Iaccarino et al., 2003), underscoring the need for rigorous validation using high-fidelity data. Consequently, data assimilation (DA) and uncertainty quantification (UQ) are essential tools for enhancing the accuracy and reliability of RANS models. Data-assimilation techniques integrate observational

reference data into numerical models to enhance their accuracy. Several DA methods have been successfully applied to RANS models, leading to a significant improvement in prediction accuracy. Adjoint-based methods, in particular, have shown remarkable performance in steady mean-flow reconstruction. Cato et al. (2023) compared different correction approaches by incorporating several corrective terms into the Spalart-Allmaras (SA) RANS equations. The authors applied data assimilation to the same flow setup considered in this study, a separated flow in a converging-diverging channel at a Reynolds number of 12,600. With dense velocity measurements, they showed that forcing applied in the momentum equations yielded the best reconstruction and accurately reproduced the flow separation from the direct numerical simulation (DNS). However, while adjoint-based methods can significantly improve accuracy, they do not inherently provide uncertainty estimates. This limitation becomes critical when dealing with sparse observational data, as it is essential to quantify the uncertainty in reconstructed flow quantities, particularly in regions that have been poorly constrained by the available reference data.

Uncertainty quantification techniques are used to enhance the reliability of model predictions, which are often compromised by model-form uncertainties (e.g., Boussinesq assumption) inherent to RANS modeling. Xiao & Cinnella (2019) provided a comprehensive review of various approaches for quantifying and reducing model-form uncertainties in RANS simulations, addressing both parametric (low-dimensional) and non-parametric (high-dimensional) UQ methods in turbulence modeling. They showed that non-parametric UQ provided a more reliable estimate by focusing directly on the modeled terms in RANS solvers. A large set of methods based on the Laplace approximation (LA), as reviewed in Ghattas & Willcox (2021), allowed the development of low-rank Hessian-based approaches for non-parametric UQ. Those methods have been successfully applied to several fields, including ocean dynamics Kalmikov & Heimbach (2014), and are particularly relevant to the present study.
In this work, adjoint-based DA is used to enhance the accuracy of the Spalart Allmaras (SA) RANS state given sparse mean velocity observations, while the low rank Hessian-based UQ is utilized to estimate the uncertainties related to the assimilated mean state, providing posterior error bounds. To this end, this study is a natural extension of our previous work (Cato et al., 2023), which explored different DA approaches using only dense observations and did not include uncertainty quantification (UQ) to assess the reconstructed mean state.

# 2 Mathematical Framework

## 2.1 Finite element discretization

In the following, we rely on the discretization provided by the finite element method so that the degree of freedom of the SA state $q = \{q_i\}^{n_q}$ and forcing $f = \{f_i\}^{n_f}$ are related to their finite element representation as

$$q(x) = \sum_{i=1}^{n_q} q_i \phi_{q_i}(x), \quad f(x) = \sum_{i=1}^{n_f} f_i \phi_{f_i}(x), \quad (1)$$

where the $\phi_i(x)$ denotes the associated basis functions. The current SA residual is enhanced by two distinct forcings $f_u(x)$ and $f_P(x)$ acting in the momentum equation and on the turbulence production respectively so that $N(q(x), f(x)) =$

$$\left( \begin{array}{l} u(x) \cdot \nabla u(x) + \nabla p(x) - \nabla \cdot (\gamma \nabla_s u(x)) - f_u(x) \\ \nabla \cdot u(x) \\ u(x) \cdot \nabla \nu(x) - \nabla \cdot (\eta \nabla \nu(x)) - s(f_P(x)) \end{array} \right) \quad (2)$$

with $q(x) = [u(x), p(x), \nu(x)]^T$ and $f(x) = \{f_u(x), f_P(x)\}$ and where the source term $s$ is defined as

$$s(f_P(x)) = (1 + f_P(x))P + D + C \quad (3)$$

including the turbulence production, destruction and cross diffusion terms. Note that in the following, the application of DA and UQ is done for one forcing at a time.

## 2.2 Problem Statement

Let us assume we have the discretized SA model $N$ which provide an estimate $q \in \mathbb{R}^{n_q}$ of a high fidelity mean state $q_y \in \mathbb{R}^{n_q}$ of a fluid system and which is a function of an uncertain forcing parameter $f \in \mathbb{R}^{n_f}$ such that

$$N(q, f) = 0. \quad (4)$$

Furthermore, suppose that the forcing $f$ has a prior mean $f_B$ and prior covariance $B_{f_B} \in \mathbb{R}^{n_f \times n_f}$. We consider also we have access to high fidelity observations $y \in \mathbb{R}^{n_y}$ through an observation matrix $Y \in \mathbb{R}^{n_y \times n_q}$ which maps the high fidelity state $q_y$ to the observations $y$ with an additive normally distributed measurement error $\eta \sim N(0, B_y)$

$$y = Yq_y + \eta. \quad (5)$$

with zero mean and covariance matrix $B_y \in \mathbb{R}^{n_y \times n_y}$.

Note that the prior covariance matrices $B_{f_B}$ and $B_y$ represent the confidence in the model $N$ and observations $y$. The prior state covariance $B_{q_B}$ is estimated by propagating $B_{f_B}$ into the prior resolvent matrix $R_B$ as

$$B_{q_B} = R_B B_{f_B} R_B^T, \quad (6)$$

where the resolvent matrix is defined as

$$R = -(\partial_q N)^{-1} \partial_f N, \quad (7)$$

which by construction relates small perturbations of forcing $\delta f$ to induced modifications of the state, $\delta q = R\delta f$.

In a first step, DA is used to determine the posterior (i.e, assimilated) mean forcing $f_A$ and mean state $q_A = q(f_A)$ given the high-fidelity observations $y$ and prior knowledge of $f$. In a second step, UQ is used to estimate the uncertainties related to the posterior means, typically characterized by the knowledge of their posterior variances $\sigma_{f_A}^2$ and $\sigma_{q_A}^2$.

## 2.3 Posterior mean

We present the methodology used to determine the posterior mean forcing $f_A$ and mean state $q_A$. The framework is derived from a general Bayesian formulation. We assume a normal prior forcing distribution $p(f)$ and a likelihood $p(y|q)$

$$p(f) \propto \exp \left( -\frac{1}{2} (f - f_B)^{\mathrm{T}} B_{f_B}^{-1} (f - f_B) \right), \quad (8)$$

$$p(y|q) \propto \exp \left( -\frac{1}{2} (y - Yq)^{\mathrm{T}} B_y^{-1} (y - Yq) \right), \quad (9)$$

so that using the Bayes theorem, we end up with the following posterior forcing distribution

$$p(f|y) \propto \exp(-J(f)), \quad (10)$$

with

$$J(f) = \frac{1}{2} (f - f_B)^{\mathrm{T}} B_{f_B}^{-1} (f - f_B) + \frac{1}{2} (y - Yq)^{\mathrm{T}} B_y^{-1} (y - Yq). \quad (11)$$

The most likely forcing, given the measurements $y$, is denoted as $f_A$; it maximizes the posterior distribution shown in Eq.(10). We remark that this translates to finding the minimum of $J(f)$, which satisfies

$$\frac{dJ}{df}|_{f_A} = \nabla J|_{f_A} = 0. \quad (12)$$

The adjoint approach is used to estimate efficiently the cost functional gradient $\nabla J|_f$. This requires constructing the Lagrangian functional

$$\mathcal{L}(q, q^\dagger, f) = J(q, f) - q^\dagger \mathrm{T} N(q, f), \quad (13)$$
where $q^\dagger$ denotes the adjoint state.

The variation of the Lagrangian $\mathcal{L}$ with respect to the forcing $f$ and state $q$ identifies the gradient $\nabla J|_f$ as

$$ \nabla J|_f = \mathbf{B}_{f_B}^{-1}(\mathbf{f} - \mathbf{f}_B) - \mathbf{R}^\mathrm{T}\mathbf{Y}^\mathrm{T}\mathbf{B}_y^{-1}(\mathbf{y} - \mathbf{Y}\mathbf{q}(\mathbf{f})). \quad (14) $$

In summary, the estimation of $\nabla J|_f$ is based on the resolution of the adjoint system where the posterior mean forcing $f_A$ and mean state $q_A$ are obtained by satisfying the following system of equations

$$ f_A = \arg \min_{f} J(f), \quad \mathcal{N}(q_A, f_A) = 0. \quad (15) $$

## 2.4 Posterior variance

We present the methodology used to determine the posterior variances $\sigma_{f_A}^2$ and $\sigma_{q_A}^2$, respectively, used to characterize the uncertainties related to $f_A$ and $q_A$. Assuming that we have identified a posterior forcing $f_A$, we now search to quantify the associated variance, which requires the ability of sampling the posterior distribution $p(f|y)$. Instead of relying on a Monte-Carlo method which is untractable in the case of high-dimensional forcings, a tractable approach may first consist in performing a second-order Taylor expansion of $J$ around the posterior forcing $f_A$.

$$ J(f) \approx J(f_A) + \frac{1}{2} (f - f_A)^{\mathrm{T}} \mathbf{H}|_{f_A} (f - f_A), \quad (16) $$

where $\mathbf{H}|_{f_A} = \frac{d\nabla J}{df}|_{f_A}$ is the Hessian of $J$ at the assimilated state of the forcing and is symmetric. Note that this approximation is valid in the limit of $f$ around $f_A$.

Accordingly, from Eq.(10), the posterior distribution $p(f|y)$ can be approximated at the assimilated state as

$$ p(f|y) \propto \exp\left(-\frac{1}{2}(f-f_A)^{\mathrm{T}}\mathbf{H}|_{f_A}(f-f_A)\right). \quad (17) $$

Note that Eq.(17) is called the Laplace approximation (LA), extensively discussed in Ghattas & Willcox (2021). This approximation is exact in the case of linear residual $\mathcal{N}$ and normal priors.

By identification and under the LA, the posterior forcing should thus follow a normal distribution with mean $f_A$ and posterior forcing covariance $\mathbf{B}_{f_A}$,

$$ \mathbf{B}_{f_A} = (\mathbf{H}|_{f_A})^{-1} = (\mathbf{B}_{f_B}^{-1} + \mathbf{H}|_{f_A}^y)^{-1}. \quad (18) $$

where $\mathbf{H}|_{f_A}^y$ is the observational Hessian.

For high dimensional forcing, a direct inversion

of Eq.(18) may become untractable. We rely on the Woodbury identity formula as in Flath et al. (2011) to perform an approximation of the posterior covariance matrix as a prior update

$$ \mathbf{B}|_{f_{A,r}} = \mathbf{B}_{f_B}^{1/2} \left( \mathbf{I}_f - \mathbf{F}_r^y \frac{\mathbf{\Lambda}_r^y}{(\mathbf{\Lambda}_r^y + \mathbf{I}_f)} (\mathbf{F}_r^y)^\mathrm{T} \right) \mathbf{B}_{f_B}^{1/2} \quad (19) $$

where $\mathbf{F}_r^y$ is an orthonormal matrix storing the $r$ dominant eigenvectors $f_i^y$ and $\mathbf{\Lambda}_r^y$ is a diagonal matrix storing the associated $r$ eigenvalues $\lambda_i^y$ of the prior pre-conditioned observational Hessian $\mathbf{B}_{f_B}^{1/2}\mathbf{H}|_{f_A}^y\mathbf{B}_{f_B}^{1/2}$. A reasonable choice of the truncation $r$ is typically obtained when

$$ r: \frac{\lambda_r^y}{1 + \lambda_r^y} < 0.01 \quad (20) $$

Regarding the posterior state covariance matrix; it translates into the following

$$ \mathbf{B}_{q_{A,r}} = \mathbf{R}_A \mathbf{B}_{f_{A,r}} \mathbf{R}_A^\mathrm{T}. \quad (21) $$

where $\mathbf{R}_A$ is the resolvent matrix evaluated at the posterior state. Note that the posterior variances $\sigma_{f_{A,r}}^2$ and $\sigma_{q_{A,r}}^2$ of the mean forcing and state respectively can be directly computed by taking the diagonal of $\mathbf{B}_{f_{A,r}}$ and $\mathbf{B}_{q_{A,r}}$ so that

$$ \sigma_{f_{A,r}}^2 = \mathrm{diag}[\mathbf{B}_{f_{A,r}}], \quad \sigma_{q_{A,r}}^2 = \mathrm{diag}[\mathbf{B}_{q_{A,r}}] \quad (22) $$

## 2.5 Prior and observational covariance matrices

So far, the formalism derived is applicable for any prior and observational covariances matrices $\mathbf{B}_{f_B}$ and $\mathbf{B}_y$. We restrict now the analysis to prior statistics which neglect spatial correlations. This amounts pre-scribing a forcing prior and observational covariance matrices which are diagonal

$$ \mathbf{B}_{f_B} = \sigma_{f_B}^2 \mathbf{I}_f, \quad (23) $$

and

$$ \mathbf{B}_y = \sigma_y^2 \mathbf{I}_y \quad (24) $$

where $\sigma_{f_B}^2$ and $\sigma_y^2$ account for the variance of the prior forcing and noise.

# 3 Application to a separated flow in a converging diverging channel

## 3.1 High-fidelity data

The high-fidelity data come from a DNS of a turbulent flow in a converging diverging channel at Reynolds number 12,600 performed by Marquillie et al. (2008). To ensure a realistic representation of the observation statistics, the data used for DA and UQ are modified by injecting a synthetic Gaussian noise $\eta$ with zero mean and variance $\sigma_y^2 = 10^{-6}$. Note that this low variance induces a negligible noise in the data, as depicted in Figure 1. Besides, we can observe the small separation of the DNS around $x \in [6, 7]$. It
is important to note that the observational data, delin-
eated by the black boxes, do not cover the DNS sepa-
ration zone.

<figure><img src="image_2.png" /><figcaption>Figure 1: DNS longitudinal velocity field $u^{\text{DNS}}(\mathbf{x}) + \eta(\mathbf{x})$. The black boxes denote the location of the data used to perform the DA and UQ.</figcaption></figure>

## 3.2 *A priori* results

In this section, we present the SA *a priori* results (i.e. setting $f_u = 0$ and $f_P = 0$). In particular, we focus on its longitudinal velocity component $u_B$, which is plotted in Figure 2. We can observe that the SA model fails to predict *a priori* the correct separation. Note that the prior error, plotted in Figure 3,

<figure><img src="image_3.png" /><figcaption>Figure 2: Prior mean longitudinal velocity field $u_B(\mathbf{x})$. Black contour lines represent the DNS solution field $u^{\text{DNS}}(\mathbf{x})$.</figcaption></figure>

is mainly located downstream the hill, because of the
wrong separation length predicted by the baseline SA
model. This error can be estimated only using the
full longitudinal velocity mean observations, which
are usually not available. Instead, we can rely on the
prior state standard deviation, which do not depend on
the high-fidelity data but only on the model and forc-
ing considered. In Figure 3, we also plot the prior lon-
gitudinal standard deviation $\sigma_{u_B}$ given by $f_P$. We
observe that the distribution of $\sigma_{u_B}$ is very similar to
the one of the prior error. Besides, regions of high
standard deviation induce the most likely area to be
effectively corrected by $f_P$.

## 3.3 A posteriori results

We present here the *a posteriori* results, namely
after having observed the sparse mean observations
highlighted in Figure 1. First, the posterior means are
estimated after minimization of the cost functional as
plotted in Figure 4. Second, the posterior variances
are approximated by a low rank decomposition of

<figure><img src="image_4.png" /><figcaption>Figure 3: Top: absolute error field $\epsilon_{u_B}(\mathbf{x})$ with respect to the DNS. Bottom: prior longitudinal standard deviation field $\sigma_{u_B}(\mathbf{x})$ for the $f_P$ forcing.</figcaption></figure>

the prior-preconditioned observational Hessian, whose
spectrum is shown in Figure 5. We observe that the as-
similated state is reached in only few optimization iter-
ations. In addition, we also notice that we need to ex-
tract more eigenmodes with the $f_u$ forcing than with
$f_P$ to approximate the posterior variances. For both
forcings, the spectrum decays quickly, mainly due to
the sparse high-fidelity data.

<figure><img src="image_5.png" /><figcaption>Figure 4: Cost functional minimization.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Evolution of $\frac{\lambda_i^y}{1+\lambda_i^y}$. The dashed horizontal line cor-
responds to the limit to ensure a reasonable ap-
proximation of the posterior variances.</figcaption></figure>

The distribution of the posterior mean forcings $f_{u_A}$
and $f_{P_A}$ obtained after minimization can be observed
in Figure 6. The corrective terms are mainly dis-
tributed near the location of the observational data.
The distribution of the posterior mean longitudinal ve-
<figure><img src="image_7.png" /><figcaption>Figure 6: Top: longitudinal component of the posterior mean forcing field $f_{u_A}(\boldsymbol{x})$. Bottom: posterior mean forcing field $f_{P_A}(\boldsymbol{x})$.</figcaption></figure>

locity $u_A$ is illustrated in Figure 7. For both forcing terms, we were able to correct the erroneous separation bubble produced by the SA model. Contrary to the $f_{u_A}$ forcing, $f_{P_A}$ managed to recover partially the small separation of the DNS around $x \in [6, 7]$. This can be better visualized with the evolution of the skin-friction coefficient plotted in Figures 10-11, where the separation and attachment points can be identified. In

<figure><img src="image_8.png" /><figcaption>Figure 7: Posterior mean longitudinal velocity field $u_A(\boldsymbol{x})$. Top: for the $f_{u_A}$ forcing. Bottom: for $f_{P_A}$ forcing. Black contour lines represent the DNS field $u^{\text{DNS}}(\boldsymbol{x})$.</figcaption></figure>

Figures 8 and 9, we compare the longitudinal-velocity posterior error field $\epsilon_{u_A}$ with the estimated posterior standard deviation $\sigma_{u_A}$ for the two forcings. Overall, both types of forcing effectively reduced the error downstream of the hill caused by the incorrect separa-

tion predicted by the SA model. This is a promising result, taking into consideration that only a few data points were used during the DA process. We can also notice that the posterior error is slightly higher near the wall for $f_{u_A}$. Note that the posterior error map cannot be computed in most of turbulence applications since the full set of reference data is usually not available. Regarding the posterior standard deviation $\sigma_{u_A}$,

<figure><img src="image_9.png" /><figcaption>Figure 8: Results for the $f_{u_A}$ forcing. Top: posterior error field $\epsilon_{u_A}(\boldsymbol{x})$. Bottom: posterior standard deviation field $\sigma_{u_A}(\boldsymbol{x})$.</figcaption></figure>

it relies only on the sparse observations used and not the full set of DNS data. We observe that the $f_{P_A}$ forcing was able to estimate a high uncertainty in the DNS separation region, which is the region of highest residual error. On the other hand, the $f_{u_A}$ forcing predicted a large error far away from the wall, despite its low residual error. Importantly, the uncertainties for other reconstructed quantities of interest can also be automatically estimated. For instance, Figures 10-11 plot the posterior pressure coefficient $C_{P_A}$ and skin-friction coefficient $C_{f_A}$ for the two forcings with their standard deviations. The pressure coefficient was improved with both forcings. However, the application of $f_{u_A}$ led to a substantial degradation of the skin friction, fully eliminating the separation present in the DNS. In contrast, the $f_{P_A}$ forcing manages to recover it partially. This highlights the limitation of $f_{u_A}$ in accurately correcting near-wall regions when data is sparse. Moreover, the UQ appears more robust with the $f_{P_A}$ forcing, which successfully estimates high uncertainties in the separation region, corresponding to areas with the largest residual errors. In contrast, the $f_{u_A}$ forcing significantly underestimates the uncertainties in these critical regions.

## 4 Conclusions

This study evaluated a data assimilation and uncertainty quantification (DA-UQ) methodology using