# Autoregressive Multiplier Bootstrap for In-situ Error Estimation and Quality Monitoring of Finite Time Averages in Turbulent Flow Simulations

Christos Papagiannis, Guillaume Balarac, Pietro Marco Congedo, Olivier P Le Maître

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

## To cite this version:

Christos Papagiannis, Guillaume Balarac, Pietro Marco Congedo, Olivier P Le Maître. Autoregressive Multiplier Bootstrap for In-situ Error Estimation and Quality Monitoring of Finite Time Averages in Turbulent Flow Simulations. Computer Methods in Applied Mechanics and Engineering, 2026, 452, pp.118664. 10.1016/j.cma.2025.118664. hal-05328929

HAL Id: hal-05328929

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

Submitted on 24 Oct 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Autoregressive Multiplier Bootstrap for In-situ Error Estimation and Quality Monitoring of Finite Time Averages in Turbulent Flow Simulations

Christos Papagiannis$^{b,c}$, Guillaume Balarac$^{d,c}$, Pietro M. Congedo$^a$, Olivier P. Le Maître$^a$

$^a$CMAP, CNRS, Inria, IPP, Palaiseau, France

$^b$CNRS, Inria, Palaiseau, France

$^c$CNRS, LEGI, Grenoble, France

$^d$INP, CNRS, Grenoble, France

## Abstract

In Computational Fluid Dynamics (CFD), and particularly within Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES), the computational cost is largely dictated by the effort required to obtain statistically converged quantities such as time-averaged fields and higher-order moments. Despite the importance of accurately quantifying statistical uncertainty in unsteady simulations, no continuous and cost-effective, on-line method currently exists for monitoring the convergence quality of such statistics during runtime. This work introduces a novel, fully on-line bootstrapping approach to estimate the variance of finite-time averages without requiring the estimation of the flows Auto-Correlation Function (ACF). Unlike existing methods that rely on ACF estimation, which are often impractical due to excessive storage demands in large-scale simulations, or require off-line processing or a priori modeling assumptions, our method operates entirely during the simulation and incurs minimal overhead. The proposed technique employs a recursive update of bootstrap replicates of the time average, using correlated random weights generated via an autoregressive model. This formulation is computationally efficient: the update cost scales linearly with the number of bootstrap replicates and the dimensionality of the flow field, and the autoregressive model is inexpensive to evaluate. The method only requires storage of a small number of fields, making it suitable for large-scale CFD applications. We demonstrate the effectiveness of the approach on synthetic data from the Ornstein-Uhlenbeck process and on two canonical LES cases: a turbulent pipe flow and a round jet. We further discuss the methods applicability to simulations with non-uniform time stepping, highlighting its flexibility and robustness.

Keywords: Uncertainty quantification, Time-averaging uncertainty, In-Situ estimators, In-Situ Bootstrapping, Turbulent flows

## 1. Introduction

Turbulence remains up to date, one of the complex physical phenomena for which a universal and robust mathematical description doesn't exist and therefore research is based mainly on numerical simulations of the governing Partial Differential Equations (PDEs) whose solutions can be unsteady in time [1]. These unsteady simulations are usually validated through comparisons to experimentally acquired time averages of the flow quantities. In both simulation and experiment time-averaging is performed on the acquired

Email addresses: christos.papagiannis@inria.fr (Christos Papagiannis), guillaume.balarac@univ-grenoble-alpes.fr (Guillaume Balarac), pietro.congedo@inria.fr (Pietro M. Congedo), olivier.le-maitre@polytechnique.edu (Olivier P. Le Maître)
data, over a finite time window. Despite the fact that the governing equations are deterministic, due to
their sensitivity to the nominal conditions for repeatability, the resulting fields are, in the most realistic
applications, chaotic [2]. The chaotic nature of both measurements and simulation calls for a statistical
treatment of the resulting fields. The flow variables are treated as random variables, whose statistics
(usually the mean and standard deviation) are of the most interest to the practitioner, eventually creating
a need for applying Uncertainty Quantification (UQ) methodologies. In general there are many sources of
uncertainty in a CFD simulation: due to the choice of spatial discretization schemes, assumptions made
on boundary/initial conditions and in our case, finite time averaging [3],[4],[5]. The results of the Error
Estimation are needed in order to either validate a simulation against experimental data (that themselves
require some UQ treatment) or to assess the level of its convergence in real time. The need for quantifying
uncertainties is included in the goals of the NASA 2030 report [6] highlighting the importance of error
estimation for the broadening of the applicability of CFD to the industrial design process. On this path to
asses the error due to finite time averaging, we develop on the ideas existing in the error estimation literature
and propose a solution using an incremental bootstrapping estimator that resolves many of the pitfalls of
current approaches. The literature, on Error Estimation for time averages in CFD, has recently started to
become richer thanks to contributions from other domains of research such as time series analysis, since
modern streaming data formats require sophisticated online algorithms. Treatment of this large volume of
data could not be performed offline or aposteriori of the simulation, due to storage limitations. As observed
in an analysis of the the Top500 supercomputers, processing power is seen to grow clearly faster than data
storage capacity [7], calling for new algorithms that can provide *online* or *in-situ* as many quantities as
possible. Additionally since we are talking about simulation of physical problems, some level of problem
specificity will always be required. This adds to the challenges of designing error estimators for a wide range
of applications. Therefore an open problem, in this domain, is proposing consistent and less problem specific
approaches to quantify in an online/in-situ way the Mean Square Error (MSE) in the time averages. Aside
from a post-processing perspective, these methods would enable the practitioner to monitor the quality of the
time averaged quantities on the fly. This is important because it can be directly connected to a termination
criterion for the simulation, since a big part of the computational overhead of the simulation comes from
the need to estimate accurate statistics. Being able to monitor the error online will enable shorter and more
efficient allocation of computational time for CFD simulations, where the run-time is directly proportional
to the power consumption and therefore its carbon footprint [8]. Therefore, in the existing literature for
error estimation of time averages in CFD, we observe a lack of methods with general applicability, that
are at the same time easy to implement and cost efficient. This is because of the inherent difficulties in
estimating the error of the averages of turbulent flow quantities. A natural data structure of the results
of either simulations or experiments, in the unsteady case, is this of *time series*. Usually in time series
data, the Mean Square Error (MSE) or variance of the estimated time-average depends on the correlation
between the samples which is quantified by the Auto Correlation Function (ACF) of the time-series [9]. The
need to estimate and use the ACF is the biggest obstacle in the application of the same MSE estimator
in various physical settings and in its in-situ implementation, especially for large scale simulations. Direct
estimation of the ACF from data can become computationally prohibitive both for performance and storage
reasons. Therefore, may approaches propose fitting simpler functions to data, but at the cost of being more
problem specific. In all cases the need to implement the methods in-situ is always desirable. To answer these
problematics all at once, we propose an Error Estimator that treats the simulation as a "black box" that
generates time-series data (therefore being really easy to implement in a large and complex solver, or even
in a completely independent field of study) and has minimal storage requirements. Our estimator is based
on a multiplier bootstrap approach and for the first time, addresses the problem of Variance Estimation in
the non-constant time-step setting.

The paper is structured as follows: in section (2) we provide some definitions and notation that will be
useful throughout the rest of the paper and then state the problem with the ACF estimation. Then we expand
on the state of the art, examined under the criteria of its in-situ applicability and ease of generalization to
many physical cases. Then the proposed approach is described in section (3) where also a general description
of bootstrapping and how it can be applied to an online setting that takes into account that the data being
resampled have a temporal autocorrelation, is given. For this a variant of the method proposed in [10] is
used, that is modified to take into account unevenly spaced time-series. We then demonstrate the accuracy of the method through 3 cases in sections (1), and (5). We use the Ornstein Uhnlenbeck process, in which we can validate the proposed estimator since we can derive an analytical expression of the desired Mean Square Error and two other LES cases, a round pipe and a round jet flow that reflect more realistic applications. In all cases, we compare the estimates using the proposed method, to the state of the art that we consider to be the direct estimation of the integral autocorrelation time through the ACF using the existing samples.

## 2. Notation, problem statement and the state of the art.

Below we will present some basic notation that will be used throughout the paper. We then underline the weak points of the existing literature, under the scope of real time monitoring and error estimation in a CFD setting. Finally, we examine how the state of the art approaches the problem: its ability to be implemented online and with a "black box" solver logic, where an implementation would not be very intrusive for a complex computational code.

### 2.1. Exact definition of the variance of the time average.

The results of an unsteady flow simulation are collected (sampled point-wise on the mesh) in the form of time series {$u_n$}$_{n=1}^N$ that can be seen as the collection of realizations of a process $u_n$ or as the discrete time history of the continuous time evolution of the variable, $u(t)$. The process $u_n$ that we have observed, has expectation $\mathbb{E}[u_n]$ that can be estimated using a sample mean estimator. Under the ergodic hypothesis, the time average, noted as $\mu$, can be computed by the expectation of the process $u_n$. Due to ergodicity, the finite time average can be calculated by the expectation,

$$ \mu_T = \frac{1}{T - t_0} \int_{t_0}^{T} u(t) \, dt = \mathbb{E}[u_N]. \quad (1) $$

The time average is defined as the infinite limit of the finite time average $\mu_T$ of the process, over a time window $T$,

$$ \mu = \lim_{T \to \infty} \mu_T = \lim_{T \to \infty} \frac{1}{T - t_0} \int_{t_0}^{T} u(t) \, dt. \quad (2) $$

The time average $\mu$ is estimated with the Sample Mean Estimator (SME) for the expectation $\mathbb{E}[u_N]$ using a finite number of samples $N$ and a constant time step:

$$ \hat{\mu} := \frac{1}{N} \sum_{i}^{N} u_i \sim \mathcal{N}(\mu, S^2), \quad (3) $$

