W N Edeling
email: wouteredeling@gmail.com

• G Iaccarino

• P Cinnella

Data-Free and Data-Driven RANS Predictions with Quantified Uncertainty

Keywords: RANS modeling, Uncertainty quantification, Bayesian inference, Return to eddy viscosity, Lag model

For the purpose of estimating the epistemic model-form uncertainty in Reynolds-Averaged Navier-Stokes closures, we propose two transport equations to locally perturb the Reynolds stress tensor of a given baseline eddy-viscosity model. The spatial structure of the perturbations is determined by the proposed transport equations, and thus does not have to be inferred from full-field reference data. Depending on a small number of model parameters and the local flow conditions, a 'return to eddy viscosity' is described, and the underlying baseline state can be recovered. In order to make predictions with quantified uncertainty, we identify two separate methods, i.e. a data-free and data-driven approach. In the former no reference data is required and computationally inexpensive intervals are computed. When reference data is available, Bayesian inference can be applied to obtained informed distributions of the model parameters and simulation output.

Introduction

Despite the increasing availability of high-performance computational resources, Reynolds-Averaged Navier-Stokes (RANS) closure models are projected to remain a widely used option for the prediction of turbulent flows [START_REF] Slotnick | Cfd vision 2030 study: a path to revolutionary computational aerosciences[END_REF]. However, it is well known that the resulting predictions are potentially sensitive to parametric and model-form uncertainty. The former concerns imperfectly known closure coefficients, while the latter deals with the uncertainty due to assumptions made in the mathematical formulation of the model itself. A user of RANS models should therefore be aware of flow-dependent performance, and expect to obtain predictions which can be trustworthy for certain regions of a flow domain, yet erroneous elsewhere.

Over the years various attempts have been made to quantify the uncertainty in RANS closures. One might consider the closure coefficients as random variables, and perform a Bayesian inference procedure to learn their posterior distribution from reference data [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF][START_REF] Edeling | Bayesian estimates of parameter variability in the k-ε turbulence model[END_REF]. To also obtain a measure for the model-form uncertainty, it is possible to combine multiple closure models in a single predictive framework using Bayesian model averaging [START_REF] Hoeting | Bayesian model averaging: a tutorial[END_REF]. Moreover, in traditional turbulence modeling practices, closure models are calibrated on a small set of canonical flow cases. The assumption is then made that these calibrated models can be applied to flows different from those in the training set. To relax this assumption, one can also average over different training scenarios as well in an attempt to broaden the predictive capability, see [START_REF] Edeling | Predictive rans simulations via bayesian model-scenario averaging[END_REF]. The downside of an ensemble approach is if all models in the chosen set are linear eddy-viscosity models, the underlying Boussinesq hypothesis is not challenged, thus resulting in a potential bias in the predictions.

We will build upon the work of [START_REF] Emory | Modeling of structural uncertainties in reynolds-averaged navierstokes closures[END_REF][START_REF] Gorlé | A framework for epistemic uncertainty quantification of turbulent scalar flux models for reynolds-averaged navier-stokes simulations[END_REF], in which perturbations in the eigenvalues of the Boussinesq anisotropy tensor are introduced. The perturbations were either homogeneous, or informed by ad-hoc marker functions. To make the perturbations representative of local flow features, we introduce two additional transport equations for linear combinations of these aforementioned eigenvalues. This limits the amount of information that must be extracted from reference data, since the spatial structure of the perturbations is now embedded in the mathematical form of the transport equations. The use of transport equations to prescribe the local perturbations also tends to impose a regularity in the solution, avoiding oscillatory Reynolds stress components. Our model form is similar to the lag model of [START_REF] Olsen | The lag model, a turbulence model for non equilibrium flows[END_REF], which can be derived from an argument involving a correction to the Boussinesq hypothesis. The location, magnitude and direction of the eigenvalue perturbations are now governed by the two model transport equations. The general behavior of our discrepancy model is determined by two coefficients, resulting in a low-dimensional Uncertainty Quantification (UQ) problem. Depending on the value of these coefficients and the local flow physics, our transport equations will revert to the underlying baseline Boussinesq state. In a manner similar to return to isotropy models, our transport equations therefore describe a 'return to eddy viscosity'.

The low-dimensional nature of our model allows for two distinct applications in a predictive setting, which we denote the data-free and data-driven predictions. If no reference data is available for the flow case at hand (or for a closely related flow), cheap uncertainty intervals can be computed. When we do have access to data, a Bayesian inference procedure can be applied in order to obtain informed probability distributions describing the uncertainty in the coefficients and Quantity-of-Interest (QoI). Although this data-driven approach is more expensive, accurate surrogate models can be readily constructed with O(100) RANS samples. Markov-chain Monte-Carlo (McMC) methods, or other common sampling techniques, are then available to obtain posterior distributions of the model coefficients. These distributions are easily interpreted from a physical point of view, as the values of the model coefficients can be directly related to a state of anisotropy. As validation cases we will consider a backward facing step and a subsonic jet, for which it is known that Boussinesq-type models struggle to accurately capture all relevant flow physics.

This paper is organized as follows. First, we will very briefly describe the RANS modeling framework. Sections 3 and 4 then deal with the eigenvalue perturbations, and describe various aspects (e.g. realizability and effect of model coefficients) of the proposed model. Next, the data-free and data-driven methods are outlined, followed by the obtained results in Section 6. Finally, we give our conclusions in Section 7.

Reynolds-Averaged Navier-Stokes Equations

The well-known Reynolds-Averaged Navier-Stokes equations [START_REF] Pope | Turbulent flows[END_REF] are given by

∂ρ ∂t + ∂ρ U i ∂x i = 0, ∂ρ U i ∂t + ∂ρ U i U j ∂x j = ∂ ∂x j σ ij -p δ ij -ρ u i u j . ( 1
)
Here, ρ is the density, U i is the mean velocity in x i direction and μ is the viscosity. The three terms on the right-hand side represent the viscous stress tensor, the isotropic stress due to the average pressure field p , and the apparent stress u i u j arising from the fluctuations in the velocity. To close (1), one must model the Reynolds stress tensor u i u j in terms of mean-field quantities. A common RANS closure is the Boussinesq hypothesis, i.e.

u i u j ≈ u i u j (bl) = 2 3 kδ ij -2ν T S ij . ( 2
)
Here, k is the turbulent kinetic energy, δ ij the Kronecker delta, and S ij is the mean strain-rate tensor, defined as

S ij := 1 2 ∂ U i /∂x j + ∂ U j /∂x i -1 3 ∂ U k /∂x k δ ij in com- pressible flows.
The superscript (bl) stands for baseline, and denotes a quantity computed using the Boussinesq hypothesis.

The scalar quantity ν T is the eddy-viscosity, which must be computed by means of some chosen turbulence model, e.g. the kε, kω or Spalart-Allmaras models, see [START_REF] Wilcox | Turbulence modeling for cfd[END_REF] for an overview. Besides the approximate nature of these transport equations, several closure coefficients and corresponding uncertainties are also present, ultimately leading to an uncertain ν T [START_REF] Edeling | Bayesian estimates of parameter variability in the k-ε turbulence model[END_REF][START_REF] Edeling | Predictive rans simulations via bayesian model-scenario averaging[END_REF]. Instead of focusing on one model-specific form of ν T , we will inject uncertainty directly into u i u j , by adding a tensorial discrepancy model D ij to the right-hand side of Eq. 2, i.e. we set

u i u j = u i u j (bl) + D ij . (3)
We therefore leave any chosen baseline turbulence model and its coefficients unperturbed, and instead focus on estimating D ij . Specifically, D ij will be specified (see Section 4) based on perturbations in the eigenvalues of the Reynolds-stress anisotropy tensor b ij . It is important to note that we do not perturb the corresponding eigenvectors. As a result, we therefore enforce that u i u j is aligned with u i u j (bl) throughout the flow domain.

To perform the simulations the commercial code ANSYS Fluent v16.1 was used. Our model was implemented into Fluent by means of a set of User-Defined Functions (UDFs). These UDFs are available from the authors upon request.

Reynolds Stress Anisotropy and the Barycentric Map

Let us define the normalized anisotropy tensor b ij as

b ij := u i u j 2k - 1 3 δ ij ; - 1 3 ≤ b αα ≤ 2 3 ; - 1 2 ≤ b αβ ≤ 1 2 ; (4)
which is a symmetric, deviatoric (zero trace) tensor. Greek subscripts are excluded from the summation convention. Banerjee et al. express the anisotropy tensor as a convex combination of three basis tensors b1c , b2c , b3c , i.e., bij = C 1c b1c + C 2c b2c + C 3c b3c .

