12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022

ESTIMATING UNCERTAINTY OF LOW- AND HIGH-ORDER TURBULENCESTATISTICS IN WALL TURBULENCE

Saleh Rezaeiravesh, Donnatella Xavier & Ricardo Vinuesa
Jie Yao
Department of Engineering Mechanics
Department of Mechanical Engineering,
KTH Royal Institute of Technology
Texas Tech University
Stockholm, Sweden
Lubbock, Texas, USA{salehr,dgxavier,rvinuesa}@mech.kth.sejie.yao@ttu.edu

Fazle Hussain
Department of Mechanical Engineering,Texas Tech University
Lubbock, Texas, USAfazle.hussain@ttu.edu

Philipp Schlatter
Department of Engineering Mechanics
KTH Royal Institute of Technology
Stockholm, Sweden
pschlatt@mech.kth.se

ABSTRACT

1. Batch-based methods, which are based on batching thetime series and acting on the batch means, examples are nonoverlapping/overlapping batch means (NOBM/OBM) methods and batch-means and batch-correlations method (BMBC)
introduced by Russo & Luchini (2017). 2. Autocorrelationfunction (ACF)-based estimators, which rely on evaluating theexact analytical expression for σ 2 (µ̂):

A framework is introduced for accurate estimation oftime-average uncertainties in various types of turbulence
statistics. A thorough set of guidelines is provided to adjust thedifferent hyperparameters for estimating uncertainty in sample
mean estimators (SMEs). For high-order turbulence statistics,
a novel approach is proposed which avoids any linearizationand preserves all relevant temporal and spatial correlations andcross-covariances between SMEs. This approach is able to accurately estimate uncertainties in any arbitrary statistical moment. The usability of the approach is demonstrated by applying it to data from direct numerical simulation (DNS) of
the turbulent flow over a periodic hill and through a straight
circular pipe.


 
n−1 
k
1
γ0 + 2 ∑ 1 −γ ,σ (µ̂) =
n
n k
k=1
2

where γk = E[(ui − µ̂)(ui+k − µ̂)] is the autocovariance for thestatistically-stationary sequence {ui }ni=1 at k = 0, 1, 2, . . . timelags. As it can be inferred from Eq. (2), for a given number ofsamples two main factors which can give rise to the variance ofan SME are the large variance of time series and the slow decay of the autocorrelation γ(k)/γ(0). These two factors can belarge either due to the physics of turbulence or numerical artifacts. An illustration of the latter can be found in Rezaeiravesh
et al. (2022), where it was shown that the spatial location along
the profile of a turbulence quantity where the sensitivity withrespect to the numerical parameters is high, the time-average
uncertainty can be large as well.
The main barrier to accurate estimation of σ 2 (µ̂) via
Eq. (2) is the oscillations in the sample-estimated autocovariances γ̂(k) at large lags. As first shown by Oliver et al.
(2014), autoregressive models (ARM) are suitable for constructing smooth models for the ACF in turbulent flows (hereafter, ARM-based uncertainty estimator). An alternative toconstructing an ARM is to fit an ad-hoc function to sample
ACF, see e.g. Gscheidle et al. (2022). In our recent study,
Xavier et al. (2022), the performance of various estimatorswith respect to associated hyperparameters was thoroughly investigated for wall turbulence. For batch-based variance estimators, the hyperparameter is the batch size, and we observedthat for the NOBM, it could be determined from the lag-1 temporal correlation of the batch means. For the BMBC estimator,
the lag-2 and lag-1 correlations of the batch means are both

TIME-AVERAGE UNCERTAINTY
One of the main goals of any turbulent flow simulation
is to compute various flow statistics after the statisticallystationary condition is achieved. Recent progress in hardware
and software has made the scale-resolving simulations of turbulent flows more accessible. However, any such simulation
can only be run for a finite time; therefore, the sample mean
estimator (SME) µ̂ = Ê[u] for computing the time-average ofa turbulent quantity u is uncertain. The SMEs are unbiased,
and by ergodicity and the central limit theorem,µ̂ ∼ Nµ, σ 2 (µ̂)

(2)

(1)