where the estimate $\hat{\mu}$ is also random (like $u_n$) and normally distributed with a mean $\mu$ and a variance $S^2$ [9]. The finite time averaging introduces the Mean Square Error $S^2$ in the sum that should reduce to zero for $N \to \infty$. Since the estimate $\hat{\mu}$ at each point in time, is also random, in order to quantify its error we seek to quantify its variance $S^2$ as $\hat{\mu} \to \mu$. The variance is expected by the Central Limit Theorem to decrease proportionally to $N^{-1}$ with $N$ being the sample size. However the samples in {$u_n$}$_n$ are correlated and the variance of the estimator of equation (3) will be proportional to the variance $\mathbb{V}[u_n]$ of the process, decreasing to zero as we acquire more samples ($\mathbb{V}[u_n] \stackrel{N}{\to} 0$), but with a correction factor $C_u$ to account for the auto-correlations in $u$,

$$ S^2 = \frac{1}{N} \mathbb{V}[u] C_u. \quad (4) $$

According to simple calculations the variance can be shown to be a function of the autocovariance $\gamma_k$ of the values of the time series and its variance $\mathbb{V}[u_n]$ [9],
$$S^2 = \mathbb{V}[\hat{\mu}] = \mathbb{V}\left[\frac{1}{N}\sum_{i}^{N} u_i\right] = \frac{1}{N}\left[\gamma_0 + 2\sum_{k=1}^{M}\left(1 - \frac{k}{N}\right)\gamma_k\right] = \frac{1}{N}\gamma_0\underbrace{\left[1 + 2\sum_{k=1}^{M}\left(1 - \frac{k}{N}\right)\rho_k\right]}_{C_u} \quad (5)$$

Above, $\gamma_k$ is the autocovariance of the process (calculated at discrete lags $k = 1, 2, \dots, M$):

$$\gamma_k = \mathbb{E}[u_n u_{n-k}] - \mathbb{E}[u_n]^2, \quad (6)$$

The autocovariance at lag $k = 0$ is the variance of the process:

$$\gamma_0 = \mathbb{E}[u_n^2] - \mathbb{E}[u_n]^2 = \mathbb{V}[u_n], \quad (7)$$

where because of the assumption of stationarity, $\mathbb{E}[u_n] = \mathbb{E}[u_{n-k}] = \mu$. The normalization by $\gamma_0$ in equation (8) gives the autocorrelation function (ACF) $\rho_k = \gamma_k/\gamma_0$. The autocovariance is estimated as

$$\hat{\gamma}_k = \frac{1}{(N-k)} \left( \sum_{i=1}^{N-k} u_i u_{i-k} - \hat{\mu}^2 \right). \quad (8)$$

In equation (8) we require the sum of the ACF over the lags $k = \{1, \dots, M\}$. This sum approximates the decorrelation time or the integral time scale of the flow in the fluid mechanics literature [11]. This is the penalty factor $C_U$ that affects the convergence of the estimator of (4). The lag time is indexed by $k$ that takes integer values in $[0, M]$, $M \le N$ and plays an important role to the estimation of the integral autocorrelation time $\tau_I = \int^\infty \rho(\tau)d\tau$. In the entire paper we will estimate the variance $S^2$ using the estimator $\hat{\sigma}_N^2$,

$$\hat{\sigma}_N^2 = \frac{1}{N} \left[ \hat{\gamma}_0 + 2 \sum_{k=1}^{M\ll N} \left( \frac{1}{N-k} \left( \sum_{i=1}^{N-k} u_i u_{i-k} - \hat{\mu}^2 \right) \right) \right] \quad (9)$$

on a sample of size $N$ and use this estimate as the reference estimate to compare with our estimator.

## 2.2. The issue with the ACF estimation and the state of the art.

We point to three major issues with trying to estimate the ACF to acquire the error of the SME for the entire field in an in-situ fashion. The first problem is common to the entire literature, where we have not seen it being addressed often and the other two are the main problems that most of the methods in the literature address. The first problem is that no methods in estimating the variability of the SME of equation 3, take into account the fact that the time step of the simulation and therefore the spacing of the time series can be variable. The estimator of equation (9) is constructed for equidistant time-series where time averages are equivalent to ensemble means,

$$\int_{t_0}^{T} u(t) dt \approx \frac{1}{N \Delta t} \sum_{i}^{N} u_i \Delta t = \frac{1}{N} \sum_{i}^{N} u_i, \quad (10)$$

which is not always the case. If the time step is not fixed, we can not ignore its contribution in the time average:

$$\hat{\mu} := \frac{1}{\sum_{i}^{N} \Delta t_i} \sum_{i}^{N} u_i \Delta t_i. \quad (11)$$

The $\Delta t$ is usually defined by the CFL number $CFL = \frac{u\Delta t}{\Delta x}$ that connects the size of the temporal discretization with the spatial one, according to the local velocity, in order to guarantee stability in codes with explicit time integration. Therefore, if the CFL is employed, the time step is also random and correlated with the velocity $u(x, t)$. This point is not frequently addressed in existing literature. It is a common
assumption that the $\Delta t$ is kept constant, as in our presented examples. More on the importance of this variability can be found in the appendix (Appendix B). In this paper we propose an alternative formulation that will be used to circumvent this issue. For the calculation of the variance $\hat{\sigma}_N^2$, we require to estimate the ACF to be able to estimate the integral autocorrelation time $\tau_I$. This can become problematic really fast. By estimating it from the existing data with equation (8), we run into the two remaining problems: how much data should we store to efficiently estimate the ACF (can we even afford to do it for every point in the field) until the lag time ($M$) that we need for the integration and what is that lag $M$. Arbitrary answers to these questions can lead to unreliable estimates of the $\tau_I$ since as we can see from figure (1) the finite length of the data used to estimate the ACF induce oscillations at its tail [11]. Therefore it is important to answer the question of the necessary length $M$ that we require for the integration.

<figure><img src="image_3.png" /><figcaption>Figure 1: Comparison of the estimated ACF of an OU process, calculated using a limited number of data vs the case of using the entire data set (M=N): (a) Truncated at $M = 100$ and (b) truncated at $M = 1000$. The shaded red area corresponds to the Correlation Coefficient $C_u$ of equation (4) that comes from the integral of $\gamma_k$</figcaption></figure>

$M$ plays an important role in the accurate estimation of the integral of $\rho_k$ and then $S^2$. In theory $M$ should be as long as possible, but in practice $M$ can be at most $M = N$. Therefore selecting the appropriate lag length would require an a posteriori analysis of the ACF which makes methods that don't automate the lag selection, inapplicable for a truly online setting. Unfortunately there is no straightforward way for the lag length selection. As we can observe from figure (1), if we increase the lags we get a different estimate of the characteristic timescale $\tau_I$ that would in turn affect the correlation coefficient $C_u$. The most general approach is to truncate the ACF in its first zero crossing [11], but the lags required for this are also unknown a priori of the simulation. Another problematic that has already been addressed in [12] is how to efficiently treat online the storage requirements for the products $u_n u_{n-k}$ in the estimation of $\hat{\gamma}_k$,

$$ \hat{\gamma}_k = \frac{1}{N-k} \left( \sum_{i=1}^{N-k} u_i u_{i-k} - \hat{\mu}^2 \right), \qquad (12) $$

since without a special treatment one would need to store at least $M$ in length of the history of $u$ for the entire field. This need for large storage is what makes the estimation of the ACF generally not feasible for the entire field. To avoid both oscillations and storage overheads, $M$ needs to be sufficiently small but large enough for the first zero crossing, which changes according to the physics of the process at each point in the flow. To circumvent the above problematics, methods proposed in the literature follow three general directions: of resampling the data to acquire non correlated time series, of block averaging methods and fitting some model for the ACF to existing data.

By resampling an existing time series one needs to estimate the resampling frequency which itself is connected to estimating a decorrelation time of the flow [13]. This results in storing a signal without autocorrelation to treat it with the conventional variance estimators. That would be through sub-sampling the flow field data at a frequency smaller than the decorrelation frequency $f_d = \frac{1}{\tau_I}$ so the resulting time series is identically and independently distributed (i.i.d.) with $\tau_I$ being the integral time scale [14]. One could see that as a reduction in the number of samples in the estimator of equation (4) where with the
new *effective* number of samples have $C_U = 1$. A more in depth review of the effective samples method for error estimation can be found in the work of Smith et al [15].

In the case of a Block Averageing approach [16], [12] one effectively "filters" the time series at some (constant or not) length. This family of methods is widely used to estimate uncertainties in Monte Carlo type estimators in various fields [17], [18] and inspired the relevant estimator proposed in the work of Grenouilloux [19]. These methods also called Batch-Means approaches can be combined also with the ACF modelling approaches as seen in the work of [16]. The main issue with these approaches is the block length that has been found to influence heavily the resulting variance estimation [20][21]. A data efficient implementation for the block-averaged methods was proposed by Russo et al. [12] and [16] but the problem of the selection of the block length persists. Block length for averaging and estimating uncertainties for correlated data sets has been explored independently, especially for resampling estimators [21]. A more extensive study on the block length of the Batch-Means approaches can be found in [21] and [16]. Using the advances in the block averaging above, the work of [22] et al. proposed to model the ACF as a weighted sum of exponentials with a set of hyperparameters that need to be learned and used the efficient in-situ update techniques for the batch averaging to enable the in-situ parameter estimation of the model. That approach was found to perform better than older models of the ACF (in challenging regions like near walls where the ACF would differ drastically compared to more turbulent regions of the flow [22]). In the same category of methods we include the broad family of ARMA (Auto Regressive Moving Average) models that are used to estimate the ACF or used to model the recorded time series as seen for the first time in CFD in [4]. However, the limitations for functional models for the ACF still persist in the existing literature. The ACF changes significantly in shape, from point to point in the mesh. Periodicity or changes in the intensity of turbulence can change the shape of the ACF, therefore rendering the expensive fitting approaches more necessary (increasing the cost in computational time and storage, if possible at all) and prevent a single model approach from being valid everywhere in the domain. In our approach we focused on addressing the problem of the case specificity that all the existing methods have as well as the constant time-step assumption. As an alternative approach to these parametric estimators for the MSE, we propose to use bootstrapping. The idea of using the bootstrap for error estimation of time averages is not new [23], but is still at its infancy [15]. However *classical* bootstrapping has notoriously high computational costs and requires the entire dataset, that has been sampled, to be stored, which in a time series setting would also make it infeasible for an online application. In this paper we propose a case independent approach based on the Auto-Regressive Multiplier Bootstrap of Palm et. al [10]. With the proposed methodology, we avoid the problems of having to estimate the *spatially dependent* [24] correlation length of the processes and the ACF, of having to change model for the ACF (each region with different physics could need a fundamentally different ACF model [25]). With this approach, given that we allocate enough of our storage resources for having enough replicas for the bootstrapping, we can, for the first time, have a truly online estimate of the MSE of a time average in the entire domain with realistic computational and storage overheads.