(5)

These basis tensors represent the three limiting states of componentality (relative strengths of components in u i u j ), i.e., they represent one-, two-and three-component turbulence. A circumflex (• ) represents the tensor (•) expressed in principal axes, e.g. bij = diag (λ 1 , λ 2 , λ 3 ) , [START_REF] Emory | Modeling of structural uncertainties in reynolds-averaged navierstokes closures[END_REF] where

λ 1 ≥ λ 2 ≥ λ 3 are the eigenvalues of b ij .
The basis tensors of Eq. 5 are given by -One-component limiting state (1C): b1c := diag 2 3 , -1 3 , -1 3 . One component (and one eigenvalue) of u i u j is nonzero. Note that a zero eigenvalue of u i u j leads to an eigenvalue of -1 3 for bij . -Two-component axi-symmetric limiting state (2C): b2c := 1 6 , 1 6 , -1 3 . Two components of u i u j are nonzero.

-Three-component limiting state (3C): b3c := (0, 0, 0). All three components of u i u j are equal, where the limiting basis tensor corresponds to isotropic turbulence.

From Eq. 5 it is clear that the coefficients C 1c , C 2c and C 3c measure how close bij is to any of the three limiting states. Since Eq. 5 is required to be a convex combination, we have by definition:

C 1c + C 2c + C 3c = 1 and C 1c ≥ 0 C 2c ≥ 0 C 3c ≥ 0. ( 7
)
C 1c = λ 1 -λ 2 λ 1 = C 1c + C 2c 2 + C 3c 3 - 1 3
C 2c = 2 (λ 2 -λ 3 ) ⇔ λ 2 = C 2c 2 + C 3c 3 - 1 3
C 3c = 3λ 3 + 1 λ 3 = C 3c 3 - 1 3 (8)
Due to Eq. 7, again only two coefficients are needed to quantify the state of anisotropy.

To visualize the nature of the anisotropy implied by these coefficients, a barycentric map can be defined as ( [START_REF] Banerjee | Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches[END_REF])

x b = C 1c x 1c + C 2c x 2c + C 3c x 3c , y b = C 1c y 1c + C 2c y 2c + C 3c y 3c . ( 9
)
Here, (x 1c , y 1c ), (x 2c , y 2c ) and (x 3c , y 3c ) are the three corner points corresponding to the limiting states of componentality. They can be chosen arbitrarily, but are commonly set to the corner points of an equilateral triangle. In Fig. 1 the barycentric map is depicted along with the variation of the coefficients along its edges. Each point x ∈ R N in the spatial flow domain has its own coefficients and thus can be mapped to a location (x b , y b ) in the barycentric map using Eq. 9. Since Eq. 5 is a convex combination of limiting states, all possible realizable states of u i u j are therefore contained within the borders of the barycentric map.

Anisotropy tensor decomposition

The eigen decomposition of the anisotropy tensor reads

b ij = v ik kl v jl , ( 10
)
where v ij are the eigenvectors and ij the diagonal matrix of eigenvalues, ordered as λ 1 ≥ λ 2 ≥ λ 3 . Hence, from Eq. 4 we see that the Reynolds stress tensor can be written as

u i u j = 2k v ik kl v jl + 1 3 δ ij (11)
Here, the turbulent kinetic energy, eigenvectors and eigenvalues represent the magnitude, orientation and shape of the Reynolds stress respectively. As noted, we keep the alignment (v ij = v (bl) ij ), and perturb the eigenvalues ij , such that

u i u j = u i u j (bl) + D ij = 2k v (bl) ik (bl) kl v (bl) jl + 1 3 δ ij + 2k v (bl) ik kl v (bl) jl , ( 12
)
where ij is the diagonal tensor containing the perturbations in the eigenvalues. To be clear, what we denote as uncertainty here is thus related to the uncertain shape of u i u j . The Boussinesq hypothesis restricts this shape to a relative small subset of all possible tensor shapes, see Fig. 2. By including ij , we aim to break away from this particular restriction of the Boussinesq model. It should also be noted that what matters to the QoIs might not be reflected in the Reynolds stresses. Uncertainties in boundary-or inflow conditions can also significantly impact the solution, and may act in locations where the epistemic model form error is relatively small.

Emory and Iaccarino [START_REF] Emory | Modeling of structural uncertainties in reynolds-averaged navierstokes closures[END_REF][START_REF] Emory | Modeling structural uncertainties in reynolds-averaged computations of shock/boundary layer interactions[END_REF] also used Eq. 11 to perturb the baseline eigenvalues. The perturbation direction and magnitude was either specified homogeneously or allowed to vary spatially through ad hoc sensors based on wall distance [START_REF] Emory | Modeling structural uncertainties in reynolds-averaged computations of shock/boundary layer interactions[END_REF] or streamline curvature a b Fig. 2 Tensor shapes of u i u j /2k visualized using super-quadratic tensor glyphs [START_REF] Kindlmann | Superquadric tensor glyphs[END_REF], figure generated with Mayavi [START_REF] Ramachandran | Mayavi: 3d visualization of scientific data[END_REF]. The 1C, 2C and 3C limiting shapes are clearly visible on the right [START_REF] Gorlé | Epistemic uncertainty quantification for rans modeling of the flow over a wavy wall[END_REF]. Instead, we will propose transport equations to specify the magnitude, direction and spatial distribution of the eigenvalue perturbations.

Return to eddy viscosity model

Since the Boussinesq hypothesis incorporates a direct proportionality between u i u j and the mean strain-rate tensor S ij (see Eq. 2), the turbulence reacts directly to changes in the mean flow. Strictly speaking, this is only a valid approximation under equilibrium conditions, for example when the time scale of the turbulence is much smaller than that of the mean flow. Our perturbation framework is based upon the lag-type models, which modulate the response through simple linear PDE's on ν t [START_REF] Olsen | The lag model, a turbulence model for non equilibrium flows[END_REF] or u i u j [START_REF] Lillard | The lagrst model: a turbulence model for non-equilibrium flows[END_REF], effectively slowing down the response when appropriate.

We prescribe the following lag equations for the coefficients of the barycentric map

DC 1c Dt = a 1c 1 τ C (bl) 1c -C 1c , DC 2c Dt = a 2c 1 τ C (bl) 2c -C 2c . ( 13
)
It is important to note that this model is first and foremost an inspiration to construct spatially dependent eigenvalue perturbations suitable for the purpose of uncertainty quantification, rather than a formal framework to build a new turbulence closure. The quantity τ -1 is an inverse turbulent time scale (modeled as ε/k or ω) and

C (bl) 1c , C (bl)
2c are the coefficients computed using the eigenvalues of the Boussinesq anisotropy tensor, i.e. from b (bl) ij = -

ν T k S ij . ( 14
)
The right-hand side of Eq. 13 is similar to that of a linear return to isotropy model, see e.g. Rotta's model [START_REF] Rotta | Statistische theorie nichthomogener turbulenz[END_REF]. However, we consider a baseline eddy-viscosity model which can be erroneous at certain locations, but which is able to yield accurate predictions outside these regions. The transport equations for C 1c and C 2c are therefore meant to describe a 'return to eddy viscosity' when the Boussinesq hypothesis ( 2) is expected to be valid. Along a streamline, a perturbed barycentric coefficient will go to its baseline value with a time scale τ . The rate of return for each C αc is independently controlled by two coefficients a αc . This enables [START_REF] Ramachandran | Mayavi: 3d visualization of scientific data[END_REF] to explore a larger area of the barycentric map, making the lag concept amenable for uncertainty quantification, see the discussion in Section 4.3.

We define the boundary conditions following an approach similar to the original lag model [START_REF] Olsen | The lag model, a turbulence model for non equilibrium flows[END_REF]. At a solid smooth wall, we set the values of C 1c and C 2c equal to the baseline values C (bl) 1c and C (bl)

2c . If the turbulence model equations are integrated down to the wall using damping functions, this essentially sets the values equal to zero. However, if wall functions are used, non zero values of C 1c and C 2c are obtained. At an inflow boundary, it is assumed that the freestream turbulence is isotropic, i.e. u i u j = 2 3 k ∞ δ ij , yielding again C 1c = C 2c = 0. An isotropic inflow condition is equivalent to the option of specifying the turbulence intensity at an inlet, which can be done in many solvers. This is not expected to be universally valid, but this is also not the only option. If one has an informed guess, from higher-fidelity simulations for instance, specifying a profile is also possible. Outflow boundary values are extrapolated from internal cells.