where N is a Gaussian distribution with true mean µ = E[u]and variance σ 2 (µ̂). The latter is a direct measure of the uncertainty in µ̂ in a single realization of a turbulent flow where theassociated µ cannot be realized. Since the time series of turbulence quantities are highly correlated, for estimating σ 2 (µ̂)
standard linear unbiased estimators typically used for independent and identically distributed (iid) samples are not accurate.
The remedy can be found in the field of uncertainty quantification (UQ) where the tools for estimating σ 2 (µ̂) for autocorrelated time samples can be divided into two main categories:

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022
ing exact expression can be derived,

needed for optimal estimation of batch size and hence accurate estimation of σ 2 (µ̂). For modeled-ACF-based estimators,
it is shown that the order of the ARM could be obtained fromthe turbulence integral time scale, without the need for any
mathematical order selection criteria. By constraining the hyperparameters using integral quantities (such as eddy turnovertimes), guidelines were derived by Xavier et al. (2022) for accurate estimation of uncertainty in SMEs of primitive variables
in turbulent flows.

cov(µ̂i , µ̂j ) =where γi j (k) := cov(uil , ujl+k ) is the cross-covariance between
the time-series associated with two quantities ui and u j . Notethat the only symmetry is γi j (k) = γ ji (−k). As a convenience,
if k > 0 in γi j (k), then k is a lead for ui and a lag for u j . Also,
the peak of γi j (k) does not necessarily happen at k = 0. Notethat for i = j, Eq. (4) is reduced to Eq. (2). Similar to Eq. (2),the sample-estimated values of γi j (±k)/γi j (0) can be oscillatory at high k, therefore, it is necessary to develop methods toconstruct models for these cross-covariances which are smoothfunctions of leads/lags. This can, in general, be more challenging compared to the modeling of γ(k) in Eq. (2). An approach
could be using an algebraic function as introduced in Gscheidle et al. (2022). A workaround, which can however be lessaccurate, is to cut off the summation in (4) at a k above which
the sample-estimated γi j (±k) are oscillatory.
To summarize, the steps of the algorithm for estimating
uncertainty in T = f (M) are the following:

HIGH-ORDER & NONLINEAR STATISTICS
The methods reviewed in the previous section can only beapplied to estimate uncertainty in a single sample mean estimator. But, turbulence statistics include higher-order moments
of primitive variables, nonlinear terms and different budget
terms in the transport equations of Reynolds-stress components and turbulent kinetic energy. As outlined in this section, the present study proposes a UQ approach to accurately
estimate time-average uncertainty in any turbulence statistics.
Using combinations of arithmetic and differentiation operators, any of such statistical terms can be written in terms ofthe moments of the primitive variables. For instance, thesecond-order central moment ⟨u′ v′ ⟩ can be written as ⟨u′ v′ ⟩ =⟨uv⟩ − ⟨u⟩⟨v⟩. Hereafter, we denote a compound statisticalterm such as ⟨u′ v′ ⟩ by T which is defined in terms of the firstorder moments M which in this case are M = {⟨u⟩, ⟨v⟩, ⟨uv⟩}.
In general, decomposition of T in terms of M should be in a
µ = {µ̂1 , µ̂2 , · · · , µ̂ p } which
way that the uncertainty in SMEs µ̂
correspond to the members of M can be estimated by any of
the approaches of the previous section.
The problem to tackle is defined as estimating uncertainty
in T = f (M) where f (·) is a given functional, from estimateduncertainty in the SMEs of the members of M. Clearly, thisis a UQ forward (uncertainty propagation) problem, (Smith,2013). For better understanding the features of the proposed
approach, a brief review of the relevant studies for estimating
uncertainty in compound statistics T is given first. The resampling methods exercised by Benedict & Gould (1996) do
not consider the fact that turbulence time-series samples areautocorrelated; therefore, its resulting uncertainty estimations
can be inaccurate. For the turbulent channel flow, Hoyas &Jiménez (2008) applied the non-overlapping batch method and
Lee & Moser (2015) employed the ARM-based uncertainty estimator. However, these studies do not provide any detail on
how uncertainties in high-order statistics are estimated.
Our proposed UQ approach has a number of distinguishing characteristics. First, all the arithmetic and differentiation
operations corresponding to f (·) are applied in a mathematically and statistically consistent way while avoiding any approximation such as linearization. As a key feature, the correlations between different SMEs appearing in a compound term
as well as the spatial correlations in the SMEs are accurately
incorporated. Representing the covariance matrix associatedµ = {µ̂1 , µ̂2 , · · · , µ̂ p } by Σ , Eq. (1) can be generwith SMEs µ̂
alized as,µ ∼ N (µµ , Σ) ,
µ̂