### 3. Proposed approach: The Auto-Regressive Multiplier Bootstrap

In the majority of the literature resampling methods and especially bootstrapping, are being used to estimate a reference statistic (like error bars and confidence intervals) due to their high level of accuracy and lack of assumptions about the underlying distribution [21]. Bootstrapping has proven to be a powerful tool for validation of other type of estimators and is usually used as a reference, when the true distribution is unknown [26]. We will start by providing an intuitive explanation of the bootstrap as its formalization is involved and is out of the scope of this paper. We focus on the use of bootstrapping for estimating the variance $V[\hat{\psi}]$ of a SME. This SME $\hat{\psi}$ is applied to a dataset $U$. To estimate the variance of the SME, we can use the sample variance estimator for a set of $B$ replicas $U^*$ of an initially available data set $U$ [27]. Schematically we can see the bootstrap as a process that "replicates" the initial data as illustrated in figure (2) and then uses the statistics of the replicas to infer the statistics of the underlying distribution.

One can see the intuition behind bootstrapping methods, as the ability to sample the distribution of the original data set, by using only the available data, which would seem like one is "trying to fly pulling on their shoelaces", hence the name, as proposed by Efron [28]. The parameters influencing the quality of
<figure><img src="image_4.png" /><figcaption>Figure 2: Schematic illustration of the idea behind the bootstrap: The initial dataset U is replicated in B replica sets $U_b^*$ and then the statistic of interest $\hat{\psi}(\cdot)$ is calculated on the replicated sets $U_b^*$.</figcaption></figure>

the bootstrapped estimates are the number of replica sets that are being created from the initial data, the number of available samples N in U and the dependencies in the data (e.g. i.i.d. data vs correlated data).

In our work, we adapt the estimator proposed in the work of Palm et al. [10] to handle non uniform time series data in the context of CFD simulations, as a cost effective and accurate tool to estimate the error of the SME of a turbulent flow variable for the entire solution at each simulation step. So far the estimators proposed for the MSE of the SME in CFD relied on the assumption that the data were conveniently sampled equidistantly in time. In appendix (Appendix B) we provide some data that suggest that this is not the case where one employs the CFL condition to define the time-step. To overcome this we apply the Autoregressive bootstrap estimator of [10] to a transformed variable v that contains the variability of the time step. Suppose we want to estimate $\mathbb{V}[\hat{\psi}]$ using bootstrapping. In our case the functional $\hat{\psi}$ is defined to be the SME:

$$ \hat{\psi}(U) = \hat{\psi} := \frac{1}{N} \sum_{i}^{N} u_i. \qquad (13) $$

However in the case of a non-constant time-step, the functional of equation (13) does not estimate the time average $\mu$. For the time average one wants to estimate:

$$ \mu_T = \frac{1}{T} \int_{0}^{T} u(t) \, dt \approx \frac{1}{\sum_{j} \Delta t_j} \sum_{i}^{N} u_i \Delta t_i \qquad (14) $$

Which is not a SME. The core idea of the paper consists of using the Autoregressive Multiplier Bootstrap on a transformed version of equation (14) that is a SME:

$$ \hat{\mu} = \frac{1}{\sum_{j} \Delta t_j} \sum_{i}^{N} u_i \Delta t_i = \frac{1}{N} \sum_{i}^{N} \underbrace{\frac{N \Delta t_i}{\sum_{j} \Delta t_j}}_{z_i} u_i \Delta t_i = \frac{1}{N} \sum_{i}^{N} v_i, \qquad (15) $$

essentially, applying the multiplier bootstrap to the product $v = zu$, where $z_i$ are weights that tend to 1 as the simulation advances. Keeping this transformation in mind, we will describe how we will use the multiplier bootstrap to estimate the variance of the time average $\mu_T$ that now is described by the ensemble average of v (equation (14)). In the following part of the paper we motivate the reader that an asymptotically equivalent online implementation of a bootstrapping estimator has various advantages in CFD (where a large volume of data is produced and needs to be processed without storage) compared to the implemented state of the art estimators that depend on explicitly estimating the ACF of the process. The ACF estimation requires assumptions about its form [22] and has large and complex to determine storage requirements, in order to be implemented efficiently online. In the following sections we will explain bootstrapping for i.i.d. realizations, extended to correlated time series and then naturally explain the construction of the (asymptotically equivalent) online version of the bootstrapping.
### 3.1. The bootstrapping for i.i.d. data

We start by describing the bootstrap for a set with i.i.d. samples. Lets assume we have $N$ measurements of some flow variable $u$, in a set $U = \{u_1, u_2, \dots, u_N\}$. Lets also assume that these data are i.i.d. realisations of an ergodic and stationary generating process. We can calculate the ensemble or empirical average of the elements in set $U$ as:

$$ \hat{\psi}(U) = \hat{\psi} = \frac{1}{N} \sum_{n}^{N} u_n. \qquad (16) $$

We wish to estimate the variance $V[\hat{\psi}(U)]$ of the empirical average $\hat{\psi}$ of $U$. For this reason we create $B$ different replica sets $\{U_b^*\}_{b \le B}$ as illustrated in figure 2. Each replica set $U_b^*$ contains samples drawn with replacement from $U$ meaning that the same element can be picked again to be in each replica set. Therefore each element has probability $1/N$ of being picked at each time we resample. This will become important later on, but here the necessity of *predetermining* $N$ becomes evident. As an example,

$$ \begin{aligned} U_1^* &= \{u_2, u_3, \dots, u_1, \dots, u_n\} \\ \vdots \\ U_B^* &= \{u_1, u_4, \dots, u_n, \dots, u_6\} \end{aligned} $$

Since $B$ replica sets are available, we can have $B$ different samples of the empirical mean of $U$:

$$ \psi_b^* = \frac{1}{N} \sum_{n=1}^{N} u_n^{(b)}, \quad 1 \le b \le B \quad u_n^{(b)} \in U_b^* \qquad (17) $$

Now if these $B$ replicas have a mean (over $B$),

$$ \tilde{\psi} = \frac{1}{B} \sum_{b=1}^{B} \hat{\psi}_b^*, \qquad (18) $$

we can estimate the variance $V[\psi^*]$ of the bootstrapped statistic as:

$$ \hat{\sigma}^2 = \frac{1}{B-1} \sum_{b=1}^{B} (\psi_b^* - \tilde{\psi})^2, \qquad (19) $$

that would provide asymptotically a consistent estimator for the variance of the sample mean. In other words if we were to use the centering on $\hat{\psi}$ we would quantify the variability of the bootstrapped mean that doesn't include the bias due to its own error from the sample mean. Therefore, the centering *includes the bias* due to the distance of the two means $(\hat{\psi} - \tilde{\psi})^2$. Asymptotically, as the bootstrapped mean converges to the sample mean, the variance of the bootstraps $\hat{\sigma}^2$ converges to the variance of the sample mean,

$$ \hat{\sigma}^2 \xrightarrow{B \to \infty} V[\hat{\psi}], \qquad (20) $$

and as $N$ increases the sampling error $V\hat{\psi}$ converges to 0. However since we are interested in the time evolution of the variance, including the bias due to $|\hat{\psi} - \tilde{\psi}|$ is important. Hence we choose to deviate from the standard literature approach and adapt the proposed approach of the "cheap bootstrap" that was proposed in [27] where consistency is proved even for a very small number of replicas $B$. The difference is that we centered the variance on $\hat{\psi}$ and not on $\tilde{\psi}$ including in our estimate the error between $\hat{\psi}, \tilde{\psi}$, since we are interested in estimating the MSE of the estimator $\hat{\psi}$. So in the end the estimator that we will use is:

$$ \hat{\sigma}^2 = \frac{1}{B-1} \sum_{b=1}^{B} (\psi_b^* - \tilde{\psi})^2. \qquad (21) $$
## 3.2. The Bootstrapping for time-series with correlation

Irregardless of what estimator we use fot the variance, the most important characteristic of the method is how we perform the resampling in a way that accounts for the autocorrelation in the data. In most applications in turbulent flows the resulting time series {$u_n$}$_{n \le N}$ are auto-correlated up to some integer index k. In this case if we use the previous notation, $U = \{u_1, u_2, ..., u_N\}$ is a set in which the ordering of the elements matters and therefore we have to perform resampling carefully. Evidently, resampling with replacement from $U$, will create replicas {$U_b^*$}$_{b \le B}$ that don't reflect the temporal correlation of the elements of $U$ and therefore the bootstrapped variance will be underestimated. To overcome this we need to treat the now ordered elements in $U$ differently. Assuming that the correlation length is not larger than the existing time history, we can partition the set $U$ into subsets $U_i^{sub}$ of $U$ that each contains $s_B$ consecutive $u_i$,

$$ U = \{ \underbrace{\{u_1 u_2 u_3\}}_{s_B} \dots \underbrace{\{u_4 u_5 u_6\}}_{s_B} \dots \underbrace{\{u_7 u_8 u_9\}}_{s_B} \} \quad (22) $$

where $s_B$ is larger than $\tau_I = \sum_k^M \rho_k$ that depends on $M$. Then we could treat these "blocks" of the time series as independent from each other. If we choose the block length to be constant and equal to $s_B$ then we get $n_B$ blocks so that $N = n_B \cdot s_B$. The sampling with replacement now applies to the entire block of $s_B$ elements. Now the replicas are themselves time-series that are created using blocks of the original one:

<figure><img src="image_5.png" /><figcaption>Figure 3: Replica time-series by sampling with replacement from the partitioned set $U$: Replicas are now time-series that retain the correlation structure but in the scale of the blocks.</figcaption></figure>