Since the C αc are such simple linear combinations of the eigenvalues (see Eq. 8), an equivalent approach would be to write lag equations for λ α , α = 1, 2. We prefer the use of the barycentric coefficients due to their intuitive properties. All 3 have realizable values in the same range ([0, 1]), and each one measures the distance from the corresponding limiting state of anisotropy, being 1 if at the limit and 0 if completely removed from it.

Finally, it should be noted that lag models have been further developed since the introduction of the lag concept, see for instance the stress-strain lag model of [START_REF] Revell | A stress-strain lag eddy viscosity model for unsteady mean flow[END_REF]. We use the original form proposed by [START_REF] Olsen | The lag model, a turbulence model for non equilibrium flows[END_REF] as a proof of concept for the use of the lag concept as a tool for uncertainty quantification.

Effect of lag coefficients

The original default values of the lag coefficients, as defined by [START_REF] Olsen | The lag model, a turbulence model for non equilibrium flows[END_REF], are

a 1c = a 2c = 0.35. ( 15
)
Let us instead investigate the model behavior for three cases that represent limiting trajectories:

-Case 1 (plane-strain): a 1c = a 2c = 0.35, -Case 2 (axi-symmetric contraction): a 1c = 0, a 2c = 0.35, -Case 3 (axi-symmetric expansion): a 1c = 0.35, a 2c = 0.

The typical resulting trajectories in the barycentric map are shown in Fig. 3. A model will follow the plane strain turbulence line when at least one eigenvalue of b ij is zero [START_REF] Banerjee | Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches[END_REF].

Fig. 3 The trajectories in the barycentric map of the three coefficient cases. The dashed line denotes the plane strain line where λ 2 = 0 Note that the Boussinesq anisotropy tensor [START_REF] Emory | Modeling structural uncertainties in reynolds-averaged computations of shock/boundary layer interactions[END_REF] will yield a zero eigenvalue when one of the S μμ equals zero. So all two-dimensional flow cases (or three-dimensional flows like the fully-developed square duct where S 11 = 0) will yield trajectories along the plane stain line when a Boussinesq model is used, and b(bl) ij = diag(λ 1 , 0, -λ 1 ). When a 1c = a 2c (case 1), both transport equations in Eq. 13 return to the Boussinesq solution at the same rate. As a result, the trajectory in the barycentric map still lies along the plane strain line (see Fig. 3). Along this line the distribution is different than that of the unperturbed model, depending on the values of a 1c and a 2c .

To observe larger deviations from the results of the baseline model, unequal lag coefficients are required. Remember from Fig. 1 that C 1c = 0 along the axi-symmetric contraction boundary. In the second limiting case we enforce a homogeneous solution C 1c = 0 by setting a 1c = 0. This therefore forces the trajectory to follow the axi-symmetric contraction border (see Fig. 3). Similarly, setting a 2c = 0 results in a trajectory along the axi-symmetric expansion border, where C 2c = 0. Setting two nonzero, unequal, values for a 1c and a 2c will result in curved trajectories.

Spatial distribution of the perturbations

The error in any given baseline model will display a complex spatial distribution, which is dependent on the local flow physics [START_REF] Ling | Evaluation of machine learning algorithms for prediction of regions of high reynolds averaged navier stokes uncertainty[END_REF]. The perturbations in the Boussinesq eigenvalues should reflect this localized behavior, and ideally be concentrated in those areas of the flow domain where a failure of the Boussinesq hypothesis might reasonably be expected. Besides determining the direction and the magnitude of the perturbations in λ (bl) α , the lag terms also dictate the spatial distribution. The exact distribution will be dependent upon the values of a 1c and a 2c . As an example, consider the contour plots of

| λ 1 | := |λ 1 -λ (bl) 1 | and | λ 2 | := |λ 2 -λ (bl)
2 |, for an axisymmetric jet with an exit Mach number of approximately 0.51 [START_REF]Axisymmetric subsonic jet[END_REF], and a backward facing step with a Reynolds number based on the step height h of Re h = 5100 in Fig. 4. Note that for the jet flow only at the nozzle exit and in the mixing layer do we find significant perturbations. Where the flow is essentially uniform (and in the potential core), no perturbations are made in λ 1 . In the backstep case, the local perturbations are made near and after the separation point. Still, when moving away from the wall, the perturbations in λ 2 approach zero.

a b

Symmetry around a 1c = a 2c

In two-dimensional flows, when the baseline state lies on the plane-strain line, certain QoIs will display symmetry in the stochastic domain around a 1c = a 2c . Specifically, -case (1):

a 1c = α, a 2c = β, -case (2): a 1c = β, a 2c = α,
for any chosen values α, β will cause the so-called eigengaps δ 1 := λ 1λ 2 and δ 2 := λ 2λ 3 to switch between case (1) and (2). The following three statements are then all equivalent δ (1)

1 = δ (2) 2 and δ (1) 2 = δ (2) 1 ,
C (1) 1c =

1 2 C (2)
2c and C

(1)

2c = 2C (2) 1c , λ (1) 1 = -λ (2)
3 , λ

(1)

2 = -λ (2) 2
and λ

(1)

3 = -λ (2) 1 . ( 16
)
Here, δ

(2)

1 denotes the first eigengap of case (2), and likewise for the other quantities. To prove the above, use the relations δ 1 = C 1c and δ 2 = C 2c /2 together with (13) to obtain the following transport equations

Dδ (1) 1 Dt = αω δ (bl) 1 -δ (1) 1
and Dδ

(2) 2 Dt = αω δ (bl) 2 -δ (2) 2 . ( 17
)
Clearly, both transport equations have the same form. In addition, the boundary conditions are equal and along the plane-strain line (where λ (bl) 2 = 0 by definition), both baseline eigengaps are always the same due to b(bl

) ij = diag(λ (bl) 1 , 0, -λ (bl)
1 ). Thus, Eq. 17 essentially displays the same right-hand side twice, implying δ (1)

1 = δ (2)
2 . A similar argument can be formed for the second equality δ (1)

2 = δ (2)
1 of Eq. 16. If Eq. 16 is satisfied, u 1 u 2 is the same for case (1) and [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF], which in turn can cause very similar solutions for other QoI. To see this, denote α 1 ∈ [-π, π] be the minimum angle from the horizontal axis to the eigenvector v 1 , which corresponds to the largest eigenvalue λ 1 . The angle α 3 is likewise defined in relation to v 3 . Since b ij is real and symmetric, any two eigenvectors from different eigenspaces are orthogonal to each other [START_REF] Lay | Linear algebra and its applications[END_REF]. Moreover, any two orthogonal eigenvectors in the x-y plane can be described by v 1 = (cos α 1 , sin α 1 , 0) T and v 3 = (cos α 3 , sin α 3 , 0) T , where α 1 = α 3 ± π/2. The remaining eigenvector must be orthogonal to the x-y plane, and thus v 2 = (0, 0, 1) T is also an eigenvector. Now we can write b ij as

b ij = ⎛ ⎝ cos α 1 0 cos α 3 sin α 1 0 sin α 3 0 1 0 ⎞ ⎠ ⎛ ⎝ λ 1 0 0 0 λ 2 0 0 0 λ 3 ⎞ ⎠ ⎛ ⎝ cos α 1 sin α 1 0 0 0 1 cos α 3 sin α 3 0 ⎞ ⎠ . ( 18
)
From Eq. 18, we can extract the following expression for the Reynolds shear-stress component after some manipulation

u 1 u 2 k = 2 (λ 3 -λ 1 ) sin (2α 3 ) = -C 1c + 1 2 C 2c sin (2α 3 ) . ( 19
)
Provided that we do not perturb the eigenvectors, u 1 u 2 /k is determined by the coefficient C 1c + C 2c /2 only. From the relations in Eq. 16) it is readily determined that this coefficient is the same for both case (1) and [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF]. As a verification exercise we randomly picked values for α and β, and compared the C 1c + C 2c /2 and u 1 u 2 profiles from the backward facing step flow of Fig. 4. We found that the results for case (1) and (2) indeed overlap.

Realizability with an unrealizable baseline model