#
"

n−1 

1
k
γij (0) + ∑ 1 −
γij (k) + γij (−k) (, 4)nnk=1µ corresponding to the SMEs of
1. Estimate uncertainty in µ̂
the members of M using any of the approaches in the previous section;
2. Estimate the off-diagonal elements of Σ from Eq. (4);µ , directly draw N joint
3. Using Σ and approximating µ by µ̂
samples from Eq. (3) and estimate T (denoted by T̂ ). Notethat this is a pure Monte Carlo (MC) sampling method for
the UQ forward problem;
4. For sufficiently large N, use sample estimator for meanand standard-deviation of T̂ . The latter is used to computea confidence interval for the mean of T̂ .
The only approximation involved in this algorithm is in step 3,µ . If the number of time samples forwhere µ is replaced by µ̂µ is sufficiently high, the potentially induced error
estimating µ̂
in the estimated variance of T̂ is negligible.
The above algorithm along with various time-average uncertainty estimators including those in the previous sections
have been implemented in UQit, (Rezaeiravesh et al., 2021)
and will be published as open-source.

VALIDATION OF THE UQ ALGORITHM
The proposed UQ algorithm aims at estimating the uncertainty of the statistics in a single realization of any dynamicalsystem when the statistically-stationary conditions hold. To
validate the algorithm, its estimated uncertainties in various
statistical terms can be compared to the corresponding ensemble uncertainties obtained from repeating a simulation (multi
realizations). If a simulation is repeated Ne times (Ne is also
called the size of ensemble) while keeping all simulations independent of each other, then the empirical ensemble mean
and variance of a sample-mean estimation µ̂ are respectively
obtained from the following two expressions:(3)

where the Gaussian distribution is p-variate. The diagonal elements of Σ are σ 2 (µ̂i ) for i = 1, 2, . . . , p and can be estimatedby any of the methods in the previous section. The estimation of the off-diagonal elements which are the covariance between the SMEs has, to our knowledge, not been previously
addressed in the literature. For these covariances the follow-µe = E[µ̂] ≈1 Ne∑ µ̂i ,
Ne i=1(5)σe2 = V[µ̂] ≈1 Ne∑ (µ̂i − µe )2 .
Ne i=1(6)

Recall that µ̂ is the sample mean estimation of a singlerealization time-series. We can design a model problem which212th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022
evaluated a-posteriori. Similarly, ⟨x′3 ⟩ can be estimated from⟨x′3 ⟩ ≈ µx3 + 2µx3 − 3µx µx2 + εµ̂x′3 with the associated uncertainty expressed by,

can be inexpensively simulated an arbitrary number of times
and for an arbitrary sample size per realization. At the sametime, the samples of such problem need to be autocorrelated
as in a statically-stationary turbulence time-series. In particular, we consider a first-order vector autoregressive model togenerate samples for two generic variables x and y:    

  
xi
a
a a
xi−1
ε
= 0 + 1 2
+ xi ,
yib0
b1 b2 yi−1εyiεµ̂x′3 = εµ̂x3 + 2εµ̂3x − 3εµ̂x εµ̂x2 + 6µx2 εµ̂x + 6µx εµ̂2x −
− 3 µx εµ̂x2 + µx2 εµ̂x ,(9)(7)µ = [µ̂x , µ̂x2 , µ̂x3 ]T . The two middle plots in Figure 1 repwith µ̂
resent the PDF of the uncertainty in ⟨x′ y′ ⟩ and ⟨x′3 ⟩ estimatedby the proposed and standard UQ methods. For ⟨x′ y′ ⟩ both
UQ techniques lead to similar PDFs which include the empirical uncertainty estimate. However, for ⟨x′3 ⟩, only the proposed
algorithm can accurately estimate the sample uncertainty. The
fact that the standard method is only accurate for one of these
considered moments is problem-dependent, but the main conclusion is that the proposed algorithm is reliable for estimating
uncertainty in various moments using single realization data.
Aside from potentially inaccurate estimates, the standard
UQ method cannotp
be applied to many compound statisticalterms, for instance ⟨x2 ⟩⟨y⟩. There is no way to evaluate thisterm by applying p
only one SME unless the term is replacedby, for instance, ⟨ x2 y⟩, which is obviously a different statistical term. In contrast, the proposed
method requires nuq
pmerically evaluating ⟨x2 ⟩⟨y⟩ ≈ (µx2 + εµ̂x2 )(µy + εµ̂y ) for