There are many options for block selection with the most notable ones being non-overlapping blocks of constant or variable size, overlapping blocks (but then one needs to investigate the optimal overlap) [29] and then random sized blocks [21],[16]. Now the variance $\hat{\sigma}^2$ is acquired in the same way as in equation (19) with the only difference being the way the replicas are created. Notably, there is a trade-off between the block size $s_B$ and the estimated variance. If the block length is very large, the replica time series will exhibit high similarities with the generating one since most of the data remain unperturbed. Thus, $\hat{\sigma}^2$ will be an underprediction of $V[\hat{\psi}]$ since replicas won't differ too much from each other. If $s_B$ is too small, the replicas will give similar variance to the classical i.i.d. bootstrap which would also be an underprediction of the true variance $V[\hat{\psi}]$ due to the lack of the correlation (the ordering of $u_i$ is in the time-series was perturbed destroying correlations). The optimal block length is the one that both captures the correlation structure (being in the scale of the correlation length of the process) but also allows for enough samples to measure the variance during the replication [21], this means having different enough replica time-series. More studies on optimal block selection for Block Bootstrap or in general Batch methods, have been extensively performed in the literature but still require a very problem-specific approach [16], rendering Batch based methods more difficult to integrate in a broad spectrum of CFD applications where the physics of the flow change drastically the autocorrelation structure from point to point.
### 3.3. Bootstrapping performed online?

So far we viewed the bootstrap as an offline procedure that needs to access all the data of the simulation, meaning that we require them to be stored in memory, something undesirable. Since we want to promote the use of bootstrapping estimators, we deem it is of interest to explain in more depth the construction of the online version of such estimators. To begin we adopt now the notation of [30] where the sampling is encoded as a vector of coefficients/weights. Now we will view the averaging as an inner product between the vectorized data $U$ and a weight vector $w$:

$$ \hat{\psi} = \frac{1}{N} U \cdot w^T = \frac{1}{N} [u_1 \ u_2 \ \dots \ u_n] \cdot [1 \ 1 \ \dots \ 1]^T. \qquad (23) $$

For the empirical average $\hat{\psi}$, $w \in \mathbb{R}^N$ is a vector of ones as seen in equation (23). We can now express the sampling procedure as drawing different vectors $w$ to avoid storing $U$ since the sum of the inner product can be calculated incrementally but weighted with the known elements of $w$:

$$ \hat{\psi}_n = \frac{n-1}{n} \hat{\psi}_{n-1} + \frac{w_n}{n} u_n. \qquad (24) $$

Therefore we don't need to store the entire time-series at each point and then resample it, we just need to store a sampling matrix $W \in \mathbb{N}^{B \times N}$ that corresponds to $B$ replicas (rows) and has $N$ elements (columns) that sum to at least $N$:

<figure><img src="image_6.png" /><figcaption>Figure 4: Schematic illustration of the process of calculating the directly the averages of the replicas without storing them but online since we already know the weights $W$.</figcaption></figure>

This would cost $O(B \times N)$ and would be prohibitive for a large simulation. It is evident that to generate the sampling vectors $w_b$ we need to know $N$ a priori, otherwise we run into two problems: 1) we don't know the probabilistic rule with which to generate the weights and 2) if we fix $N$ to solve 1) we need to wait until $N$ steps have been performed, to report an unbiased average that corresponds to $N$ steps. Thus we could implement the bootstrapping in a "quasi online" fashion where we predetermine a number of iterations $N$ up to which we collect data and then perform the bootstrap. This is not ideal for our scenario where we seek a truly online estimator of the variance of the time average to use as a monitoring tool of the simulation. The solution to this problem comes from writing the update of each one of the $B$ replica averages $\psi_b^*$ as an autoregression:

$$ \psi_b^{*(n)} = \frac{n-1}{n} \psi_b^{*(n-1)} + \frac{1}{n} w_b^{(n)} u_n, \qquad (25) $$

where the new data are weighted by a weight $w_b^{(n)}$ exactly like in the case of equation (23) but now the weights are not predefined. The fact that the weight depends on the iteration $n$ is what gives the flexibility to design specific weighting rules to have an asymptotically consistent estimator at each step. The weights are drawn from some regression that has, vaguely speaking, desirable properties. Now that we described in general how to approach the issue of implementing the bootstrapping in-situ, we address the second important question: How do we select the weights to account for a correlated time-series? The answer
to this was given in the method of Palm et al. [10] where the suggested auto regressive weights sequence provides an asymptotically consistent bootstrapping estimator. Avoiding the technical details and definitions we focus on the implementation of the method. Assuming that our data are generated from a stochastic process $(u_t)_{n \le N}$ and satisfies certain conditions regarding its expectation and correlation structure [10], we can generate the weights $w^{(b)}$, for each different replica set, from the following Auto Regression:

$$w_b^{(n)} = 1 + \rho_n(w_b^{(n-1)} - 1) + \sqrt{1 - \rho_n^2} \zeta_b, \quad \rho_n = 1 - n^{-\beta}, \quad \zeta^{(b)} \sim \mathcal{N}(0, 1) \quad (26)$$

In the above, the weight $\rho$ that is a function of the iteration counter $n$ and the parameter $\beta$, controls the strength of dependence of the Auto Regressive (AR) weights. As $n$ grows larger the weighting tends to the one of the Multiplier bootstrap of Van der Vaart and Wellner [31]. The principle idea behind the method is that since we assume a decaying correlation between consecutive samples, the weights should exhibit similar behavior to prioritize more correlated samples of the initial iterations ($n <<$). The parameter $\beta \in (0, \frac{1}{2})$ is important as it controls the convergence rate of the estimator. It can be proven that $\beta$ has an optimal value at $\beta = \sqrt{2} - 1$, for which the bootstrapping has a balanced convergence rate [10] between its bias and variance. The updates that need to be performed are described by algorithm (1):

**Algorithm 1** The updating Auto-Regressive algorithm of Palm et al. [10]

1: **Initialize**: $\psi_b^* = 0, w_b^{(0)} = 0, \bar{w}_b^{(0)} = 0, b = 1, \dots, B$.

2: **For**: times $n = 1, 2, \dots, N$

3: Calculate $\rho = 1 - n^{-\beta}$.

4: Observe new data $u_n$.

5: **For**: $b = 1, \dots, B$

6: Simulate $\zeta_b \sim \mathcal{N}(0, 1)$.

7: Update (in this order):

$$w_b^{(n)} \leftarrow 1 + \rho_t (w_b^{(n-1)} - 1) + \sqrt{1 - \rho_n^2} \zeta_b,$$

$$\psi_b^{*(n)} \leftarrow \frac{(n-1) \bar{w}_b^{(n)} \psi_b^{*(n-1)} + u_n w_b^{(n-1)}}{(n-1) \bar{w}_b^{(n)} + w_b^{(n)}},$$

$$\bar{w}_b^{(n)} \leftarrow \frac{(n-1) \bar{w}_b^{(n-1)} + w_b^{(n)}}{n}.$$

8: Compute empirical variance of $\{\hat{\psi}_b^{*(n)}, \dots, \hat{\psi}_B^{*(n)}\}$ to quantify MSE of $\hat{\psi}^{(n)}$ using:

$$\hat{\sigma}_{BAR}^2 = \frac{1}{B-1} \sum_{b}^{B} (\psi_b^{*(n)} - \hat{\psi}^{(n)})^2 \quad (27)$$

In equation (27) we use the estimator for the variance that was proposed in [27] (by centering the variance on the sample estimate and not the bootstrapped estimate) where using a small number of replicas $B$ is justified theoretically. Our aim is to propose a method that can provide the MSE estimates in the entire field with minimal computational overreads and minimal model assumptions. We opt to compare our method to the state of the art (SotA), on the basis of accuracy, keeping in mind the method's ease of implementation. Our main motive for proposing the use of this method is implementation ease and the ability of the method to perform equally well in different physical scenarios without assumptions or modelling. We remind the reader that the Auto-Regressive Bootstrap is dependent only on one user defined parameter: the number of replicas. On the contrary, in any method that models the ACF two main questions arise: what form of the ACF should we assume and how many lags should we select according to the point in the domain where we want to estimate the ACF (physics dependence). Evidently, previous methods revolve around finding
workarounds to the previous problematics. The proposed method does not suffer from these problems.

## 4. Verification on the Ornstein Uhlenbeck process.

We begin the demonstration of the BAR method using the example of the Ornstein Uhlenbeck (OU) process. This model was selected because we can derive analytically the variance of the process and the variance $V[E[u_t=T]]$ of the of the expectation $E[u_t=T]$. We can then use the analytical results to compare them to the two estimation methods under examination: the ACF based approach and the BAR method.

### 4.1. Analytical results of the time evolution of the variance of the OU process

The OU process is the solution of the SDE:

$$du_t = -\theta(u_t - \mu)dt + \sigma dW_t \quad (28)$$

where $\mu \in \mathbb{R}$ is an arbitrary selected asymptotic mean. The parameter $\theta \in \mathbb{R}$ controls the speed of the reversion of the process to the asymptotic mean $\mu$ and $\sigma \in \mathbb{R}$ controls the intensity of the fluctuations of the process. By $W_t$ we denote the Weiner process driving the stochastic part of the dynamics. From equation (28) we can derive the conditional (on an initial condition $u_0 = 0$) expectation of the process $u_t$:

$$\mathbb{E}[u_t] = \mu_T = \mu(1 - e^{-\theta t}) \quad (29)$$

The covariance of the process as a function of lag time $\tau$ is:

$$\gamma(\tau) = \mathrm{Cov}[u_t, u_{\tau}] = \frac{\sigma^2}{2\theta} \left( e^{-\theta|t-\tau|} - e^{-\theta(t+\tau)} \right) \quad (30)$$

The reason we select this model as a reference, is that we can derive the analytical expression for the time evolution of the variance $\sigma_{\mu_T}^2$ of the time average $\mu_T$. The variance of the integral $\mu_T$ is given by,

$$\sigma_{\mu_T}^2 = V[\mu_T] = \mathbb{E}[(\mu_T)^2] - \mathbb{E}[\mu_T]^2 = \frac{1}{T^2} \int_0^T \int_0^T \mathrm{Cov}[u_t, u_{\tau}] dtd\tau \quad (31)$$

$$= \frac{\sigma^2}{2\theta^2 T} \left( 2 + \frac{e^{-2\theta T} - 1}{\theta T} \right). \quad (32)$$