To obtain physically feasible values for the components of the Reynolds stress tensor, the realizability constraints first defined by [START_REF] Schumann | Realizability of reynolds-stress turbulence models[END_REF] must hold, i.e.

u μ u μ ≥ 0 u μ u ν 2 ≤ u μ u μ u ν u ν det u i u j ≥ 0, μ,ν ∈ {1, 2, 3}. (20)
Weak realizability conditions can be derived for second-order closures [START_REF] Durbin | Realizability of second-moment closure via stochastic analysis[END_REF], and for Eq. 13 as well [START_REF] Edeling | A return to eddy viscosity model for epistemic uq in rans closures[END_REF]. We now investigate weak realizability of our perturbation model on the backstep case, when the chosen baseline model is unrealizable. In this case C (bl) 1c > 1 3 along the plane-strain line for locations of unrealizable turbulence. This surely can occur, see for instance the trajectory of the baseline kε model for the backward-facing step flow in Fig. 5a, and the corresponding spatial locations of unrealizable turbulence in Fig. 6. We selected the kε model precisely because it shows such a pronounced violation of realizability for this particular flow case, making it a good model to showcase the ability of the lag equations to retrieve a realizable state. For the predictions of Section 5 we used the k -ω SST model. Luckily, we can control the similarity to the Boussinesq solution through the lag coefficients, and limiting the values of a 1c and a 2c can prove sufficient to not follow an unrealizable baseline model across the two-component boundary. Hence, the approach does not only perturb the assumptions in the Boussinesq hypothesis, but also provides an easy means to correct violations of realizability. In Fig. 6 we show how much we have to increase the lag coefficients in order to obtain an unrealizable u i u j for the same backward-facing a b 13) with the default lag coefficient values yields a realizable result, even with the highly unrealizable baseline model. We have to increase a 1c , a 2c by an order of magnitude in order to force the result across the two-component boundary.

Predictions

The final goal of all preceding analysis, is to use Eq. 13 to make predictions with quantified uncertainty. We choose our method of prediction based on the availability of reference data. If no data is available, we will not attempt to extrapolate from other calibrated cases to the prediction scenario at hand. Instead, we compute conservative intervals, based on the method described in Section 5.1. Otherwise, we use the Bayesian inference procedure outlined in Section 5.2.

Data-free prediction

In order to compute the data-free intervals, let us determine the possibility of a priori identifying the lag coefficients a 1c and a 2c that will yield the maximum and minimum bound on the perturbed eigenvalues that can be achieved with model [START_REF] Ramachandran | Mayavi: 3d visualization of scientific data[END_REF]. As noted in Eq. 12, we decompose the anisotropy tensor as

b ij = b (bl) ij + v (bl) ik kl v (bl) jl . ( 21
)
Where ij := diag ( λ 1 , λ 2 , λ 3 ), and λ α is the perturbation such that λ α = λ (bl) α + λ α . Assuming again that the baseline model will follow the plane strain line in the barycentric map, ij can be written as

ij = ⎛ ⎜ ⎜ ⎝ 2 3 C 1c -3 2 C (bl) 1c + 1 6 C 2c 0 0 0 -1 3 C 1c + 1 6 C 2c 0 0 0 -1 3 C 1c -3C (bl) 1c -1 3 C 2c ⎞ ⎟ ⎟ ⎠ . ( 22

) By fixing the baseline coefficient to an arbitrary value C (bl)

1c ∈ [0, 1/3], we can plot the isocontours of the λ α in order to identify the directions of maximum and minimum perturbation. The chosen value of C (bl) 1c will affect the scale but not the shape of the isocontours. The results are shown in Fig. 7. Clearly, λ 1 is maximized in the direction of 1C, and minimized toward 3C. The perturbation λ 2 is maximized at 2C and minimized Fig. 7 Isocontours of eigenvalue perturbations for a given value of C (bl) 1c , which this case is set to zero at 1C. Finally, the maximum of λ 3 is located at 3C, and its minimum is found anywhere along the 2C boundary. Thus, straight lines in the direction of minimum/maximum λ α are the axi-symmetric expansion, contraction and the plane-strain line. Trajectories along these lines can be obtained with a priori known values for a 1c and a 2c (see Fig. 3). The maximum values we use for a αc are discussed in Section 6.1.2.

Data-driven prediction

Bayesian calibration

In a Bayesian approach all sources of uncertainty (parametric, model form and experimental) are condensed into a single expression, i.e Bayes' rule [START_REF] Gelman | Bayesian data analysis[END_REF], given by

p (θ| z) = p(z|θ)p(θ) p(z) . ( 23
)
Here, θ ∈ R C is a vector containing the model coefficients to be inferred, z ∈ R N is a vector of high-fidelity reference data, and p (θ) is the prior distribution.

The likelihood function p(z|θ) is found by specifying a relation between the true process ζ , the corresponding QoI q from the code, and the data z as [START_REF] Kennedy | Bayesian calibration of computer models[END_REF]

z = ζ(x) + e. (24)
Here, x are the spatial coordinates and e represents the uncertainty in z due to measurement noise, commonly modeled as a zero-mean Gaussian e ∼ N (0, ), := diag(λ i ). In [START_REF] Oliver | Bayesian uncertainty quantification applied to rans turbulence models[END_REF], the true process was modeled as

ζ(x) = (1 + δ m (x; γ )) q(x; θ), ( 25
)
but other forms are also possible. The quantity δ m is a multiplicative model error term superposed on the code output, again distributed as a zero-mean Gaussian, and dependent upon a set of hyper parameters γ . Note that without this term we would state that the true process equals our RANS code, which is unlikely to be an accurate assumption. In turn, this yields a Gaussian likelihood such that z|θ ∼ N (q, K). The covariance of the likelihood is K = + K m , where the entries of K m are given by the function

k m (x , x) = q(x ; θ) δ m (x ; γ )δ m (b; γ ) q(x; θ). ( 26
)
Here, δ m (b ; γ )δ m (b; γ ) is the covariance of the model error term, which remains to be specified.

Since Eq. 13 already injects perturbations in the Boussinesq assumption, it can be considered as an internal model error term. However, it is still dependent upon a 1c and a 2c . To prevent the posterior distribution of these coefficients from becoming overly confident, we use a simple spatially uncorrelated model δ m (b ; γ )δ m (b; γ ) := σ 2 δ(xx) when performing the Bayesian inference. Here, σ is a hyper parameter that is calibrated alongside a 1c and a 2c , and δ(xx) is the Dirac delta function. The parameter σ has the effect of broadening the posterior p(θ|z).

The Markov-chain Monte Carlo (McMC) method [START_REF] Chib | Understanding the metropolis-hastings algorithm[END_REF] is often used to draw samples θ i ; i = 1, • • • , N, from the posterior density [START_REF] Durbin | Realizability of second-moment closure via stochastic analysis[END_REF]. The posterior density itself is therefore not inferred directly, and can only be estimated using a kernel density estimation after sufficient samples have been drawn, making McMC an expensive method to use in a RANS context.

A common practice is therefore to construct a polynomial approximation of the forward problem over the support of the prior distribution [START_REF] Marzouk | A stochastic collocation approach to bayesian inference in inverse problems[END_REF]. This approximation then replaces the full RANS code in the likelihood function, yielding a surrogate posterior distribution which can be sampled at minimal cost. We use the Stochastic Collocation (SC) expansion, formed by approximating the forward problem as the sum of N p (multidimensional) Lagrange interpolation polynomials [START_REF] Witteveen | Uncertainty quantification of the transonic flow around the rae 2822 airfoil[END_REF]. This technique is well known and not new, and many alternatives to the SC expansion exist, see e.g. [START_REF] Edeling | Simplex-stochastic collocation method with improved scalability[END_REF][START_REF] Ghanem | Polynomial chaos in stochastic finite elements[END_REF][START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF][START_REF] Dongbin | The wiener-askey polynomial chaos for stochastic differential equations[END_REF]. We therefore refer to [START_REF] Witteveen | Uncertainty quantification of the transonic flow around the rae 2822 airfoil[END_REF] for details.

Constructing a SC surrogate requires N p := m 1 × m 2 evaluations of the RANS code, which in our case corresponds to running the simulation on the collocation points of the tensor grid

a i 1 1c ⊗ a i 2 2c , i 1 = 1, • • • , m 1 , i 2 = 1, • • • , m 2 .
Let qN p denote the corresponding surrogate for some RANS output q. The hierarchical surplus error for the backstep case