where the associated noise samples are chosen to be Gaussianand correlated:[εx , εy ]T ∼ N (0, Cε ) ,a23 ρε a3 b3.
2
ρε a3 b3 b3

Cε =To ensure the stationarity condition, the model parameters a1 ,
a2 , b1 , and b2 should be chosen such that the eigenvalues of
the coefficient matrix in Eq. (7) lie within the unit circle. The
experiments to produce Figure 1 are carried out using a0 to a2
equal to 0.1 , 0.5 , 0.2 and b0 to b2 equal to 0.3 , 0.4 , 0.6. Moreover, a3 = 1, b3 = 2, and ρε = 0.5. Each realization of Eq. (7)starts from a set of random initial samples for x and y taken
jointly from U [0, 1]2 , and a total of 2n samples are generated
where the first n samples are discarded for burn-in. Combining
moments of x and y, any arbitrary form of statistical terms canbe created for which the accuracy of the proposed UQ algorithm can be examined.
First, the sample-estimated uncertainty of the SME of thefirst-order moment ⟨x⟩ is considered. Figure 1(top) shows the
PDF (probability density function) of σ̂ (µ̂x ) estimated by the
ARM-based uncertainty estimator which uses Eq. (2). The
PDF contains the empirical σe of Eq. (6), and since the modeof the PDF is close to σe , the consistency of the ARM-based
uncertainty estimator for ⟨x⟩ is confirmed.
Next, estimation of uncertainty in the pcentral mo⟨x2 ⟩⟨y⟩, isments ⟨x′ y′ ⟩ and ⟨x′3 ⟩ and a nonlinear term
considered. For the first two terms, we can compare theperformance of the proposed UQ method with the standard
method. In the standard method any of the uncertainty estimation techniques in the previous section is applied to the samples of x′ y′ = (x − µ̂x )(y − µ̂y ) and x′3 = (x − µ̂x )3 . This way,
the uncertainty in µ̂x and µ̂y is ignored. In contrast, the proposed UQ algorithm requires decomposing the SME of ⟨x′ y′ ⟩as ⟨x′ y′ ⟩ ≈ µxy − µx µy + εµ̂x′ y′ , where the uncertainty term is,a sufficiently large number of MC samples. This is straightforward, after accurately estimating cov(µ̂x2 , µ̂y ). For this term,
the bottom plot in Figure 1 confirms the accuracy of the uncertainty estimated by the proposed method.