### 4.2. Estimation of the RMSE using the ACF and the BAR methods

To demonstrate the effectiveness of our method we compare the estimation techniques involving the ACF with the Auto-Regressive Bootstrap on three test cases. The time evolution of the OU process and two more complex LES cases. The comparison will be structured as follows. We consider that all existing methods rely on estimating the integral time scale from a model or direct estimate of the ACF:

$$\hat{\sigma}_{ACF}^2 = \frac{1}{N} \left[ \gamma_0 + 2 \sum_{k=1}^{M} \left( 1 - \frac{k}{N} \right) \gamma_k \right] \stackrel{M \ll N}{=} \frac{1}{N} \left[ \gamma_0 + 2 \sum_{k=1}^{M} \gamma_k \right], \quad \hat{\gamma}_k = \frac{1}{N-k} \sum_{i=1}^{N-k} U_i U_{i-k} - \hat{\mu}^2. \quad (33)$$

On this estimator we take $M \ll N$ so that we don't include the large lag oscillations in the integral time-scale estimation, which would lead to an overestimation of the variance. The selection of $M$ will be done after visually inspecting each ACF estimate using $M \gg$ and then deciding the amount of necessary lags. In the literature, usually the ACF is integrated until the first zero crossing ($\gamma(\tau_0) = 0$) and then $M = \tau_0/\Delta t$ [11]. This is a rather conservative choice as we would prefer to overestimate $\tau_i$, that would give positive bias to the MSE rather than underestimate the error. We name the Root Mean Square Error of the mean (RMSE) using the ACF based approach $\hat{\sigma}_{ACF}^2$. We compare this estimator with the Bootstrap estimator with Auto-Regressive weights: the *Bootstrap with Auto Regressive weights* BAR method ($\hat{\sigma}_{BAR}^2$).
In the present paper, we compare the time evolution of the RMSE based on the ACF estimator (33) with the estimation performed with the Bootstrapping with Auto Regressive weights (BAR method). Then the two estimation methods will be compared to the analytically derived RMSE $\sqrt{\sigma_{\mu_T}^2}$. For simplicity of notation, we will write the analytical RMSE as $\sigma_{ref}$ and the estimated RMSE as $\hat{\sigma}_{BAR}, \hat{\sigma}_{ACF}$. The parameter $\beta$ of the method remains fixed at $\beta = \sqrt{2} - 1$ following [10]. In Appendix A we will explore the selection process of $\beta$ in more depth. The OU is integrated using $\theta = 1$ and $\sigma = 2$ where the ACF is integrated until $M = 100$ lags. For a given realisation of the time series $\{u_n\}_{n \le N}$ the estimator based on the ACF $\hat{\sigma}_{ACF}$ is dependent only on the number of lags $M$ selected. Repeating the calculation with the same $M$ would give the same time evolution of RMSE. However, the BAR method has a sampling error due to only $B$ replicas being and by depending on the random variable $\zeta_b$ each time evolution is unique (as seen in algorithm (1)) and as one can observe from figure (5). By repeating the estimation for a fixed $B$, we end up quantifying its sampling

<figure><img src="image_7.png" /><figcaption>Figure 5: Estimation of $\hat{\sigma}_{BAR}$ using the BAR method, repeated 2 times: each time the sequence of $\{\zeta_n\}$ used for the auto regression, are different resulting in different time evolutions but using the same number of replicas $B$ maintaining the average behavior.</figcaption></figure>

error due to finite $B$. Therefore, in order to draw more robust conclusions, we report the mean behavior of the estimation of the BAR method for $B$ replicas, using the sample average,

$$ \mu_{\hat{\sigma}_{BAR}} = \frac{1}{R} \sum_{r=1}^{R} \hat{\sigma}_{BAR}^{(r)} \quad (34) $$

and Standard Deviation from repeating the estimation for $R = 100$ times,

$$ \sigma_{\hat{\sigma}_{BAR}} = \sqrt{\frac{1}{R-1} \sum_{r=1}^{R} \left( \left( \hat{\sigma}_{BAR}^{(r)} \right)^2 - \hat{\mu}(\hat{\sigma}_{BAR}^2) \right)^2} \quad (35) $$

### 4.3. Experimental Results

The time series data produced by a simulation of the OU process will be used to verify the validity of the estimator, as the time evolution of the RMSE of the SME is derived analytically. The time average is acquired over a time window of $T = 8000[tu]$ (in time units) or equivalently $N = 8000/0.05 = 160000$ steps with a $\Delta t = 0.05$. The time evolution of the different estimators for the RMSE can be seen in figure (6). The BAR estimator uses $B = 20$ replicas. The error bars are constructed using $\sigma_{\hat{\sigma}_{BAR}}$ from equation (35).

$$ \epsilon_{ref} = \frac{|\hat{\sigma}_{BAR} - \hat{\sigma}_{ref}|}{\hat{\sigma}_{ref}} \quad (36) $$

For aiding comparison, the relative error is defined in equation (36).

We observe that the BAR, on average, has an error $\epsilon_{ref}$ that diminishes from 24% to 8% as the integration progresses. This error evolution is specific to the $B = 20$ replicas, as the error is measured from the average of $R = 100$ realizations. If the number of replicas were to increase, the average estimation would become more precise (including the bias of the BAR estimator) and the variability would diminish. By increasing the number of replicas, we decrease the noise in the repeated estimations of the RMSE. In figure (7) the
<figure><img src="image_8.png" /><figcaption>(a) Time evolution of the different RMSE estimators for the OU process. The exact time evolution is plotted in solid black and the ACF based estimator of equation (33) in dashed, serves as a reference to the SotA. Error bars using $\hat{\sigma}_{\hat{\sigma}_{BAR}}$ from equation (35) shown in red.</figcaption></figure>

<caption>(b) Time evolution of the error $\epsilon_{ref}$ between the average (in dashed-dot black) BAR estimate and the analytic RMSE.</caption>

Figure 6: Comparison of the $\hat{\sigma}_{BAR}$, $\hat{\sigma}_{ACF}$ RMSE estimates using a fixed $B = 20$ number of replicas and $M = 100$ for the lags of the ACF (left) and time evolution of the error $\epsilon_{ref}$.

BAR estimator is applied using $B = 100$ replicas. Since the BAR estimator is asymptotic, it is expected to be more accurate as $N$ grows. Therefore, in the following, the first 4000 steps (200[tu]) of the simulation are truncated, to not display the estimators in their transient phase.

<figure><img src="image_9.png" /><figcaption>(a) BAR method using $B = 20$ replica sets for the estimation of the RMSE.</figcaption></figure>

<caption>(b) BAR method using $B = 100$ replica sets for the estimation of the RMSE.</caption>

<caption>(c) Time evolution of the error $\epsilon_{ref}$ between the average (in dashed-dot black) BAR estimate and the analytic RMSE for $B = 1000$.</caption>

Figure 7: Comparison of the BAR RMSE estimates after increasing the number of replicas. The average behavior of the RMSE estimation does not get significantly affected by the number of replicas used, although the noise in the estimation is significantly reduced. The error between the average $\epsilon_{ref}$ diminishes in a similar rate to the case of $B = 20$.

The relative error $\epsilon_{ref}$ on the average, remains in comparable levels, from 20% to 11%, however the underestimation seems persistent. The negative bias of the $\hat{\sigma}_{BAR}$ could be attributed to the miss-match of the autocorrelation structure of the process, by the bootstrapping weights, using the specific $\beta = \sqrt{2} - 1$. More on the importance of $\beta$ can be found in Appendix A. The underestimation is compensated from the fact that the reference value is largely covered in a $1\hat{\sigma}_{\hat{\sigma}_{BAR}}$ error and at most times in a $2\hat{\sigma}_{\hat{\sigma}_{BAR}}$ error band. Finally, in figure (8), We compare the RMSE $\hat{\sigma}_{BAR}$ to the RMSE $\hat{\sigma}_{ACF}$ of equation 33 using $M=100$ and $M=1000$ lags respectively to demonstrate the sensitivity in the lag length $M$, of the ACF based method:

The $\hat{\sigma}_{ACF}$ that includes tail noise in the $\tau_I$ deviates from the reference as much as the $\hat{\sigma}_{BAR}$. This underlines the problem that in a realistic scenario, where $M$ would be unknown, the process of tuning $M$ in an ACF based estimation setting, would render the classical estimator infeasible in a large CFD simulation. In such a case, re-performing the simulation would be costly and storing every snapshot of the field to investigate the structure of the ACF is generally undesired.

# 5. Application in Error Monitoring and Error Estimation in Large Eddy Simulation

The motivation for the implementation of the current work, comes from the need of the scientific and industrial community to have an in-situ measure of convergence of a large scale simulation and also be able to monitor the quality of the time averaged quantities directly. To these ends, we developed and integrated
<figure><img src="image_10.png" /><figcaption>Figure 8: Comparison of ACF estimates with different lag lengths. On the left, *M* = 100 that is close to the first zero crossing (most accurate estimate of the decorrelation time) and on the right *M* = 1000 which includes the tail noise, making the decorrelation time estimate uncertain.</figcaption></figure>

in the YALES2 solver [32] the online bootstrapping method to produce fields of the MSE of the velocity field time average available at **each** iteration of the simulation. The entire process is done *online*. The simulation code YALES2 is a massively parallel LES code utilizing unstructured meshes. The code solves the incompressible low-Mach number Navier-Stokes equations for turbulent flows using a projection method to handle the pressure velocity coupling. The code is based on central finite volume schemes that relieve on efficient linear solvers [33]. The temporal integration of the convective terms is performed explicitly using a fourth order Runge-Kutta scheme while for the diffusive terms a semi implicit scheme is used. For each simulation presented in this work, an adequately adapted unstructured mesh is used, following the methodologies of Grenouilloux et al. [34] and the general mesh adaptation methodology of Benard et al. [35]. In this process the LES equations are integrated in intervals to converge statistics of the flow (mean and fluctuating fields). These statistics are used to directly control the local mesh size $\Delta$ through two criteria to guarantee: (i) an accurate resolution of the velocity mean field (ii) that enough turbulent scales are adequately resolved [35]. In the following, two flow configurations are being considered, a turbulent pipe flow and a turbulent round jet.