E surplus = qN p a i 1 1c , a i 2 2c -qN p-1 a i 1 1c , a i 2 2c 2 , ( 27
)
evaluated at the collocation points of qN p , was O(10 -3 ) for mean velocity profiles and O(10 -5 ) for u i u j components, when using N p = 144 and N p-1 = 121. Note that a SC surrogate is exact ( qN p = q), when evaluated at its collocation points. We selected the surrogate with N p = 144 for use in the likelihood function.

Results

Application to a subsonic jet

Flow case description

This flow a well-known NASA validation case for which the settings and computational grids can be found in [START_REF]Axisymmetric subsonic jet[END_REF]. The jet exit Mach number is approximately 0.51. For computational stability, the jet does not flow into quiescent air, and a small reference Mach number M ref = 0.01 is applied to the farfield. No-slip boundary conditions are applied to the walls of the nozzle, and at the nozzle inlet the total-to reference pressure ratio is p t /p ref = 1.19671, which yields the correct jet exit Mach number. Experimental Particle Image Velocimetry (PIV) data for several quantities, including uncertainty estimates on the data, are available for this flow case from [START_REF] Bridges | The nasa subsonic jet particle image velocimetry (piv) dataset[END_REF]. A selection of experimental results, normalized by the centerline jet-exit velocity U jet and nozzle diameter D jet , are depicted in Fig. 8.

Data-free prediction

Let us start by displaying the baseline results as computed with the k-ω SST closure model. Figure 9 shows the normalized axial velocity U /U jet and Reynolds stress component

u 1 u 1 /U 2 jet .
Note that for the axial velocity, the results agree fairly well with the reference data at x/D jet = 2. However, the accuracy of the baseline prediction deteriorates when we move away from the nozzle exit. In the case of the Reynolds stress component, clear discrepancies can be observed along all profiles.

We now compute the intervals due to the 4 code solves that are predicted to yield the maximum and minimum eigenvalues. i.e. the solves corresponding to the tensor product [a 1c,min , a 1c,max ] ⊗ [a 2c,min , a 2c,max ]. We select the support of the lag coefficients based on two limiting cases. The first case is the straightforward choice of isotropic turbulence, and thus a 1c,min = a 2c,min = 0. Second, for a 1c = a 2c 0 the baseline state is approached. It makes sense that the baseline state should also be included in the interval. We therefore a b select a maximum cutoff value for a αc for which a quantity of interest is sufficiently close to the Boussinesq result. Already for a 1c,max = a 2c,max = 1.0, the velocity components are almost indistinguishable from those obtained with the unperturbed SST model, and only minor differences are observed in the Reynolds stresses, see Fig. 10. Hence we compute the results for the combinations shown in Table 1, and form an envelope from the corresponding outputs. The results for the normalized axial velocity are shown in Fig. 11. Almost everywhere the envelope intersects the region of experimental uncertainty, and where there is no overlap the discrepancy is small. Moreover, as we move away from the nozzle exit (where the baseline model becomes less accurate, see Fig. 9a), the uncertainty bounds increase. We also show the results for the normal Reynolds stress component u 1 u 1 /U 2 jet in Fig. 12. The width of the intervals is more pronounced than in the case of the velocity profile, and the experimental data are completely enveloped. While it is certainly true that the Reynolds stress components have a very large interval, the width of velocity intervals is less pronounced. Thus, one way to interpret the Reynolds stress intervals is that, in the absence of reference data, large intervals can be advantageous in order to capture the model inadequacy in other quantities which are less sensitive to u i u j , at least for this particular flow case.

Data-driven predictions

We will compute the posterior distribution of θ = [a 1c , a 2c ] T , using Eq. 23. As mentioned, the surrogate q will replace the full RANS code in the likelihood function such that we have: z|θ ∼ N ( q, K). Here, K = + K m , and z will be experimental U 1 data from [START_REF] Bridges | The nasa subsonic jet particle image velocimetry (piv) dataset[END_REF]. The reported experimental accuracy is used to specify the diagonal . K m is given by Eq. 26, where δ m (b ; γ )δ m (b; γ ) = σ 2 δ(xx). The hyper parameter σ (with prior σ ∼ U [0, 3]) is calibrated alongside θ. Since in general it will be difficult to a-priori say which coefficients will be most likely, we use a uniform prior p(θ) which is defined on the same range as displayed in Table 1.

We performed a single calibration using the U 1 PIV data at x/D jet = {2, 12, 22, 32, 42}. The resulting posteriors are depicted in Fig. 13. The axisymmetric state is preferred by the data. As a result, p(α 1c |z) resembles a delta distribution centered at 0. The posterior U 1 results are shown in Fig. 14, the variance of which is greatly reduced compared to the interval results of Fig. 11.

Fig. 12 The uncertainty intervals PIV data of u 1 u 1 /U 2 jet along 4 centerline-normal directions. The squares represent the PIV data plus the experimental error, and the shaded area is the RANS uncertainty interval Fig. [START_REF] Ramachandran | Mayavi: 3d visualization of scientific data[END_REF] The posterior distributions of the jet case

Cost overhead

Finally, we use the jet case to briefly give an indication of the cost overhead of our approach. Table 2 displays the relative cost overhead of the four cases of Table 1, compared to the runtime of the jet case using just the baseline SST model. All cases were run for 6000 iterations, and roughly add a 30 % increase to the baseline runtime.

Generally speaking, no convergence issues were encountered for all considered flows and a 1c , a 2c combinations, with the exception of a 1c = a 2c = 0. Depending on the flow case, perturbing all points to the 3C corner can require an increase in the under-relaxation factor of Eq. 13 in order to reach a converged solution. If reducing the under-relaxation is not enough to stabilize the solution, one might employ u i u j (bl) + α( u i u ju i u j (bl) )

with α < 1 as a momentum source term, after applying the divergence theorem. Since this is not true under relaxation (due to the u i u j (bl) terms), it is important to keep γ as close to 1 as possible. 6.2 Application to a backward facing step

Flow case description

This is a well-known backward-facing step case with Re h = 5100, where h represents the step height. DNS data from [START_REF] Le | Direct numerical simulation of turbulent flow over a backward-facing step[END_REF] is available for this flow, as well as experiments performed by [START_REF] Jovic | Backward-facing step measurements at low reynolds number, re h = 5000[END_REF]. At the bottom wall no-slip boundary conditions are applied, and the top boundary a parallel slip condition is prescribed in order to match the DNS setup. The computational domain extends to 40 step heights after the backstep, which is located at x/ h = 0.

Reattachment length

Eddy-viscosity models are known to provide inaccurate predictions of the reattachment location, see e.g. [START_REF] Thangam | Turbulent flow past a backward-facing step-a critical evaluation of twoequation models[END_REF]. The mean calculated reattachment length from the DNS data is l = 6.28h [START_REF] Le | Direct numerical simulation of turbulent flow over a backward-facing step[END_REF]. We now examine the range in l when using our return to eddy viscosity framework. The baseline model is the SST model, and we estimate the reattachment length by the location for which we have a zero skin-friction coefficient, i.e. c f = 0. Figure 15 shows a contour plot of the reattachment length l versus a 1c and a 2c , as well as various skin-friction profiles. From this figure it becomes clear that two different states exist. Fairly similar results are obtained throughout most of the [a 1c , a 2c ] domain. However, close to the axi-symmetric borders, distinctly different results are obtained. Here, the reattachment length is highly over predicted, leading to significantly larger recirculation bubbles. The maximum occurs at the 3C corner, with a reattachment length close to x/ h = 10. For comparison, we show the streamlines in Fig. 16 of the velocity field for representative samples of the aforementioned two states. The transition between the two states is very sharp. Let us examine this transition by visualizing the difference between trajectories for a 2c = and a 2c = 0, where is a small parameter and a 1c > 0 is kept constant. Figure 17 shows the trajectories in the barycentric map for = 10 -3 and = 10 -4 , computed from the flow in wall-normal direction at x/ h = 4. For a 2c = = 10 -3 , the source term for DC 2c /Dt is still active, and the corresponding trajectory is different from that of the full axi-symmetric state, which is also depicted in Fig. 17. This effect is still present for = 10 -4 , although much less pronounced. Hence, a sharp gradient exists in the space of the a αc coefficients as an axi-symmetric boundary is approached.

Data-free prediction