APPLICATION TO TURBULENT FLOWS
The proposed UQ algorithm is general and can be applied
to any turbulence statistics. Here, the application to the timeseries obtained from two canonical wall-bounded turbulent
flows, the flow over a periodic hill and turbulent pipe flow, isconsidered. For these flows, there are respectively one and two
spatially homogeneous directions at which periodic conditions
are applied at each time step during simulation. Representing spatially-averaged quantities by an overbar, a central moment, for instance ⟨u′ v′ ⟩, can be approximately decomposed as⟨uv − u v⟩ = ⟨uv⟩ − ⟨u v⟩ which is not, in general, equal to thestandard decomposition ⟨uv⟩ − ⟨u⟩⟨v⟩, noting that the covariance (in time) between u and v is not generally zero. This simplification leads to reducing the number of primitive SMEs andhence the size of Σ when computing the uncertainty in central
moments. For the considered flow cases, this approximation isvalid because of having relatively long periodic directions withlarge number of grid points which result in low uncertainty inspatially-averaged quantities. Note that the standard decomposition can also be treated with the proposed approach which
is however not shown here due to brevity.
In all cases discussed below, the uncertainty in the SMEsof the primitive variables is estimated by an ARM method, see
Oliver et al. (2014); Xavier et al. (2022). The steps to estimate
uncertainty are the same as those explained in the previous sections and applied to the model problem. A main characteristicof turbulence signals is, however, that the autocorrelations aswell as cross-correlations decay much slower than the model
problem considered above. This reflects the long history effects in turbulent flows. Moreover, such correlations dependon the variables and spatial location. Thus, careful monitoringof the modeling of the autocorrelations and cross-covariances
when using Eqs. (2) and (4), is essential.

εµ̂x′ y′ ∼ (N (0, σµ̂2xy ) − N (0, σµ̂2x )N (0, σµ̂2y )− µx N (0, σµ̂2y ) − µy N (0, σµ̂2x )) .(8)µ = {µ̂x , µ̂y , µ̂xy } in Eq. (3). Having estimated elTherefore, µ̂
ements of Σ, expression (8) can be evaluated at any Monte
Carlo sample taken from the multivariate Gaussian distribution (3). The result of subtraction/addition of correlated Gaussian distributions is still Gaussian. But multiplication of correlated Gaussian distributions as for instance in Eq. (8) is not
Gaussian. Therefore, the uncertainty in the sample estimationof a compound term such as ⟨x′ y′ ⟩ can, in general, be nonGaussian. Although it may happen that the nonlinear termin Eq. (8) is much smaller than the others and hence the resulting distribution of εµ̂x′ y′ remains Gaussian. This meansthat the distribution of the uncertainty of a compound statistical term depends on the problem at hand and has to be3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022
800

ARM method
Empirical estimate

PDF

600
400
200
0

0.022

0.024

0.026
[]

PDF

150

0.028

0.030

Proposed method
Standard method
Empirical estimate

100
50
0.07

0.08
[]

0.10

Proposed method
Standard method
Empirical estimate

30
PDF

0.09

y/h

0

3.0
2.5
2.0
1.5
1.0
0.5
0.0

20

x/h = 2x/h = 4x/h = 6

0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2< u > /Ub

10
0

0.1

0.2
[]

Proposed method
Empirical estimate

300
PDF

0.3

200
100
0

0.045

0.050

[]

0.055

0.060

y/h

Figure 1. From top to bottom: PDFs
p of sample-estimateduncertainty of ⟨x⟩, ⟨x′ y′ ⟩, ⟨x′3 ⟩, and ⟨x2 ⟩⟨y⟩ obtained fromthe proposed and standard methods compared to the associatedempirical uncertainty by Eq. (6) (vertical line). The uncertainties are computed based on Ne = 104 repetitions of system (7)with sample size n = 50000 per realization. The uncertainty inthe primitive SMEs (diagonal elements of Σ ) are estimated by
the ARM-based estimator. The off-diagonal elements of Σ areestimated from Eq. (4) using modeled γi j (±k).

3.0
2.5
2.0
1.5
1.0
0.5
0.0

x/h = 2x/h = 4x/h = 60.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
< u0u0 > /Ub2

Figure 2. Contours and profiles of ⟨u⟩/Ub and ⟨u′ u′ ⟩/Ub2 of
the flow over a periodic hill at Reb = 2800. The profiles aretaken at x/h = 2, 4, 6 (vertical lines in the contour plots) where
solid lines represent the expected value and shaded areas show
95% confidence intervals (CIs). Note that the CIs of ⟨u⟩/Ub
are multiplied by 3 for better visibility. The profiles with thedotted lines represent the DNS of Breuer et al. (2009).