### 5.1. Turbulent Pipe flow

The proposed methodology is first applied to a turbulent pipe flow with a Reynolds number based on the bulk velocity $u_b$ and the pipe diameter $D$ defined from $Re = \frac{U_b D}{\nu} = 5300$. The flow configuration is similar to the case performed by Grenouilloux et al. [34]. We enforce periodic boundary conditions and a forcing term to ensure a constant flow rate. The adapted mesh used for our simulation consists of approximately $N_{ele} = 4.7$ million elements. The DNS data of Wu and Moin [36] are used to verify the validity of the statistics. In figure (9) we show the instantaneous field at the end of the simulation (after $T=600[D/u_b]$) and the acquired mean field who's RMSE we want to quantify.

<figure><img src="image_11.png" /><figcaption>Figure 9: Pipe flow velocity field (instantaneous top and mean bottom) with probe location at ($\frac{D}{2}, \frac{D}{2}, \frac{7.5D}{2}$)</figcaption></figure>

We start by describing the implementation of the BAR method for the estimation of the MSE $\hat{\sigma}_{BAR}$ of
the estimated mean field $\bar{u}_N$ by the SME. We directly implement the steps described in algorithm (1) in the code YALES2, since the implementation of the method is quite simple and non intrusive to the main solver. More specifically, we need to define just the B replica fields for each quantity we seek to analyse. Every step is performed online and at the end of each iteration $n$, we get a field of the RMSE $\hat{\sigma}_{BAR}$ as seen in figure (10). At each step, we loop over the entire mesh and at each node loop over the B replicas. First we generate a different $\zeta \sim N(0, 1)$ for each replica, update the weights $w^{(b)}$ in every node according to the autoregression in step 7, then update the B replicas for each node of the mesh and finally we update the mean weight $\hat{w}^{(b)}$. After completing the loop over B at a specific node, we can calculate explicitly the RMSE of $\{\hat{U}_{(b)}^{*}\}_{b \le B}$ as:

$$ \hat{\sigma}_{BAR} = \sqrt{\frac{1}{B-1} \sum_{b}^{B} (U_{(b)}^{*} - \langle U(t) \rangle)^2}, \quad (37) $$

where the $\langle U(t) \rangle$ is the time average as computed by the code. Then we repeat the process for each node. The result of this is an entire field of the RMSE of the SME of the velocity, as seen in figure (10).

<figure><img src="image_12.png" /><figcaption>Figure 10: Field of error of the SME at $T = 600[D/u_b]$. The entire field of the RMSE of the mean velocity is available at each iteration $n$ of the simulation.</figcaption></figure>

We select $B = 20$ and we advance the simulation for a total of $N = 6 \cdot 10^5$ iterations each with a time step $\Delta t = 0.01$. The parameter $\beta = \sqrt{2} - 1$ is selected as before. We then record from a probe, as indicated in figure (9), the time evolution of the velocity in a specific location of the pipe. We then intend to demonstrate the application of the BAR method, in comparison with the *off line* version of the ACF estimation method in order to assess its performance, taking into account the ease of the implementation of the method. We also provide error bars on the BAR estimate for a fixed number of replicas $B$ as before.

The following are performed *offline*. We use this stored time series, for an posteriori analysis of the BAR method. We will perform two different calculations using the stored time series. First we want to provide an *error bar* for the estimated $\hat{\sigma}_{BAR}$. The BAR method is a variance estimation that comes from the variance of replicated values of the mean field that are constructed using a random weight $w$ following the auto-regression of equation (26). Therefore, each realization of the BAR method is independent and is performed by using the same time series recorded the point in the domain. We repeat the BAR calculation for a total of $R = 100$ realizations using the same number of replicas $B = 20$ and then calculate the mean and STD on these realizations.

Now we use the same recorded time series data from the probe, to estimate the variance of the SME using the classic ACF based approach. As in the OU process, the time series have an autocorrelation structure that corresponds to an integral time scale $\tau_I$. This is the quantity that the state of the art methods implicitly or explicitly estimate, in order to estimate the variance of the SME. This time is given by the integral of the ACF, theoretically over the entire lag time. However since this is impossible due to the finite size of the sample set, one must select carefully the number of lags that contribute physical information to the $\tau_I$. In terms of integer lags, this is equivalent to inquiring until which point do we need to integrate the autocorrelation function. This is where the problem with the state of the art (that resorts to models) lies: the estimation of the necessary lag length has to be done posteriori and also depends on the physics (i.e. the location of the probe in the domain) and the time step. As a solution to the problem of the lag selection we seek the first zero crossing of the ACF [11]. We initially select an arbitrary number of lags (M=1000)
and compare the ACF against the ACF where all the available data are considered. This way we verify our estimator ($\hat{\sigma}_{ACF}$) is as accurate as possible, since truncating the number of lags used for the estimation only truncates the ACF.

<figure><img src="image_13.png" /><figcaption>(a) Time evolution of $\hat{\sigma}_{ACF}$, $\hat{\sigma}_{BAR}$, $\sigma_N$. For the BAR estimator, a $3\sigma$ error bar is presented along with the mean, by evaluating the variance from $R = 50$ realizations of the BAR method in post processing for $B = 20$. In black is the estimate of the RMSE produced by YALES2.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Autocorrelation function of the velocity time series sampled from the point $\mathbf{p} = (0, 0, 7.5)$ using $m = 1000$, first zero crossing and ACF calculated using the entire timeseries used for the estimation of the $\hat{\sigma}_{ACF}$ in figure (11a).</figcaption></figure>

<caption>Figure 11: Time evolution of the different RMSE estimators for the pipe flow, from data sampled from</caption>

From figure (11b) we observe that the estimator used for the $\hat{\sigma}_{ACF}$ doesn't include too much of the tail noise of the ACF in the estimation of the integral time scale, therefore we can use that as a conservative estimate (overestimation) of the "true" variance of the estimator. We verify the necessary lag length, by detecting the first zero crossing of the estimated ACF, at $M = 160$ as can be seen from figure (11b). We now will compare the time evolution of $\hat{\sigma}_{BAR}$, $\hat{\sigma}_{ACF}$. We also provide the mean RMSE with $3\sigma$ error bars, from the BAR method calculated during post processing $\hat{\sigma}_{BAR,OF}$ using $R = 100$ realizations of the BAR method. (Figure 11a). We observe that the $\hat{\sigma}_{BAR}$ predicted by the BAR method is consistently smaller than the one estimated using the ACF. However, they can only be meaningfully compared by taking into account some physical characteristic of the flow, like the degree of isotropy of the turbulence. The $\hat{\sigma}_{ACF}$ contains itself a bias stemming from the ACF integration, i.e. the local decorrelation time that is dependent on the physics of the flow. Therefore, a coverage of the ACF in the $2\sigma$ error of the BAR estimator, is deemed satisfactory for the application of the BAR method to a complex case. We can use this estimated error of the SME to monitor the convergence of a profile of the mean in any location of the flow. We provide comparison to DNS data given by Wu and Moin [36] for a line at $z = L/2$ of the domain.

<figure><img src="image_15.png" /><figcaption>(a) Profile of the mean velocity field at a location $z = 5$ after $10s$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Profile of the mean velocity field at a location $z = 5$ after $600s$</figcaption></figure>

<caption>Figure 12: Mean velocity field profiles with time averaging error bounds at different total simulation times</caption>

As the integration time increases, naturally the error due to time integration diminishes and the deviations from the DNS can be more clearly attributed to other sources like the spatial discretisation or the choice of subgrid scale model. This kind of visual quality inspection of a simulation can be now performed easily, directly on any location of the domain, opening the way for a more streamlined workflow for CFD
engineers. The BAR method is an asymptotic bootstrap estimator, therefore there are two "costs" that we pay for the easy implementation and integration in any flow solver. The first is the number of replica sets B whose increase has a positive correlation with the accuracy of the method. Evidently, the available computational resources should be spent to select a B as large as possible. The other one, is the fact that the BAR estimator should become progressively more accurate as we acquire more samples (time steps) from the simulation, being more accurate, later in the run. However these are small prices to pay for a general model that does not require any significant user knowledge of the flow and is flexible enough to be applied to different physical settings, as is. Finally, even the B = 20 number of replicas give a satisfactory coverage of the $\hat{\sigma}_{ACF}$ estimate.

## 5.2. MSE of time averages of LES simulations: Jet flow

As a last demonstration case, we apply our methodology to a round jet flow similar to the works of Grenouilloux et al. [34]. This configuration is defined by its inlet profile. The streamwise component of the velocity is given by a hyperbolic tangent profile, according to the work of [37], that depends on the radius of the jet $R$, a centerline velocity $U_1$ and a co-flow velocity $U_2$. We fix the Reynolds number at $Re = (U_1 - U_2)\frac{R}{\nu} = 4,000$. The computational domain has a size of $30R \times 28R \times 28R$ in the streamwise and the two transverse directions. We run the simulation for a duration of $600R/(U_1 - U_2)$ and with a constant time step of $\Delta t = 0.05$. The methodology followed for the application of the BAR method is similar to the pipe flow case. We select $B = 20$ and perform a simulation with $N = 120000$ iterations. In this simulation we update at every step the RMSE of the SME. For this reason we allocate just B fields for the replicas $\hat{U}_{(b)}^*$. The fields $\hat{U}_{(b)}^*$ are updated at every iteration $n$ and are available for inspection by the use, where the coefficients $w^{(b)}, \bar{w}^{(b)}$ are being used for the entire snapshot of the flow.

<figure><img src="image_17.png" /><figcaption>Figure 13: Instantaneous, mean and RMSE fields for different transverse locations in the round jet configuration.</figcaption></figure>

In areas near the core of the jet where the velocity field has smaller turbulent intensity, we observe smaller RMSE whereas around the core and on the shear layer where the turbulence is more intense, after the jet break-down due to Kelvin Helmholtz instabilities, the RMSE is higher.