Let us briefly discuss the data-free results for the backward facing step, which were computed in the same way as discussed in Section 6.1.2. The results for several QoI are displayed in Fig. 18. For U 2 /U ref , there is a region where there is no overlap between our interval Fig. 15 Left: Skin-friction profiles and the computed reattachment length l. Note the different shape of the axi-symmetric and isotropic profiles compared to those computed with a 1c > 0 and a 2c > 0. Right: a contour plot of l and the experimental error bars. This indicates that, given the current SST baseline model, other uncertainties besides the shape of the Reynolds stresses must be taken into account if U 2 /U ref is considered important. For the other QoI the overlap is reasonable, although it occurs at the edges of the computed intervals. We can therefore already expect the posterior distributions to be skewed toward the edge of the prior domain.

Data-driven predictions

For this case, z will be experimental u 1 u 2 data from [START_REF] Jovic | Backward-facing step measurements at low reynolds number, re h = 5000[END_REF]. We performed 5 calibrations, at x/ h ∈ {4, 6, 10, 15, 19}, and the posterior distributions are plotted in Fig. 19. Despite the fact that we calibrated using data from both in-and outside the recirculation bubble, the posterior distributions are very similar for all x/ h. This is a desirable result, since we do not want the most-likely value of the a αc to change from one station to the next. Figure 19 shows that regarding u 1 u 2 , the baseline state is favored by the data since the joint posteriors are located at the top-right corner of the stochastic domain, see also Table 1. This is an indication that the use of a simple eddy-viscosity model is a reasonable choice for this particular flow, and that similar results might be obtained by only calibrating the closure coefficients of the a b QoIs is skewed toward the data. However, discrepancies between the validation data and the posterior distribution of U 2 /U ref can be observed, which comes as no surprise in light of the interval results of Section 6.2.3. Since we take care to include the baseline state within our prior, a better baseline model would be able to capture the data with less bias. Thus there can be an influence of the chosen baseline model, especially if the baseline state is preferred by the data. Given that the SST model is not perfect, other uncertainty measures might help to envelop the data better. We reiterate that we considered perturbations in the eigenvalues of b ij only. A more general UQ framework will also perturb the orientation of the Reynolds stress tensor. A study which introduced uniform perturbations in the eigenvectors can be found in [START_REF] Iaccarino | Eigenspace perturbations for uncertainty estimation of single-point turbulence closures[END_REF][START_REF] Mishra | Rans predictions fo high-speed flows using enveloping models[END_REF]. Introducing non-uniformity through transport equations is the subject of future work.

Conclusion

We introduced two model transport equations for the coefficients of the barycentric map, for the purpose of epistemic model-form uncertainty quantification of Boussinesq type eddy-viscosity closures. The model effectively adds a tensorial discrepancy function to the Boussinesq hypothesis, based on perturbations made in the eigenvalues of the anisotropy tensor. The discrepancy tensor is assumed to be transported in a similar fashion as other turbulent transport equations such as k and ω, yielding spatially varying perturbations. Through the inclusion of so-called lag terms, these eigenvalue perturbations are made with respect to a baseline eddy-viscosity model. Depending on the local flow physics, our model allows for the baseline state to be recovered.

The global behavior of the perturbation model is determined by two lag coefficients. These control the amount of perturbation from the baseline state, as well as the direction of perturbation in the barycentric map. This results in a low-dimensional inference problem, which is amenable to physical reasoning. This low-dimensional structure also provides us with two distinct UQ strategies, i.e. data-free and data-driven prediction. In the former we employ bounds on the eigenvalue perturbations to compute cheap intervals on the quantity of interest. In the data-driven approach, we use Bayesian inference to infer the posterior distribution of the two lag coefficients from any available data. Although more expensive, again due to the low dimensional nature of the UQ problem we can easily construct accurate surrogate models which can replace of the full RANS code in the likelihood function. With this in place, informed distributions of the lag coefficients are readily available.

Fig. 1

 1 Fig.1The barycentric map with the one-, two-and three-component corners, and the variation of its coefficients. Arrows indicates directions along which the corresponding coefficient linearly decreases from one toward zero

Fig. 4

 4 Fig. 4 Examples of | λ α | := |λ αλ (eq) α | for a subsonic jet (left), and a backward facing step (right)

Fig. 5 Fig. 6

 56 Fig. 5 Unrealizable turbulence. a The barycentric map for the baseline kε model when simulating the backwards-facing step flow of [25] (Re h = 5100, where h is the step height). b The spatial locations of unrealizable turbulence (C 3c < 0) for the kε model

Fig. 8

 8 Fig.8The experimentally-determined normalized axial velocity (top) and Reynolds stress component (bottom). Source:[START_REF] Bridges | The nasa subsonic jet particle image velocimetry (piv) dataset[END_REF]

Fig.Fig. 10 A

 10 Fig. Several SST baseline profiles (solid lines) along directions normal to the centerline, compared to the corresponding PIV data

Fig. 11

 11 Fig. 11 Experimental axial PIV velocity (squares), and the uncertainty intervals (shaded)

Fig. 14

 14 Fig. 14 The posterior U 1 /U jet trace of the jet case

Fig. 16 Fig. 17 Fig. 18 4 Fig. 19 Fig. 20

 16171841920 Fig. 16 Two distinctly different recirculation bubbles

  Compliance with Ethical StandardsThis investigation was funded by the United States Department of Energy's (DoE) National Nuclear Security Administration (NNSA) under the Predictive Science Academic Alliance Program II (PSAAP II) at Stanford University.

Table 1

 1

	Combinations of
	parameters used for computing	[a 1c , a 2c ]	state
	the intervals in the jet case
		[0, 0]	isotropic
		[1, 1]	baseline
		[0, 1]	max axisymmetric contraction
		[1, 0]	max axisymmetric expansion,

Table 2

 2

	The relative
	computational overhead t	[a 1c , a 2c ]	t
	compared to the runtime of the
	baseline SST model	[0, 0]	+ 32.3%
		[0, 1]	+ 30.8%
		[1, 0]	+ 28.1%
		[1, 1]	+ 32.6%

Conflict of interests

The authors declare that they have no conflict of interest.


Martial Briand
email: martial.briand@inra.fr

Mariam Bouzid
email: mariam.bouzid@icloud.com

Gilles Hunault
email: gilles.hunault@univ-angers.fr

Marc Legeay
email: legeay.marc@free.fr

Marion Le Saux
email: marion.le-saux@inra.fr

Matthieu Barret
email: matthieu.barret@inra.fr

Fischer-Le Saux

A rapid and simple method for assessing and representing genome sequences relatedness

Keywords: ANI, k-mers, circle packing, Pseudomonas, metagenome

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Background

Species is the unit of biological diversity. Species delineation of Bacteria and Archaea historically relies on a polyphasic approach based on a range of genotypic, phenotypic and chemo-taxonomic (e.g. fatty acid profiles) data of cultured specimen. According to the List of Prokaryotic Names with Standing in Nomenclature (LPSN), approximately 15,500 bacterial species names have been currently validated within this theoretical framework [START_REF] Parte | LPSN -List of Prokaryotic names with Standing in Nomenclature (bacterio.net), 20 years on[END_REF]. Since the number of bacterial species inhabiting planet Earth is predicted to range between 10 7 to 10 12 species according to different estimates [START_REF] Amann | After All, Only Millions?[END_REF][START_REF] Locey | Scaling laws predict global microbial diversity[END_REF], the genomics revolution provides an opportunity to accelerate the pace of species description. Prokaryotic species are primarily described as cohesive genomic groups and approaches based on similarity of whole genome sequence, also known as overall genome relatedness indices (OGRI), have been proposed for delineating species. Average nucleotidic identity (ANI) is nowadays the mostly acknowledged OGRI for assessing relatedness between genomic sequences. Distinct ANI algorithms such as ANI based on BLAST (ANIb [START_REF] Goris | DNA-DNA hybridization values and their relationship to whole-genome sequence similarities[END_REF]), ANI based on MUMmer (ANIm [START_REF] Richter | Shifting the genomic gold standard for the prokaryotic species definition[END_REF]) or ANI based on orthologous gene (OrthoANIb [START_REF] Lee | OrthoANI: An improved algorithm and software for calculating average nucleotide identity[END_REF]; OrthoANIu [START_REF] Yoon | A large-scale evaluation of algorithms to calculate average nucleotide identity[END_REF]; gANI,AF [START_REF] Varghese | Microbial species delineation using whole genome sequences[END_REF]), which differ in their precisions but more importantly on their calculation times [START_REF] Yoon | A large-scale evaluation of algorithms to calculate average nucleotide identity[END_REF], have been developed. Indeed, improvement of calculation time for whole genomic comparison of large datasets is an essential parameter. As of November 2018, the total number of prokaryotic genome sequences publicly available in the NBCI database is 170,728. Considering an average time of 1 second for calculating ANI values of one pair of genome sequence, it would take approximately 1,000 years for obtaining ANI values for all pairwise comparisons.