Flow Over a Periodic Hill
The periodic-hill case is chosen because of its complexity owing to the inherent physical uncertainty in the separated
region behind the first hill. The bulk Reynolds number is
Reb = Ub h/ν = 2800 where Ub is the bulk velocity, h is theheight of the hill, and ν is the kinematic viscosity. The simulation is performed using the high-order spectral-element solver,
Nek5000 (Fischer et al., 2008), based on 9th -order polynomial basis approximations per element in each of the spatial
directions. Figure 2 shows the 95% confidence intervals forthe time-averaged streamwise velocity ⟨u⟩ and the Reynoldsstress component ⟨u′ u′ ⟩ at three streamwise locations, x/h = 2,
4, and 6, for averaging over 44 flow-through times. Reference
data by Breuer et al. (2009) are also included. Interestingly,
the region with the separated flow i.e. x/h = 2 has a lower uncertainty in ⟨u⟩ compared to the downstream locations. Thisis because the autoregressive model can exactly model finite
temporal correlations embedded within the time series of thelarge separation bubble occurring behind the first hill. The
uncertainty in the high-order statistics also reflects this observation. At each of the three streamwise locations, the largest
uncertainty in ⟨u′ u′ ⟩ is observed at the peak near the bottom
wall. This location is where the shear flow originating fromthe top of the first hill is developed.

Turbulent Pipe Flow
The second case we considered is turbulent flow in asmooth straight pipe of circular cross-section. The direct numerical simulation (DNS) of the Navier–Stokes equations isconducted using the open-source code “OPENPIPEFLOW”,

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022
Table 1. Grid resolution parameters for DNS of turbulent
pipe flow at Reb = 5300.

Nr × Nθ × Nz∆r+∆(Rθ )+∆z+

192 × 256 × 1024

(0.02-2)

4.4

5.5

for i, j = 1, 2, 3, where the successive terms on the right-handside of this equation denote production (P), turbulent diffusion (TD), pressure diffusion (PD), pressure strain (PS), viscous diffusion (VD), and dissipation (D) of ⟨u′i u′j ⟩. The definition of these terms can be found, for instance in Moser &Moin (1984). Estimating the uncertainty in the budget terms
is straightforward after decomposing each term into the moments to which the UQ estimators described above can be applied. The only point is about terms having spatial derivativesof moments. For the specific case of turbulent pipe flow, where
samples are averaged over z and θ before time-averaging, estimation of uncertainty in d⟨M⟩/dr is straightforward, noting
the differentiation in r can commute with the time-averaging.
Therefore, it is enough to collect time samples of elements M
during the simulation. This method has the benefit of retaining
the numerical accuracy of the spatial derivatives that is important especially in higher-order flow solvers. An alternative
would be extending Σ to account for covariances in space (associated to the points in a stencil used for computing spatial
derivatives) between primitive SMEs. Figure 5 shows the uncertainty in various terms in Eq. (11) for turbulent kinetic energy (TKE) and ⟨u′r u′z ⟩. For the latter, the pressure strain anddiffusion exhibit larger uncertainty compared to other budgets
especially at small y+ . At these locations, the expected valueof the mentioned terms also deviate more from the DNS of El
Khoury et al. (2013).