Every step for the update of the BAR method is performed online. As in the pipe flow in section (5.1), we perform an offline calculation of the BAR method, for all sampled locations of the flow. We sample from the flow in 5 different points in the domain to demonstrate the change of the ACF due to the differences in the physics of the solution. The first two, near the shear layer at the entrance of the jet in the domain and the rest are scattered in regions where the jet has broken down. We perform the offline calculations on this specific set of probes to provide errors on the BAR method estimate and to reproduce the classical ACF based method that would not be possible online. In the jet case, the variability in the integral time scales $\tau_I$ is more evident due to the larger variability in the auto covariance of the flow throughout the domain due to changes in turbulent intensity and isotropy. Specifically, the ACFs from points 1 and 2, differ significantly from the rest, due to persisting coherent structures. Initially we estimate the ACF using $M = 1000$ lags. Then for the ACF based estimation of the MSE, we will use again the truncation criterion of the "first zero crossing". As we can see from figure (15), in the first two points, the ACF the truncation criterion of the first zero crossing is omitting information from the oscillations in the tail of the ACF before the numerical
<figure><img src="image_18.png" /><figcaption>Figure 14: Instantaneous, mean and RMSE fields for a single slice of the computational domain (a) and the location of the 5 probes used to collect data from the simulation.</figcaption></figure>

noise. In the shear layer the flow is displaying large coherent structures, which influence the ACF, that are a result of the transitional instabilities.

<figure><img src="image_19.png" /><figcaption>Figure 15: Autocorrelation functions (ACFs) of flow velocity sampled at various points in the domain.</figcaption></figure>

In the rest of the points in the domain, the ACF doesn't display such intensely oscillating tails, but the first zero crossing appears in larger lags. At each point, we compare the RMSE estimated using the ACF (entirely in post processing) and using the BAR method (completely online). The realization of the BAR method as recorded by the code ($\hat{\sigma}_{YALES}$) is shown by a dashed black line, in the dot-marked line the $\hat{\sigma}_{ACF}$ and in black dash-dotted line is the mean RMSE of $R = 100$ realizations of the method offline, in order to quantify its error due to the finite number of replicas $B = 20$.

A smaller distance between $\hat{\sigma}_{ACF}$, $\hat{\sigma}_{BAR}$ is observed around the shear layer, in points 1 and 2 from figure (14b), in contrast to where the jet has transitioned to more isotropic turbulent mixing. This could be attributed to the discrepancies in estimating the integral time scale in each case as well as inadequate correlation in the weights $w_b$. Early truncation, by the first zero crossing rule, in the shear layer probes where the ACFs display periodic tails can be less severe compared to early truncation in the jet breakdown region where correlations persist for longer. We remind the reader, that in the case of the LES runs, we can not have an absolute reference, since also the $\hat{\sigma}_{ACF}$ is biased in an *unknown* direction, therefore allowing
<figure><img src="image_20.png" /><figcaption>Figure 16: Time evolution of the different estimators ($\hat{\sigma}_{ACF}$, $\hat{\sigma}_{BAR}$, $\hat{\sigma}_{YALES}$) for the RMSE at the 5 points of the jet domain.</figcaption></figure>

only for qualitative comparisons.

### 5.3. *Implementation and computational cost of the method*

To assess the computational overheads of the method, we ran a series of pipe flow simulations with different number of replicas *B*. The method requires only storage of the *B* replica mean fields. Since we are not calculating spatial correlations, we can imagine the same $w^{(b)}$, $\bar{w}^{(b)}$ being used for the entire snapshot of the field. Therefore the only real cost is the *B* replicas $\hat{U}^*$ that we need for equation (27). We

<figure><img src="image_21.png" /><figcaption>Figure 17: Overall computational overhead factors for the pipe flow simulation. Values are normalized by the cost without using the method.</figcaption></figure>

extract the time spent for a single iteration of the simulation, on average from the simulation logs. We observe an overhead of maximum 1.7 times of the iteration's cost, in the scenario when one uses *B* = 80 which could be considered excessive for our applications. On average the increase in *B* seems to affect the
computational cost logarithmically. For the memory scaling we observe a totally linear scaling, making it evident that the first and possibly only limiting factor for the application of our method is a large number of replicas. Evidently, the quality of the estimation is analogous to the number of replicas used and therefore the practitioner should allocate the available resources to storing more replicas. However as observed already by our experiments, using just $B = 20$ replicas suffices for a $2 \ 3\sigma$ coverage of the ACF based estimator, making larger $B$ possibly unnecessary.

## 6. Conclusions and discussion

In figure (18) we observe the spatial distribution of the error of the SME at $T = 600$ for the pipe flow and the jet flow cases. This is, in our knowledge, the first time we can extract an entire field of error of the SME at *any* time of the simulation without any significant computational overhead, post-processing or a-priori flow physics knowledge (that would be required for batch based methods or for selecting subsampling times for example) enabling the use of this field as a monitoring quantity to take decisions regarding the state of the simulation *in a completely dynamic* way.

<figure><img src="image_22.png" /><figcaption>Figure 18: Fields of the estimated $\hat{\sigma}_{BAR}$ for the pipe and jet flow at $T = 600[\text{tu}]$.</figcaption></figure>

In this paper we propose the use of a newly developed online bootstrapping scheme to estimate errors in time averaged turbulent flow quantities. We demonstrate the equivalent performance with the state of the art that involves estimating the ACF. Compared to the state of the art, our proposed method requires no assumptions of the flow quantity it is used on and compared to the existing methods, so far, has a highly simple and easy implementation in large codes, making it an attractive alternative to current ACF estimation/fitting based methods. We presented the time evolution of the estimated standard deviation of the time-average of a turbulent flow variable and of the time evolution of the OU process. The RMSE was estimated directly, by estimating the ACF and using it's integral to correct the sample Standard Deviation. This direct calculation of the Autocorrelation coefficient was compared to the method we propose implementing, an On-Line Bootstrap with Auto-Regressive weights (BAR). The method was tested on synthetic data for which a reference $\hat{\sigma}$ is available for comparison and then applied it to two LES cases, a pipe and a jet flow. The behaviour of predicting, in general, smaller values than the ACF based method, is observed throughout all the tests performed. The systematic estimation of a lower value of the error of the SME can be attributed to the method only being asymptotically consistent meaning that the accuracy of the method increases with simulation time and to the fact that the optimal value for the "correlation persistence" $\beta$ is not capturing accurately enough the temporal autocorrelation structure of our time series. Furthermore, the increase in $B$ should reduce the estimator's sampling error, however the decrease is small for relative large increases in $B$ (five-fold), justifying even a small number of replicas. Overall, we deem this underestimation not significant, given the method's ease of implementation and robustness to the number of replicas. Finally, we successfully provide an extension to the non constant time step case, that has, to our knowledge, never been previously explored. Gathering the above results, we want to stress the importance of the difference between our proposed method and the existing ACF estimation methods that are based on a
specific time series lag length. With the proposed method, the calculation of the variance is dependent only on the number of replicas and not on any parameter that requires an a posteriori treatment to select. This a posteriori treatment was necessary up to some degree to all the currently proposed methods in the literature, even the closest to a truly online setting [22]. The motivation behind proposing a method to easily estimate the error of a time averaged turbulent quantity is having a way to monitor the quality of the simulation at every point of the domain. With these error estimates we plan to propose decision making strategies for deciding when to perform crucial actions on a simulation like a mesh adaptation or the termination of the simulation. Both of these, contribute positively to having a more informed way of distributing allocated computational time, enabling the broadening of the applications of advanced simulations of turbulent flows, like LES.

## Acknowledgements