The number of words of length k (k-mers) shared between read sets [START_REF] Benoit | Multiple Comparative Metagenomics using Multiset k-mer Counting[END_REF] or genomic sequences [START_REF] Déraspe | Phenetic Comparison of Prokaryotic Genomes Using k-mers[END_REF] is an alignment-free alternative for assessing the dis(similarities) between entities. Methods based on k-mers counts, such as SIMKA [START_REF] Benoit | Multiple Comparative Metagenomics using Multiset k-mer Counting[END_REF], can quickly compute pairwise comparison of multiple metagenome read sets with high accuracy. In addition, specific kmers profiles are now routinely employed by multiple read classifiers for estimating the taxonomic structure of metagenome read sets [START_REF] Wood | Kraken: ultrafast metagenomic sequence classification using exact alignments[END_REF][START_REF] Kim | Centrifuge: rapid and sensitive classification of metagenomic sequences[END_REF][START_REF] Ounit | CLARK: fast and accurate classification of metagenomic and genomic sequences using discriminative k-mers[END_REF]. While these k-mers based classifiers differ in term of sensitivity and specificity [START_REF] Sczyrba | Critical Assessment of Metagenome Interpretation -a benchmark of computational metagenomics software[END_REF], they rely on accurate genome databases for affiliating read to a taxonomic rank.

The objective of the current work was to evaluate an alternative method based on kmers counts to study species delimitation on extensive genome datasets. We therefore decided to employ k-mers counts for assessing similarity between genome sequences belonging to the Pseudomonas genus. Indeed, this genus contains an important diversity of species (n = 207), whose taxonomic affiliation is under constant evolution [START_REF] Peix | Historical evolution and current status of the taxonomy of genus Pseudomonas[END_REF][START_REF] Berge | A user's guide to a data base of the diversity of Pseudomonas syringae and its application to classifying strains in this phylogenetic complex[END_REF][START_REF] Gomila | Clarification of Taxonomic Status within the Pseudomonas syringae Species Group Based on a Phylogenomic Analysis[END_REF][START_REF] Gomila | Phylogenomics and systematics in Pseudomonas[END_REF][START_REF] Peix | The current status on the taxonomy of Pseudomonas revisited: An update[END_REF][START_REF] Garrido-Sanz | Genomic and Genetic Diversity within the Pseudomonas fluorescens Complex[END_REF][START_REF] Hesse | Genome-based evolutionary history of Pseudomonas spp[END_REF] was selected for this comparison. This subset was composed of genome sequences containing less than 150 scaffolds.

Development of KI-S tool

An integrative tool named KI-S was developed. The number of shared k-mers between genome sequences is first calculated with Simka [START_REF] Benoit | Multiple Comparative Metagenomics using Multiset k-mer Counting[END_REF]. A custom R script is then employed to cluster the genome sequences according to their connected components at different selected threshold (e.g. 50% of shared 15-mers).

Taxonomic inference of metagenomic read sets

The taxonomic profiles of 9 metagenome read sets derived from seed, germinating seeds and seedlings of common bean (Phaseolus vulgaris var. Flavert) were estimated with Clark version 1.2.4 [START_REF] Ounit | CLARK: fast and accurate classification of metagenomic and genomic sequences using discriminative k-mers[END_REF]. These metagenome datasets were selected because of the high relative abundance of reads affiliated to Pseudomonas [START_REF] Torres-Cortés | Functional Microbial Features Driving Community Assembly During Seed Germination and Emergence[END_REF]. The following Clark default parametersk 31 -t <minFreqTarget> 0 and -o <minFreqtObject> 0 were used for the taxonomic profiling.

Results

Selection of optimal k-mers size and percentage of shared k-mers

Using the percentage of shared k-mers as an OGRI for species delineation first required to determine the optimal k-mers size. This was performed by comparing the percentage of shared k-mers to a widely acknowledged OGRI, ANIb [START_REF] Goris | DNA-DNA hybridization values and their relationship to whole-genome sequence similarities[END_REF], between 934 Pseudomonas genome sequences. Since species delineation threshold was initially proposed following the observation of a gap in the distribution of pairwise comparison values [START_REF] Patrick | Use of DNA reassociation in bacterial classification[END_REF], the distribution profiles obtained with k-mers lengths ranging from 10 to 20 were compared to ANIb values.

Short k-mers (k < 12) were evenly shared by most strains and then not discriminative (Fig. 1). As the size of the k-mers increased, a multimodal distribution based on four peaks were observed (Fig. 1). The first peak is related to genomes sequences that do not belong to the same species. Then, depending on k length, the second and third peaks (e.g. 50% and 80%

for k = 15) corresponded to genome sequences associated to the same species and subspecies, respectively. The fourth peak at 100% of shared k-mers was related to identical genome sequences.

Fifty percent of 15-mers is closed to ANIb value of 0.95 (Fig. 2), a threshold commonly employed for delineating bacterial species level [START_REF] Goris | DNA-DNA hybridization values and their relationship to whole-genome sequence similarities[END_REF]. More precisely the median percentage of shared 15-mers is 49% [34%-66%] for ANIb value ranging from 0.94 to 0.96. In addition, 15-mers allows the investigation of inter-and infra-specific relationship at lower and higher percentage of shared 15-mers, respectively.

Computing time of 15-mers for 934 genome sequences was 4 hours on a DELL Power Edge R510 server, while it took approximately 3 months for obtaining all ANIb pairwise comparisons (500-fold decrease of computing time).

Classification of Pseudomonas genomes

The percentage of shared 15-mers was then used to investigate relatedness between 3,623