which is a specific solver developed for pipe flows using
Fourier-finite-difference methods (Willis, 2017). The bulk
Reynolds number Reb = 2Ub R/ν is chosen to be 5300, andthe corresponding friction Reynolds number (Reτ ≡ ⟨uτ ⟩R/νwith ⟨uτ ⟩ = (⟨τ w ⟩/ρ)1/2 the friction velocity) is 181.95,
where ρ, ν, and R are the fluid density, kinematic viscosity,and pipe radius, respectively. In what follows, the radial, azimuthal, and axial directions are represented by r, θ , and z,
respectively, and · denotes averaging over θ and z. The pipe
length Lz /R is 10π and the details of the spatial resolutions are
summarized in Table 1. Note that this low Reynolds number ischosen here only for illustration purposes, and we have alreadyobtained data at higher Reynolds numbers up to Reτ = 5200
(Yao et al., 2021), where the current UQ analyses are applied
as well.
Figure 3 shows the the 95% confidence intervals for the
first- and second-order velocity moments. All the quantities
are inner scaled using the time-averaged wall friction velocity ⟨uτ ⟩ and the viscous length scale ν/⟨uτ ⟩. One of the strong
points of the proposed UQ approach is the possibility of accurately including the uncertainty in moments of wall-friction
velocity uτ when computing uncertainty in inner-scaled quantities. For instance, an inner-scaled Reynolds-stress component reads as,⟨u′i u′j ⟩+ =µu′i u′j + εµ̂u′ u′i j≈.
⟨uτ ⟩2(µuτ + εµ̂uτ )2⟨u′i u′j ⟩CONCLUSIONS
A framework is introduced for accurately estimating timeaverage uncertainty in sample-estimated statistics of turbulent
flow simulations. A set of guidelines connected to the turbulence physics is introduced to adjust the hyperparameterswhich control the accuracy of various uncertainty estimators.
For quantifying uncertainties in high-order and general turbulence statistics, a novel UQ algorithm is proposed which takes
into account autocorrelation and cross-covariances in turbulence time-series and associated sample-mean estimators. The
approach shows an excellent accuracy when its sample uncertainty estimations are validated against the ensemble empiricaluncertainties for a model problem. The algorithm is demonstrated to be general and accurate for various types of statistics, some of which cannot be dealt with using standard uncertainty estimation methods. To illustrate the application in turbulent flow simulations, DNS of two canonical wall-bounded
turbulent flows are considered. A next development step will
include transferring the workflow to the in-situ framework of
Gscheidle et al. (2022). This will result in a set of efficient andreliable tools for monitoring convergence of various statistics
in large-scale simulations of turbulent flows independent of theflow solver.

(10)µ = {µ̂u′i u′j , µ̂uτ } and associated covariance matrix ΣHaving µ̂
estimated, the right-hand-side of the above expression canbe evaluated. As expected, in Figure 3 the uncertainty for⟨uz ⟩+ is large near the pipe center. In addition, it turns outthat compared to other Reynolds-stress components, the uncertainty is relatively large for ⟨u′z u′z ⟩+ , particularly near the
peak (y+ ≈ 15). Also note that the uncertainty varies between
the Reynolds-stress components.
In a similar way, the uncertainty in higher-order velocity
moments and any compound term can be estimated. Figure 4
illustrates the uncertainty in the profiles of skewness(uz ) =′2 3/2 and kurtosis(u ) = ⟨u′ 4 ⟩/⟨u′ 2 ⟩2 . As expected,⟨u′3zz ⟩/⟨uz ⟩zz
by increasing the order of a central moment, the uncertainty
can increase, compare the uncertainties in kurtosis(uz ) and the
Reynolds-stress components. In fact, for y+ ≲ 3, the uncertainty for kurtosis(uz ) is considerable. For both terms in Figure 4, an increase in the uncertainty is captured at y+ ≈ 140.
The contribution of each primitive moment to the uncertainty
of a compound term can be quantified in the future using elaborate UQ techniques.
The uncertainty in various budget terms in the transport
equations of the Reynolds stresses and turbulent kinetic energy can also be computed. The transport equation for the
Reynolds-stress component ⟨u′i u′j ⟩ in a periodic pipe can be
written as,D⟨u′i u′j ⟩Dt

ACKNOWLEDGMENTS
This work has been supported by: i) EXCELLERAT
project which has received funding from the European Union’s
Horizon 2020 research and innovation programme under grant
agreement No 823691 (SR, PS), and ii) TTU Distinguished
Chair funding (JY, FH). The flow simulations were performedon the resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC (KTH) and Frontera at
Texas Advanced Computing Center.

REFERENCES
Benedict, L. H. & Gould, R. D. 1996 Towards better uncertainty estimates for turbulence statistics. Exp. Fluids 22 (2),
129–136.
Breuer, M., Peller, N., Rapp, Ch. & Manhart, M. 2009 Flowover periodic hills – numerical and experimental study in= Pi j + T Di j + PDi j + PSi j +V Di j + Di j , (11)5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19–22, 2022

uz +

20

uz+rmsur+rmsu +rmsu0 r u0 z +
0

2.5

15

2.0

10

1.5

0
0

1.0

5

0.5

0
10 1

100

y+

101

0.0
10 1

102

100

101