This work was granted access to the HPC resources of CINES/TGCC/IDRIS under the project 2A00611 made by GENCI. Part of this work has been initiated during the Extreme CFD Workshop Hackathon (https://ecfd.coria-cfd.fr).

## Funding

The authors gratefully acknowledge support from NETHUNS project under grant ANR-21-CHIN-0001-01.

## Appendix A. Sensitivity of the BAR method to hyperparameters

The two hyperparameters of the BAR method are the number of replica sets B and a "correlation decay/persistence" coefficient β that appears in equation (26) in $ρ_n = 1 - n^{-β}$. The parameter β controls the speed of convergence of the MSE of the estimated RMSE through the BAR method, by controlling how correlated the sequence {$w_n^{(b)}$}$_{n≤N}$ is. The choice of $β = \sqrt{2} - 1$ comes from balancing the order of convergence of the bias of the estimation and the variance of the estimate. Specifically as shown in [10], if the target infinite time variance is

$$ \sigma_{\infty}^{2} = \lim_{n \to \infty} \mathbb{V} \left[ \frac{1}{\sqrt{n}} \sum_{i}^{n} u_{i} \right], \qquad (A.1) $$

then the mean square error of the estimation scales as:

$$ \mathcal{O} \left( \mathbb{E} \left[ \mathbb{V} \left[ \frac{1}{\sqrt{n}} \sum_{i}^{N} w_i u_i \right] - \sigma_{\infty}^2 \right]^2 \right) = \underbrace{\mathcal{O}(n^{-\frac{\beta}{\beta+1}})}_{\text{bias}} + \underbrace{\mathcal{O}(n^{\beta-1})}_{\text{Variance}} \qquad (A.2) $$

selecting $β_{opt} = \sqrt{2} - 1$ ensures the bias term of (A.2) converges in the same rate as the variance term of the MSE between the BAR estimate and $σ_{\infty}^2$. Selecting a small $β$ implies a very slow convergence of the bias term, or in other words very severe underestimation of the RMSE (something undesirable). We study the sensitivity of the method on $β$ on the OU process, where we re-perform the estimation for a small value of $β = 0.01$. We again provide statistics for the estimation performed on $R = 50$ repetitions.

From figure (A.19) we can see that the effect of decreasing $β$ from its optimal value, slows down the convergence of the bias term of equation (A.2) as now the order becomes $\underbrace{\mathcal{O}(n^{-\frac{0.01}{1.01}})}_{\text{bias}} + \underbrace{\mathcal{O}(n^{-0.99})}_{\text{Variance}}$. However for larger times, the effect of $β$ should be increasingly less strong but the damage from the underestimation should persist. Hence, we propose to use the $β = \sqrt{2} - 1$ as suggested initially, to have the same rate of convergence between the bias-variance of the estimate
<figure><img src="image_23.png" /><figcaption>Figure A.19: Time evolution of the BAR estimator against the BAR estimator with $\beta = 0.01$.</figcaption></figure>

## Appendix B. The importance of time-step variability in the estimation process.

It is often assumed in the existing literature that the spacing between the elements in the collected time series, coming from CFD simulations, are equal i.e. the data is sampled at a constant rate. However, when the simulation code is using an explicit time integrator, the time step is non constant. The time average is now the random integral:

$$ \mathcal{I}_T = \frac{1}{\sum_i \Delta t_i} \sum_i^N u_i \Delta t_i, \qquad (B.1) $$

where $Cor(u, \Delta_t) = 0$ but the $\Delta t \sim N(\bar{\Delta}t, V_{\Delta t})$. After calculations, we can show that the variance of the time average will be of the form of equation (4) plus a term related to the non uniform time step.

$$ \mathbb{V}[\mathcal{I}_T] = \mathbb{V} \left[ \sum_{i}^{N} \frac{u_i \Delta t_i}{\sum_{i}^{N} \Delta t_i} \right] = \frac{1}{N} (C_U + C_{\Delta t, U} \rho_{\Delta t}^2) \mathbb{V}[U] \qquad (B.2) $$

where $\rho_{\Delta t} = \frac{\sqrt{\Delta t}}{\Delta t}$ and $C_{\Delta t, U}$ is the correlation coefficient of the velocity with the time-step as also seen in equation (4). Therefore compared to the first case of the non varying time step, we expect an extra term due to the time step variation, contributing to the uncertainty of the time average. We can clearly see this in figure (B.20) where there is a clear bias to the estimate of the MSE due to the variability of the time step. In figure (B.20) we compare the same pipe flow problem simulated for $T = 600$ using a constant $\Delta t = 0.01$ and using a CFL condition. In each case we use our transformation of equation (15) to provide aposteriori error bands for the BAR estimation. The weight factor $z$ for $\Delta t = \text{const}$ becomes $z = \frac{N\Delta t}{N\Delta t} = 1$. The result is the clear bias, showing the importance of the extra correlation term in equation (B.2). The average $\mu_{\sigma_{BAR}}^{CFL}$ is clearly larger than $\mu_{\sigma_{BAR}}^{\text{const}}$, denoting the delay due to the extra correlations $C_{\Delta t, U} \rho_{\Delta t}^2$.

<figure><img src="image_24.png" /><figcaption>Figure B.20: Time evolution of the predicted mean RMSE with error bands calculated from $R = 50$ repetitions of the estimation. The estimation is repeated for the non constant $\Delta t$ case.</figcaption></figure>
In such a case, a comparison with the sample ACF estimator of (33) would be impossible. Firstly, because the two estimators have different sample sizes and secondly because in the non equidistant time series, consecutive points are not sampled in the same physical times as two consecutive points in the equidistant case, therefore an estimator that works for both should account for this lag.

## References

[1] S. Pope, Ten questions concerning the large-eddy simulation of turbulent flows, New Journal of Physics 6 (2004). doi:10.1088/1367-2630/6/1/035.

[2] S. B. Pope, Turbulent Flows, Cambridge University Press, 2000.

[3] J. Meyers, P. Sagaut, B. J. Geurts, Quality and reliability of large-eddy simulations, in: ERCOFTAC Series, volume 12, 2008.

[4] T. Oliver, N. Malaya, R. Ulerich, R. Moser, Estimating uncertainties in statistics computed from dns (2013).

[5] S. Rezaeiravesh, R. Vinuesa, P. Schlatter, An uncertainty-quantification framework for assessing accuracy, sensitivity, and robustness in computational fluid dynamics, Journal of Computational Science 62 (2022) 101688. doi:10.1016/J.JOCS. 2022.101688.

[6] J. P. Slotnick, A. Khodadoust, J. J. Alonso, D. L. Darmofal, W. Gropp, E. A. Lurie, D. J. Mavriplis, Cfd vision 2030 study: A path to revolutionary computational aerosciences, 2014. URL: https://api.semanticscholar.org/CorpusID: 215813979.

[7] A. Khan, H. Sim, S. S. Vazhkudai, A. R. Butt, Y. Kim, An analysis of system balance and architectural trends based on top500 supercomputers, in: The International Conference on High Performance Computing in Asia-Pacific Region, ACM, 2021, pp. 11-22. doi:10.1145/3432261.3432263.

[8] X. I. A. Yang, W. Zhang, M. Abkar, W. Anderson, Computational fluid dynamics: its carbon footprint and role in carbon emission reduction (2024).

[9] W. Wei, Time Series Analysis: Univariate and Multivariate Methods, volume 33, 1989. doi:10.2307/2289741.

[10] N. Palm, T. Nagler, An online bootstrap for time series (2023).

[11] P. O'Neill, D. Nicolaides, D. Honnery, J. Soria, Autocorrelation functions and the determination of integral length with reference to experimental and numerical data (2004).

[12] S. Russo, P. Luchini, A fast algorithm for the estimation of statistical error in dns (or experimental) time averages, Journal of Computational Physics 347 (2017) 328-340. doi:10.1016/J.JCP.2017.07.005.

[13] P. Beyhaghi, S. Alimohammadi, T. Bewley, Uncertainty quantification of the time averaging of a statistics computed from numerical simulation of turbulent flow (2018).

[14] W. K. George, Lectures in Turbulence for the 21st Century, Lecture notes, 2013.

[15] B. L. Smith, D. R. Neal, M. A. Feero, G. Richards, Assessing the limitations of effective number of samples for finding the uncertainty of the mean of correlated data, Measurement Science and Technology 29 (2018) 125304. doi:10.1088/ 1361-6501/aae91d.

[16] D. Xavier, Uncertainty quantification for time varying quantities in turbulent flows (2024). URL: https://urn.kb.se/ resolve?urn=urn:bn:se:kth:diva-342785.

[17] M. Wallerberger, Efficient estimation of autocorrelation spectra (2018).

[18] M. Papageorge, J. A. Sutton, Statistical processing and convergence of finite-record-length time-series measurements from turbulent flows, Experiments in Fluids 57 (2016). doi:10.1007/s00348-016-2211-z.

[19] A. Grenouilloux, High-Fidelity simulation of the aerothermal performances of a turbofan Thrust Reverser, Theses, Normandie Université, 2023. URL: https://theses.hal.science/tel-04538595.

[20] G. Treviño, E. L. Andreas, Dynamical implications of block averaging, Boundary-Layer Meteorology 120 (2006) 497-508. doi:10.1007/s10546-006-9060-x.

[21] D. N. Politis, The impact of bootstrap methods on time series analysis, Statistical Science 18 (2003) 219-230. URL: http://www.jstor.org/stable/3182852.

[22] S. Rezaeiravesh, G. Scheidle, A. Peplinski, J. Garcke, P. Schlatter, In-situ estimation of time-averaging uncertainties in turbulent flow simulations (2023).

[23] L. H. Benedict, R. D. Gould, Towards better uncertainty estimates for turbulence statistics, Experiments in Fluids 22 (1996) 129-136. doi:10.1007/s003480050030.

[24] P. Sagaut, Y.-T. Lee, Large eddy simulation for incompressible flows: An introduction. scientific computation series, Applied Mechanics Reviews 55 (2002) 115-. doi:10.1115/1.1508154.

[25] S. T. Salesky, M. Chamecki, N. L. Dias, Estimating the random error in eddy-covariance based fluxes and other turbulence statistics: The filtering method, Boundary-Layer Meteorology 144 (2012) 113-135. doi:10.1007/s10546-012-9710-0.

[26] J. L. Horowitz, Bootstrap methods in econometrics, Annual Review of Economics 11 (2019) 193-224. doi:10.1146/ annurev-economics-080218-025651.

[27] H. Lam, A cheap bootstrap method for fast inference (2022).

[28] B. Efron, Bootstrap methods: Another look at the jackknife, The Annals of Statistics 7 (1979). doi:10.1214/ aos/ 1176344552.

[29] S. Gonçalves, D. Politis, Discussion: Bootstrap methods for dependent data: A review, 2011. doi:10.1016/j.jkss.2011. 07.003.

[30] E. Antal, Y. Tillé, Simple random sampling with over-replacement, Journal of Statistical Planning and Inference 141 (2011) 597-601. doi:10.1016/j.jspi.2010.06.029.
[31] A. W. van der Vaart, J. A. Wellner, Weak Convergence and Empirical Processes, Springer International Publishing, 2023.
doi:10.1007/978-3-031-29040-4.

[32] V. Moureau, P. Domingo, L. Vervisch, Design of a massively parallel cfd code for complex geometries, Comptes Rendus
Mécanique 339 (2010) 141–148. doi:10.1016/j.crme.2010.12.001.

[33] M. Malandain, N. Maheu, V. Moureau, Optimization of the deflated conjugate gradient algorithm for the solving of elliptic equations on massively parallel machines, Journal of Computational Physics 238 (2013) 32–47. doi:10.1016/j.JCP.2012.11.046.

[34] A. Grenouilloux, G. Balarac, J. Leparoux, V. Moureau, G. Lartigue, P. Bénard, R. Mercier, P. Ferrey, On the use of kinetic-energy balance for the feature-based mesh adaptation applied to large-eddy simulation in complex geometries, in: Volume 10C: Turbomachinery Design Methods and CFD Modeling for Turbomachinery; Ducts, Noise, and Component Interactions, American Society of Mechanical Engineers, 2022. doi:10.1115/GT2022-80315.

[35] P. Bénard, G. Balarac, V. Moureau, C. Dobrzynski, G. Lartigue, Y. DAngelo, Mesh adaptation for large-eddy simulations
in complex geometries, International Journal for Numerical Methods in Fluids 81 (2015) n/a–n/a. doi:10.1002/fld.4204.

[36] X. WU, P. MOIN, A direct numerical simulation study on the mean velocity characteristics in turbulent pipe flow, Journal
of Fluid Mechanics 608 (2008) 81–112. doi:10.1017/S0022112008002085.

[37] A. Michalke, G. Hermann, On the inviscid instability of a circular jet with external flow, Journal of Fluid Mechanics 114
(1982) 343. doi:10.1017/S0022112082000196.