Pseudomonas genome sequences publicly available. At a threshold of 50% of 15-mers, we identified 350 cliques. The clique containing the most abundant number of genome sequences was by far related to P. aeruginosa species (n = 2,341), followed by the phylogroups PG1 (n = 111), PG3 (n = 92) and PG2 (n = 74) of P. syringae species complex ( [START_REF] Berge | A user's guide to a data base of the diversity of Pseudomonas syringae and its application to classifying strains in this phylogenetic complex[END_REF]; Table S1). At the clustering threshold employed, 185 cliques were composed of a single genome sequence, therefore highlighting the high Pseudomonas strain diversity.

Moreover, according to Chao1 index, Pseudomonas species richness is estimated at 629 cliques [+ 57], which indicates that additional strain isolations and sequencing effort are needed to cover the whole diversity of this bacterial genus. Graphical representation of hierarchical clustering by dendrogram for a large dataset is generally not optimal. Here we employed Zoomable circle packing as an alternative to dendrogram for representing similarity between genome sequences (Fig. 3 and FigS1.html). The different clustering thresholds that can be superimposed on the same graphical representation allow the investigation of interand intra-groups relationships (Fig. 3 and FigS1.html). This is useful for affiliating specific clique to a group or subgroup of Pseudomonas species.

Improvement of taxonomic affiliation of metagenomic read sets.

The taxonomic composition of metagenome read sets is frequently estimated with k-mers based classifiers. While these k-mers based classifiers differ in term of sensitivity and specificity, they all rely on accurate genome databases for affiliating read to a taxonomic rank. Here, we investigated the impact of database content and curation on taxonomic affiliation. Using Clark [START_REF] Ounit | CLARK: fast and accurate classification of metagenomic and genomic sequences using discriminative k-mers[END_REF] as a taxonomic profiler with the original Clark database, we classified metagenome read sets derived from bean seeds, germinating seeds and seedlings [START_REF] Torres-Cortés | Functional Microbial Features Driving Community Assembly During Seed Germination and Emergence[END_REF]. Adding the 3,623 Pseudomonas genomes with their original taxonomic affiliation from NCBI to the original Clark database did not increase the percentage of classified reads (Fig. 4). However, adding the same genome sequences reclassified in cliques according to their percentage of shared k-mers (k=15; threshold= 50%) increased 1.4-fold on average the number of classified reads (Fig. 4).

Discussion

Classification of bacterial strains on the basis on their genome sequences similarities has emerged since a decade as an alternative to the cumbersome DNA-DNA hybridizations [START_REF] Rosselló-Móra | Past and future species definitions for Bacteria and Archaea[END_REF].

Although ANIb is the current gold-standard method for investigating these genomic relatedness, its intensive computational time prohibited its used for comparing large genome datasets [START_REF] Yoon | A large-scale evaluation of algorithms to calculate average nucleotide identity[END_REF]. In contrast, investigating the percentage of shared k-mers is scalable for comparing thousands of genome sequences.

In a method based on k-mers counts, choosing the length of k is a compromise between accuracy and speed. The distribution of shared k-mers values between genome sequences is impacted by k length. For k = 15, four peaks were observed at 15%, 50%, 80% and 100% of shared k-mers. The second peak is closed to ANIb value of 0.95 and falls in the so called grey or fuzzy zone [START_REF] Rosselló-Móra | Past and future species definitions for Bacteria and Archaea[END_REF] where taxonomists might decide to split or merge species.

Hence, according to our working dataset, it seems that 50% of 15-mers is a good proxy for estimating Pseudomonas clique. Despite the diverse range of habitats colonized by different

Pseudomonas populations [START_REF] Peix | The current status on the taxonomy of Pseudomonas revisited: An update[END_REF], it is likely that the percentage of shared k-mers has to be adapted when investigating other bacterial genera. Indeed, since population dynamics, lifestyle and location impact molecular evolution, it is somewhat illusory to define a fixed threshold for species delineation [START_REF] Bromham | Why do species vary in their rate of molecular evolution?[END_REF]. While 15-mers is a good starting point for investigating infra-specific to infra-generic relationships between genome sequences, the computational speed of KI-S offers the possibility to perform large scale genomic comparisons at different k sizes to select the most appropriate threshold.

Genomic relatedness using whole genome sequences becomes a standard for bacterial strain identification and bacterial taxonomy [START_REF] Rosselló-Móra | Past and future species definitions for Bacteria and Archaea[END_REF][START_REF] Parks | A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life[END_REF]. This proposition is primarily motivated by fast and inexpensive sequencing of bacterial genome together with the limited availability of cultured specimen for performing classical polyphasic approach. Whether full genome sequences should represent the basis of taxonomic classification is an ongoing debate between systematicians [START_REF] Garrity | A New Genomics-Driven Taxonomy of Bacteria and Archaea: Are We There Yet?[END_REF]. While this consideration is well beyond the objectives of this work, obtaining a classification of bacterial genome sequences into coherent groups is of general interest. Indeed, the number of misidentified genomes sequences is exponentially growing in public databases. A number of initiatives such as Digital Protologue Database

(DPD [START_REF] Rossello-Mora | Reflections on the introduction of the Digital Protologue Database -A partial success?[END_REF]), Microbial Genomes Atlas (MiGA [START_REF] Rodriguez-R Lm | The Microbial Genomes Atlas (MiGA) webserver: taxonomic and gene diversity analysis of Archaea and Bacteria at the whole genome level[END_REF]), Life Identification Numbers database (LINbase [START_REF] Vinatzer | A proposal for a portal to make earth's microbial diversity easily accessible and searchable[END_REF]) or the Genome Taxonomy Database (GTDB [START_REF] Parks | A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life[END_REF]) proposed services to classify and rename bacterial strains based ANIb values or single copy marker proteins.

Using the percentage of shared k-mers between unknown bacterial genome sequences and reference genome sequences associated to these databases could provide a rapid complementary approach for bacterial classification. Moreover, KI-S tool, provides a friendly visualization interface that could help systematicians to curate whole genome databases.

Indeed, zoomable circle packing could be employed for highlighting (i) misidentified strains, (ii) bacterial taxa that possess representative type strains or (iii) bacterial taxa that contain few genomes sequences.

Association between a taxonomic group and its distribution across a range of habitats is useful for inferring the role of this taxa on its host or environment. For instance, community profiling approaches based on molecular marker such as hypervariable regions of 16S rRNA gene have been helpful for highlighting correlations between host fitness and microbiome composition. Finer-grained taxonomic resolution of microbiome composition could be achieved with metagenomics through k-mers based classification of reads. In this study we demonstrate that employing a database with a classification of strains reflecting their genomic relatedness greatly improve taxonomic assignments of reads. Therefore, investigating relationship between bacterial genome sequences not only benefits bacterial taxonomy but also deserves microbial ecology. mers. Each dot represents a genome sequence, which is colored according to its group of species [START_REF] Berge | A user's guide to a data base of the diversity of Pseudomonas syringae and its application to classifying strains in this phylogenetic complex[END_REF][START_REF] Hesse | Genome-based evolutionary history of Pseudomonas spp[END_REF]. These genome sequences have been grouped at three distinct thresholds for assessing infraspecific (0.75), species-specific (0.5) and interspecies relationships (0.25). TableS1.csv : Pseudomonas cliques. Description of the 350 cliques obtained after clustering at 50% of shared 15-mers. For each clique, the Pseudomonas group [START_REF] Hesse | Genome-based evolutionary history of Pseudomonas spp[END_REF] and subgroup [START_REF] Berge | A user's guide to a data base of the diversity of Pseudomonas syringae and its application to classifying strains in this phylogenetic complex[END_REF][START_REF] Hesse | Genome-based evolutionary history of Pseudomonas spp[END_REF] are displayed.

Figures and Supplemental files

FigureS1.html: Zoomable circle packing representation of Pseudomonas genome sequences. Similarities between genome sequences were assessed by comparing the percentage of shared 15-mers. Each dot represents a genome sequence, which is colored according to its group of species [START_REF] Berge | A user's guide to a data base of the diversity of Pseudomonas syringae and its application to classifying strains in this phylogenetic complex[END_REF][START_REF] Hesse | Genome-based evolutionary history of Pseudomonas spp[END_REF]. These genome sequences have been grouped at three distinct thresholds for assessing infraspecific (0.75), species-specific (0.5) and interspecies relationships (0.25).

  Three distinct Clark databases were employed: (i) the original Clark database from NCBI/RefSeq at the species level (ii) the original Clark database supplemented with the 3,623 Pseudomonas genome sequences and their original NCBI taxonomic affiliation (iii) the original Clark database supplemented with the 3,623 Pseudomonas genome sequences whose taxonomic affiliation was corrected according to the reclassification based on the number of shared k-mers. For this third database, genome sequences were clustered at >50% of 15-mers.

Figure 1 :

 1 Figure 1: Distribution of shared k-mers values. Relatedness between genome sequences were estimated with ANIb (green) or shared k-mers (blue). The x axis represents ANIb or percentage of shared k-mers while the y axis represents the number of values by class in the subset of 934 Pseudomonas genomic comparison.

Figure 2 :

 2 Figure 2: Comparison of various k-mers length and ANIb values. Pairwise similarities between genome sequences were assessed with average nucleotidic identity based on BLAST (ANIb, x-axis) and percentage of shared k-mers of length 10 (A), 15 (B) and 20 (C).The red line corresponds to ANIb of 0.95, a threshold commonly employed for delineating species level.

Figure 3 :

 3 Figure 3: Hierarchical clustering of Pseudomonas genome sequences. Zoomable circle packing representation of Pseudomonas genome sequences (n = 3,623). Similarities between genome sequences were assessed by comparing the percentage of shared 15-

Figure 4 :

 4 Figure 4: Percentage of classified reads. Classification of metagenome read sets derived from bean seeds, germinating seeds and seedlings with Clark [13]. Three distinct databases were employed for read classification: the original Clark database (red), Clark database supplemented with 3,623 Pseudomonas genome sequences (green) and the Clark database supplemented with 3,623 Pseudomonas genome sequences that were classified according to their percentage of shared k-mers (blue).

  The clustering result is visualized with Zoomable

	Circle	Packing	representation	with	the	D3.js	JavaScript	library

(https://gist.github.com/mbostock/7607535). The source code of the KI-S tool is available at the following address: https://sourcesup.renater.fr/projects/ki-s/. A wrapper for accessing KI-S in a user-friendly Galaxy tool is also available at the following address: https://iris.angers.inra.fr/galaxypub-cfbp.

Acknowledgements

The authors wish to thank Claire Lemaitre and Guillaume Rizk for their assistance with the SIMKA software.

Funding

This research was supported in part by grant awarded by the Region des Pays de la Loire (metaSEED, 2013 10080).

Competing interests

The authors declare that they have neither competing interests nor conflict of interest.

KI-S Circlepack