y+

102

Figure 3. Profiles of inner-scaled (left) mean axial velocity and (right) Reynolds-stress components of turbulent pipe at Reb = 5300.
Solid lines represent the expected values and shaded areas show 95% CIs due to time-average uncertainties. Note that the CIs are
multiplied by 5 for better visibility. The uncertainty in ⟨uτ ⟩ is ignored for y+ = ⟨uτ ⟩(R − r)/ν.
1.00

4.5
4.0

0.50

Kurtosis(uz)Skewness(uz)

0.75
0.25
0.00
0.25
10 1

0.2
0.1
0.0

100

y+

101

10 1

102

100

y+

101

102

Profiles of (left) skewness and (right) kurtosis of the streamwise velocity uz of the turbulent pipe flow at Reb = 5300.

P
D
PS
PD
TD
VD

0.075
0.050
0.025
0.000
0.025
0.050
0.075

0.1
10 1

3.0
2.5

0.50

Figure 4.

3.5

100

y+

101

10 1

102

P
D
PS
PD
TD
VD
100

y+

101

102

Figure 5. Inner-scaled budget of (left) turbulent kinetic energy and (right) ⟨u′r u′z ⟩. The 95% CIs of the turbulent kinetic energy are
multiplied by 5 for better visibility. The dashed lines represent the DNS of El Khoury et al. (2013).

a wide range of Reynolds numbers. Comput. Fluids 38 (2),433–457.
El Khoury, G. K., Schlatter, P., Noorani, A., Fischer, P. F.,
Brethouwer, G. & Johansson, A. V. 2013 Direct numerical simulation of turbulent pipe flow at moderately highReynolds numbers. Flow Turbul. Combust. 91 (3), 475–495.
Fischer, P. F., Lottes, J. W. & Kerkemeier, S. G. 2008
NEK5000: Open source spectral element CFD solver.
Available at: http://nek5000.mcs.anl.gov.
Gscheidle, C., Rezaeiravesh, S., Garcke, J. & Schlatter, P. 2022
In-situ estimation of time-averaging uncertainties in turbulent flow simulations. In Workshop MMCP, HPC Asia 2022,
Jan 12-14, Kobe, Japan.
Hoyas, S. & Jiménez, J. 2008 Reynolds number effects on the
Reynolds-stress budgets in turbulent channels. Phys. Fluids
20 (10), 101511.
Lee, M. K. & Moser, R. D. 2015 Direct numerical simulationof turbulent channel flow up to Reτ ≈ 5200. J. Fluid Mech.
774, 395–415.
Moser, R. D. & Moin, P. 1984 Direct numerical simulation of
curved turbulent channel flow. NASA Technical Memorandum 85974.
Oliver, T. A., Malaya, N., Ulerich, R. & Moser, R. D. 2014

Estimating uncertainties in statistics computed from direct
numerical simulation. Phys. Fluids 26 (3), 035101.
Rezaeiravesh, S., Vinuesa, R. & Schlatter, P. 2021 UQit: A
Python package for UQ in CFD. JOSS 6 (60), 2871.
Rezaeiravesh, S., Vinuesa, R. & Schlatter, P. 2022 An
uncertainty-quantification framework for assessing accuracy, sensitivity, and robustness in computational fluid dynamics. J. Comput. Sci. 62, 101688.
Russo, S. & Luchini, P. 2017 A fast algorithm for the estimation of statistical error in DNS (or experimental) time
averages. J. Comput. Phys. 347, 328 – 340.
Smith, R. C. 2013 Uncertainty Quantification: Theory, Implementation, and Applications. Philadelphia: SIAM.
Willis, A. P. 2017 The Openpipeflow Navier–Stokes solver.
SoftwareX 6, 124–127.
Xavier, D., Rezaeiravesh, S., Vinuesa, R. & Schlatter, P. 2022
Reliable quantification of time-averaging uncertainties inturbulent flow simulations. To be submitted .
Yao, J., Schlatter, P., Rezaeiravesh, S. & Hussain, F. 2021
Direct numerical simulation of turbulent pipe flow up toReτ =5200. Bulletin of the American Physical Society .

6