Valentin Resseguier
email: valentin.resseguier@scalian.com

Agustin M Picard

Etienne Mémin

Bertrand Chapron

Quantifying truncation-related uncertainties in unsteady fluid dynamics reduced order models *

Keywords: Fluid dynamics, reduced order model, uncertainty quantification, stochastic closure, proper orthogonal decomposition AMS subject classifications. 60H35, 65M60, 65M75, 76M35, 93B11

In this paper, we present a new method to quantify the uncertainty introduced by the drastic dimensionality reduction commonly practiced in the field of computational fluid dynamics, the ultimate goal being to simulate accurate priors for real-time data assimilation. Our key ingredient is a stochastic Navier-Stokes closure mechanism that arises by assuming random unresolved flow components. This decomposition is carried out through Galerkin projection with a Proper Orthogonal Decomposition (POD-Galerkin) basis. The residual velocity fields, model structure and evolution of coefficients of the reduced order's solutions are used to compute the resulting multiplicative and additive noise's correlations. The low computational cost of these consistent correlation estimators makes them applicable to the study of complex fluid flows. This stochastic POD-ROM is applied to 2D and 3D DNS wake flows at Reynolds 100 and 300, respectively, with Uncertainty Quantification (UQ) and forecasting outside the learning interval being the main focus. The proposed stochastic POD-ROM approach is shown to stabilize the unstable temporal coefficients and to maintain their variability under control, while exhibiting an impressively accurate predictive capability.

1. Introduction. The industrial application of partial differential equations (PDE)driven processes -fluid dynamics, for instance -can be a daunting task, mainly due to the computational complexity involved with its resolution. This computational burden becomes even excruciatingly difficult to handle when it comes to achieving real-time simulations. To tackle this difficulty, reduced order models (ROM) are commonly employed to speed up deterministic and stochastic design simulations [START_REF] Boyaval | Reduced basis techniques for stochastic problems[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF], or optimal control problems [START_REF] Chen | Adaptive reduced-order-model-based control-law design for active flutter suppression[END_REF][START_REF] Gunzburger | Optimal control of stochastic flow over a backwardfacing step using reduced-order modeling[END_REF][START_REF] Leclercq | Linear iterative method for closed-loop control of quasiperiodic flows[END_REF]. The interested community has proposed a plethora of different algorithms to reduce the computational cost of stochastic PDEs for uncertainty quantification (UQ) applications [START_REF] Chen | A weighted reduced basis method for elliptic partial differential equations with random input data[END_REF][START_REF] Chen | Reduced basis methods for uncertainty quantification[END_REF][START_REF] Le Maitre | A stochastic projection method for fluid flow. II. random process[END_REF][START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF][START_REF] Nouy | Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations[END_REF][START_REF] Sapsis | Blending modified Gaussian closure and non-Gaussian reduced subspace methods for turbulent dynamical systems[END_REF][START_REF] Torlo | Stabilized weighted reduced basis methods for parametrized advection dominated problems with random inputs[END_REF]. In particular, when some parameters of deterministic or stochastic PDEs are random, the weighted reduced basis (wRB) method [START_REF] Chen | A weighted reduced basis method for elliptic partial differential equations with random input data[END_REF][START_REF] Chen | Reduced basis methods for uncertainty quantification[END_REF] helps choosing the best parameters' value for each new full-order simulation. With the addition of stabilizing terms, [START_REF] Torlo | Stabilized weighted reduced basis methods for parametrized advection dominated problems with random inputs[END_REF] applied this method to linear advection-dominated problems, including advection by a random velocity field.

In turbulent fluid dynamics, the system's energy usually spreads out over many degrees of freedom. This prevents low-dimensional approximations of the state vector from being sufficiently accurate, but such a rough approximations can theoretically be sufficient for specific industrial applications, especially when the quantity of interest (QoI) is a spatial average (e.g. lift and drag). However, that multiscale property prevents, in general, the constitution of accurate reduced systems for the dynamics of those low-dimensional state approximations. Indeed, severe modal truncations often usually end up destabilizing the system and overdamping some of the stable coefficients of the reduced order solution [START_REF] Sapsis | Blending modified Gaussian closure and non-Gaussian reduced subspace methods for turbulent dynamical systems[END_REF].

Consequently, to stabilize the ROM, authors commonly introduce an additional deterministic term (typically an eddy viscosity term) [START_REF] Aubry | The dynamics of coherent structures in the wall region of a turbulent boundary layer[END_REF][START_REF] Cazemier | Proper orthogonal decomposition and low-dimensional models for driven cavity flows[END_REF][START_REF] Wang | Proper orthogonal decomposition closure models for turbulent flows: a numerical comparison[END_REF], along with a possible calibration on available data [START_REF] Buffoni | Low-dimensional modelling of a confined three-dimensional wake flow[END_REF][START_REF] Cordier | Identification strategies for model-based control[END_REF][START_REF] Xie | Data-driven filtered reduced order modeling of fluid flows[END_REF]. This calibration procedure can be extended to the complete set of the ROM's parameters [START_REF] Artana | Strong and weak constraint variational assimilations for reduced order fluid flow modeling[END_REF][START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Perret | Polynomial identification of POD based low-order dynamical system[END_REF][START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. ROM performances have mainly been evaluated on low Reynolds-number flows (say Re 100) [e.g. [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Stabile | POD-Galerkin reduced order methods for CFD using finite volume discretisation: vortex shedding around a circular cylinder[END_REF][START_REF] Xie | Data-driven filtered reduced order modeling of fluid flows[END_REF], while for flows at much greater Reynolds numbers, evaluations close to the learning time interval and reduced dynamics for two-dimensional flow observables (e.g. particle image velocimetry) have usually been the norm [e.g. [START_REF] Artana | Strong and weak constraint variational assimilations for reduced order fluid flow modeling[END_REF][START_REF] Buffoni | Low-dimensional modelling of a confined three-dimensional wake flow[END_REF][START_REF] Cordier | Identification strategies for model-based control[END_REF][START_REF] Stabile | POD-Galerkin reduced order methods for CFD using finite volume discretisation: vortex shedding around a circular cylinder[END_REF]. Long time ROM predictions have been otherwise performed with a substantial number of modes [e.g. [START_REF] Carlberg | Efficient non-linear model reduction via a least-squares Petrov-Galerkin projection and compressive tensor approximations[END_REF][START_REF] Stabile | A reduced order variational multiscale approach for turbulent flows[END_REF]]. Yet, turbulent flow ROMs remain inexact and uncontrolled in the long run, owing to their intrinsic chaotic nature and the growth of accumulated error along time. Outside the learning time interval, predictions becomes less and less accurate. Moreover, these predictions necessitate accurate initial conditions, forcing and geometric information, which are often poorly known.

Ensemble-based data assimilation -such as particle filters or ensemble Kalman filters -can alleviate these issues by forecasting an ensemble of likely future states of the system, while sequentially constraining them with on-coming measurements [START_REF] Doucet | A tutorial on particle filtering and smoothing: Fifteen years later[END_REF][START_REF] Evensen | Data assimilation: The ensemble Kalman filter[END_REF]. Still, this necessitates accurate quantification of their associated simulation's uncertainties. In this context, the aim is not to reduce the dimensionality of UQ, but rather to quantify the uncertainty introduced by the dimensionality reduction.

Note, UQ is a recurrent issue in applied fluid dynamics, and many strategies have been proposed for incorporating randomness in the physical models through some of their parameters [e.g. [START_REF] Le Maitre | A stochastic projection method for fluid flow. II. random process[END_REF][START_REF] Lucarini | Mathematical and physical ideas for climate science[END_REF]. However, the error introduced by these noisy parameters is not a priori tied to dimensionality reduction and to the contribution of unresolved components.

In particular, in fluid dynamics, random initial conditions have first been widely used for both UQ and predictability studies [e.g . 51]. It was later demonstrated this yields underdispersive quantification, i.e. it has a tendency to underestimate the error associated to the dimensionality reduction [START_REF] Berner | Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations[END_REF][START_REF] Franzke | Stochastic climate theory and modeling[END_REF][START_REF] Gottwald | The role of additive and multiplicative noise in fil-tering complex dynamical systems[END_REF][START_REF] Mitchell | Data assimilation in slow-fast systems using homogenized climate models[END_REF].

Alternatively, authors have considered the introduction of additive noises, most likely beginning with the introduction of EDQNM [START_REF] Leith | Atmospheric predictability and two-dimensional turbulence[END_REF][START_REF] Orszag | Analytical theories of turbulence[END_REF]. Without special care, such a strategy often leads to energy conservation loss, stability issues and radical changes of the underlying attractor [START_REF] Chapron | Large-scale flows under location uncertainty: a consistent stochastic framework[END_REF].

The Modified Quasilinear Gaussian (MQG) method [START_REF] Sapsis | Blending modified Gaussian closure and non-Gaussian reduced subspace methods for turbulent dynamical systems[END_REF][START_REF] Sapsis | Statistically accurate low-order models for uncertainty quantification in turbulent dynamical systems[END_REF][START_REF] Sapsis | A statistically accurate modified quasilinear Gaussian closure for uncertainty quantification in turbulent dynamical systems[END_REF] approximates the thirdorder moment to help redistribute the right amount of energy between the coefficients of the reduced solution. Finally, several techniques related to averaging and homogenization theory exploit a time-scale separation hypothesis [START_REF] Gottwald | Stochastic climate theory[END_REF][START_REF] Kurtz | A limit theorem for perturbed operator semigroups with applications to random evolutions[END_REF][START_REF] Mitchell | Data assimilation in slow-fast systems using homogenized climate models[END_REF][START_REF] Papanicolaou | Asymptotic theory of mixing stochastic ordinary differential equations[END_REF], one of the most notable being the MTV model [START_REF] Majda | A mathematical framework for stochastic climate models[END_REF]. This latter approach can reproduce intermittency and extreme events alike, thanks to its correlated additive and multiplicative noises. However, besides potential energy conservation issues, the noise covariance is often not explicit enough, and has to be simplified and estimated using the available data. Interested readers can refer to [START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF], and references herein, for more detailed reviews on model error specification in coarse-scale computational fluid dynamics (CFD). For ROM UQ, [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-dimensional and high-dimensional nonlinear models[END_REF] propose distributions and efficient sampling methods for the projection matrices in Galerkin-projection-based dimensionality reduction methods. Although apparently more pertinent than methods based on randomized parameters, the relation of this technique to errors associated to mode truncation and turbulent chaotic behaviour remains unclear.

In this paper, we propose to quantify the uncertainty introduced by modal dimensionality reduction through the so-called dynamics under location uncertainty (LU) [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, part I random transport and general models[END_REF]. Specifically, we adapt the aforementioned stochastic closure to the Galerkin-projectionbased ROM. Inspired from the theoretical works of [START_REF] Brzeźniak | Stochastic partial differential equations and turbulence[END_REF][START_REF] Mikulevicius | Stochastic Navier-Stokes equations for turbulent flows[END_REF], the LU closure relies on the stochastic transport of the flow variables, together with a decorrelation assumption of the unresolved fluctuations with respect to the resolved slow/large scales. More precisely, the residual velocity -i.e. the difference between the usual Navier-Stokes solution v and some large-scale velocity component w -is assumed to be time-uncorrelated at the characteristic time of the large-scale processes. This residual velocity is informally denoted v = σ Ḃ = σdB t /dt where t → σB t is a Q-Wiener process [START_REF] Da Prato | Stochastic Equations in Infinite Dimensions, Encyclopedia of Mathematics and its Applications[END_REF][START_REF] Prévôt | A concise course on stochastic partial differential equations[END_REF], and hence, Gaussian in nature. Note that this apparent simplified Gaussian assumption leads to, as we will see, a non-Gaussian multiplicative noise in the dynamics. Spatial correlations of the residual velocity are then specified through the Hilbert-Schmidt integral operator σ with a C 2 kernel σ:

Ω 2 → R d×d : σ(x)dB t = Ω σ(x, z)dB t (z)dz ∀(x, t) ∈ Ω × [0, T ]. (1.1)
This operator -or equivalently the spatial covariance of the residual velocity -can be modeled or learned from data. The review [START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF] describes some of the many choices that have been explored in this vein. This includes for instance parametric models based on fluid velocity self-similarity or brute-force non-parametric covariance estimation from highresolution datasets. Here, B t is an I d -cylindrical Wiener process and dB t /dt plays the role of spatio-temporal white noise. The above definition enables us to characterize the way physical quantities are transported by the stochastic flow:

dX t = w(X t , t)dt + σ(X t )dB t ∀t ∈ [0, T ], (1.2)
This resembles the expression for transport in classical fluid dynamics. Material derivatives and other differential operations of fluid dynamics are then derived through the use of stochastic PDEs (SPDE), in particular, by applying the Itō-Wentzell formula [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF] and a stochastic version of the Reynolds transport theorem [START_REF] Harouna | Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling[END_REF][START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, part I random transport and general models[END_REF]. As such, LU models can be applied to model error quantification [START_REF] Chapron | Large-scale flows under location uncertainty: a consistent stochastic framework[END_REF][START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, part II quasi-geostrophy and efficient ensemble spreading[END_REF][START_REF] Resseguier | Data-driven versus self-similar parameterizations for stochastic advection by lie transport and location uncertainty[END_REF], to improve large-scale simulations [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF][START_REF] Bauer | Stochastic representation of mesoscale eddy effects in coarse-resolution barotropic models[END_REF][START_REF] Harouna | Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, part II quasi-geostrophy and efficient ensemble spreading[END_REF][START_REF] Resseguier | Data-driven versus self-similar parameterizations for stochastic advection by lie transport and location uncertainty[END_REF], to reduced order modeling and data analysis [START_REF] Resseguier | Stochastic modelling and diffusion modes for proper orthogonal decomposition models and small-scale flow analysis[END_REF] or to data assimilation purposes [START_REF] Chandramouli | Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty[END_REF][START_REF] Yang | High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar[END_REF] in geophysical fluid dynamics and CFD. To note, in the geometric mechanics community [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF][START_REF] Cotter | Numerically modeling stochastic lie transport in fluid dynamics[END_REF], the Stochastic Advection by Lie Transport (SALT) method has been derived for large-scale modeling and data assimilation [START_REF] Cotter | Sequential monte carlo for stochastic advection by Lie transport (SALT): A case study for the damped and forced incompressible 2d stochastic euler equation[END_REF]. Both frameworks have been compared, numerically [START_REF] Bauer | Stochastic representation of mesoscale eddy effects in coarse-resolution barotropic models[END_REF][START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF] and conceptually, [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF][START_REF] Resseguier | Data-driven versus self-similar parameterizations for stochastic advection by lie transport and location uncertainty[END_REF], with LU and SALT exhibiting different conservation properties, namely energy preservation and circulation conservation, respectively. Applied to a barotropic Quasi-Geostrophic model, LU leads to improved accuracy when compared to a classical large-scale deterministic framework or to the circulation conservation stochastic setup [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF][START_REF] Bauer | Stochastic representation of mesoscale eddy effects in coarse-resolution barotropic models[END_REF]. The LU setting also fully captures the structural deformation of the large-scale flow component by the spatial inhomogeneity of the small-scales [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF]. It is important to note that these properties are independent of the choice of stochastic integral.

In this paper, our focus is to analyse the LU setting to help define efficient, highly reduced order models for real time data assimilation or control applications. For the sake of simplicity, we will deal with Proper Orthogonal Decomposition (POD) [START_REF] Lumley | Coherent structures in turbulence[END_REF] as dimensionality reduction technique -where time and space dependency are separatedbut the proposed methodology applies to all types of modal decompositions.

In section 2, we introduce the stochastic fluid dynamics closure we will employ throughout the paper, as well as its algebraic structure, followed by the Galerkin projection of this SPDE to derive our stochastic ROM, with a brief recall of the principle behind POD-ROMs in section 3. Several estimations needed to complete our stochastic POD-ROM are detailed in section 4, as well as the efficient and consistent estimators we exploit to rely both on data and on the closure's physical grounding. In section 5, we discuss the conservative properties of LU dynamics and of its reduced versions, and finally, section 6 is dedicated to the numerical evaluation of the UQ capabilities of our ROM.

T ∈ R * + and L 2 (Ω) (resp. L 2 = L 2 (Ω) d
) the space of square-integrable scalar (resp.

vector) fields on Ω. In the case of incompressible fluids, the LU Navier-Stokes equations on Ω × [0, T ] read:

D t w Stochastic transport = -∇ (p dt + p σ dB t )
Pressure forcing

+ 1 Re ∆(wdt + σdB t ) Molecular viscous dissipation , (2.1) 0 = ∇ • (w * dt + σdB t ) Mass conservation , (2.2)
where p + p σ dBt dt informally represents the pressure also decomposed into a large-scale and a small-scale time-uncorrelated component, Re stands for the Reynolds number, and for every smooth-enough function q : Ω × [0, T ] → R d , we denote:

(D t q) k = d t q k = q k (x,t+dt)-q k (x,t) Time increment + (w * dt + σdB t ) • ∇q k Advection -∇ • 1 2 a∇q k Turbulent diffusion
dt, (2.3)

w * = w -1 2 (∇ • a) T , (2.4) a = E (σdB t ) (σdB t ) T /dt = σ σ, (2.5)
where operator is defined later on in this section.

Throughout this paper, we consider various linear integral operators defined on some subspace D ⊂ L 2 . Denoted η, each of those operators can be defined through a matrix kernel η of L 2 Ω 2 d×d as follows:

(ηξ) (x) = Ω η(x, z) kernel of η ξ(z)dz ∀ξ ∈ D ⊂ L 2 , ∀x ∈ Ω. (2.6)
Two first examples of such operator are σ encoding the small-scale velocity component spatial correlations (defined in (1.1)) and p σ encoding the small-scale pressure component spatial correlations. For two such linear integral operators η and θ, we note:

(η θ)(x) = Ω η(x, z) θ T (x, z)dz ∀x ∈ Ω. (2.7)
Theoretical foundations to analyze SPDEs are outside the scope of this paper. Interested readers may refer to [START_REF] Flandoli | The interaction between noise and transport mechanisms in PDEs[END_REF][START_REF] Mikulevicius | Stochastic Navier-Stokes equations for turbulent flows[END_REF] for deeper insights. Here, we focus on the stochastic

Physical meaning Notation

Full-order term Property ROM term Molecular viscous is generic, and depending on the application, forces and boundary conditions, it may be modified accordingly to adjust to any incompressible flow configurations. For compressible stochastic flows, some new terms appear. Interested readers may refer to [START_REF] Resseguier | Geophysical flows under location uncertainty, part I random transport and general models[END_REF].

L 1 Re ∆ symmetric, < 0 l dissipation Turbulent F dif ∇ • 1 2 a∇• symmetric, < 0 fdif diffusion Advecting velocity F adv 1 2 (∇ • a) ∇ skew-symmetric fadv correction (if ∇ • (∇ • a) T = 0) Usual C(w, •) -(w • ∇) skew-symmetric c advection (if ∇ • w = 0)
In addition to the classical physical assumptions pertaining to the establishment of the Navier-Stokes equations for incompressible fluids, the main assumption of the LU setting is to consider the unresolved fluctuation velocity component uncorrelated in time. To note, the incompressible character of the random fluctuations can be relaxed at the price of additional terms in the transport operator [START_REF] Resseguier | Geophysical flows under location uncertainty, part I random transport and general models[END_REF] 2.2. Algebraic structure of the model. The algebraic properties of the different terms can be quickly described. We can formally rewrite the velocity evolution law (2.1) for t ∈ [0, T ] as follows:

d t w = (dM)(w) = ((L + F )(w) + C(w, w)) dt + (GdB t )(w) + (HdB t ) + dP , (2.8)
where table 1 details each term. Operator L represents the molecular viscosity term, F = F dif + F adv , the turbulent diffusion plus the advecting velocity correction, and GdB t , the advection by the random residual velocity. All of them are linear differential operators, while C, the term representing the usual non-linear advection effect, is a bilinear differential operator. The additive noise HdB t corresponds to the molecular viscous dissipation of the time-uncorrelated velocity component, σ Ḃ, while the last term on the right-hand side, dP , is the pressure forcing.

Additionally, under suitable boundary conditions, the algebraic structures of the different operators can be further specified. For instance, L and F dif are symmetric negative operators. For the other terms, additional incompressibility conditions are needed. The mass conservation equation (2.2) implies ∇ • σ = 0 1 , which is actually the case in practice, either because of the parametric model used for σ or because of the incompressibility of the velocity data used to estimate σ. Accordingly, GdB t is a skew-symmetric operator.

As a matter of fact, if (ζ, ξ) = Ω ζ • ξ denotes the scalar product of L 2 , for every ζ and ξ in the Sobolev space

H 1 0 (Ω) = {f ∈ L 2 : Ω ∇f T 2 < ∞; f | ∂Ω = 0} (with partial
derivative taken in the weak sense and • stands for the Euclidean norm on R d×d ), an integration by parts gives:

(ζ, (GdB t ) (ξ)) = - d k=1 Ω ζ k (σdB t • ∇) ξ k = - d k=1 Ω ζ k ∇ • (σdB t ξ k ) , (2.9) = d k=1 Ω (σdB t • ∇ζ k ) ξ k = -((GdB t ) (ζ), ξ). (2.10)
The mass conservation equation also implies ∇ • w * = 0 which may or may not be the case in practice. If ∇ • w = 0, operator C is skew-symmetric with respect to the second

argument -i.e. g → C(f , g) is skew-symmetric -whereas F adv is skew-symmetric if the drift correction w * -w = -1 2 (∇ • a) T is divergence-free.
Moreover, the turbulent diffusion F dif is related to the random skew-symmetric operator GdB t . Indeed, for every process ξ in

H 2 (Ω) = {f ∈ L 2 : ∂f ∂x i , ∂∂f ∂x i ∂x j ∈ L 2 , 1 i, j d}, (F dif (ξ)) k = ∇ • 1 2 a∇ξ k = ∇ • 1 2 (σ σ)∇ξ k = 1 2 (σ T ∇) T (σ T ∇) T ξ k , (2.11) = 1 2 G ( G(ξ k )) = 1 2 G ( G(ξ)) k = -1 2 G * ( G(ξ)) k , (2.12)
where G * denotes the adjoint of G. It also shows that F dif (w)dt = 1 2 d G(w), B t , where g, h denotes the quadratic covariation of any functions g and h. The diffusion term explicitly appears when working with the Itō stochastic integral and is only implicitly taken into account with Stratonovich integral [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF][START_REF] Resseguier | Data-driven versus self-similar parameterizations for stochastic advection by lie transport and location uncertainty[END_REF].

As discussed in more detail in section 5, these algebraic properties make the LU Navier-Stokes model -and to a certain extent, its reduced order versions -conservative (up to molecular viscosity and boundary conditions effects).

Galerkin projection.

To sample good priors for future Bayesian estimation algorithms, we aim at deriving a computationally efficient fluid dynamics ROM able to quantify its own errors with respect to the true fluid dynamics (i.e. the Navier-Stokes equations).

As previously mentioned, standard Galerkin techniques -even with the best deterministic closures -are hardly capable of such a goal as they were not originally designed for it. Hence, we propose to perform Galerkin projections on the LU Navier-Stokes model instead, and to study its appropriateness for this sort of tasks.

3.1.

A ROM with correlated additive and multiplicative noise. Let v be the real velocity field (i.e. the Navier-Stokes equations' solution) and φ 0 , a background velocity field, typically the velocity temporal mean v = 1 T T 0 v. To reduce the state space dimension, we project the fluid velocity anomaly, vφ 0 , in a subspace spanned by a number of

1 i.e. ∇x • σ(x, z) = 0, ∀(x, z) ∈ Ω 2 orthonormal spatial modes (φ i ) 1 i n . v(x, t, ω) = n i=0 b i (t, ω)φ i (x) =w R + Residual =v ∀(x, t) ∈ Ω × [0, T ], (3.1)
where b 0 = 1 and φ 0 = v by convention. The associated temporal coefficients, b i , are possibly random and depend on a realization ω of a sample space Ω whereas the reduced basis functions φ i , are assumed to be deterministic and stationary.

As typically done while working with ROMs, we aim at specifying the evolution of the projected velocity field w R . A standard technique for that is the Galerkin projection of the physical PDE -here, the Navier-Stokes equations -onto the reduced basis' functions φ i , where the resolved component w R is approximated by the solution of these projected equations. Mode truncation can create many problems. For moderately turbulent to turbulent flows and small dimension n, a closure model to handle the truncated modes is unavoidable. With the LU setting -described in section 2 -an elegant stochastic alternative for this problem can be derived. Accordingly, we will (i) assume that the residual velocity v is time-uncorrelated and is denoted σdB t /dt and (ii) approximate the resolved component w R to a realization of the solution of the Galerkin projection of the stochastic Navier-Stokes representation (2.1)-(2.2). The former hypothesis is a debatable choice with respect to some ROM's applications whose pertinence and limitations will be discussed further in section 4.5.

To obtain this stochastic ROM, the SPDE (2.8) is first projected onto the divergencefree function space through the non-local Leray operator P = I d -∇∇ T ∆ -1 . This projection, which requires the resolution of a Poisson equation is used to simplify the system by removing the pressure term. Since divergence-free solution is considered (see section 3.2), the resolved velocity component w R is naturally incompressible and we get

Pd t w R = d t w R , PL(w R ) = L(w R ), PdP = 0. (3.2)
Moreover, P (HdB t ) = (HdB t ) because of the necessary incompressibility of the Brownian term in the continuity equation (2.2)2 . Then, the resulting SPDE is projected onto each of the reduced basis' functions:

db i = (φ i , d t w R ) = (dM R i )(b) = φ i , P(dM)(w R ) , 1 i n, (3.3) where b = (b i ) 0 i n and (dM R i )(b) = n p=0 φ i , L(φ p ) =l pi b p dt + n p=0 φ i , PF (φ p ) = fpi b p dt + n p,q=0 φ i , PC(φ p , φ q ) =c pqi b p b q dt + (φ i , (HdB t )) =(θ i• dBt) + n p=0 φ i , P(GdB t )(φ p ) =(α pi• dBt) b p , 1 i n. (3.4)
The terms (α pi• dB t ) 1 p,i n and ((θ i• + α 0i• )dB t ) 1 i n correspond to a Gaussian skewsymmetric matrix and a Gaussian vector respectively, both with correlated coefficients.

If the functions φ i were spatial Fourier modes associated with small wave-numbers, the ROM (3.3) would represent a (stochastic) Large-Eddy-Simulation-like model expressed in Fourier space and b would be the set of Fourier coefficients of the solution. Hereafter, we will solely focus on the POD.

Proper Orthogonal Decomposition.

In the POD framework, the reduced basis' functions are computed through a set of velocity snapshots (v obs (•, t i )) 0 i N -1 . More precisely, they are a solution to the constrained optimization problem:

Maximize (φ i ) 1 i n n i=1 (φ i , v obs -v) 2 subject to φ i , φ j = δ ij , 1 i, j n. (3.5)
Reduced basis functions φ i are thus the n orthonormal functions which can best explain the snapshots' temporal variability. Similarly to a principal component analysis (PCA), the solutions of this optimization problem are the eigenfunctions of the velocity anomalies' (v obs -v) spatial covariance. Numerically, this matrix is extremely large, and we generally opt to solve the dual problem instead: we seek the eigenvalue decomposition of the velocity anomalies' temporal covariance:

C v ij = ((v obs -v) (•, t i ), (v obs -v) (•, t j )) , 0 i, j N -1. (3.6)
This method is often referred to as the snapshot method. In the POD framework, the mode φ 0 = v is set to the time averaged velocity, and the temporal coefficients energies b 2 i are denoted λ i . Furthermore, if the snapshots describe a divergence-free velocity field v obs , the spatial bases φ i are divergence-free as well.

4. Estimations of subgrid terms. To close our stochastic ROM (3.3)-(3.4), we need to estimate the variance tensor a (involved in the ROM matrix f = fdif + fadv defined in table 1 and in equation (3.4) under the braces) as well as the ROM noise variances and correlations. Firstly, if we recall that b 0 = 1, we note that:

(α •i• dB t ) T b + (θ i• dB t ) = n k=1 (α ki• dB t ) b k + ((θ i• + α 0i• )dB t ) , 1 i n. (4.1)
The multiplicative and additive noises of the ROM correspond to the first and second term of the right-hand side, respectively. To simplify notations, we write:

αpi• = α pi• + δ p0 θ i• , 1 i n, 0 p n. (4.2)
To fully specify the ROM, the following correlations must be estimated: involves the two-point quadratic cross-variation tensor of the small-scale velocity:

a(x) = E (σ(x)dB t ) (σ(x)dB t ) T /dt ∀x ∈ Ω, Σ α pi,qj = E {( αpi• dB t ) ( αqj• dB t )} /dt, 1 i, j n, 0 p, q n.
Q(x, y) = E (σ(x)dB t ) (σ(y)dB t ) T /dt ∀x, y ∈ Ω. (4.4)
The coefficients (i, j) of the covariance matrix of the Gaussian vector (θdB t ) can be expressed as follows:

E {(θ i• dB t ) (θ j• dB t )} = E φ i , ν∇ 2 σdB t ν∇ 2 σdB t , φ j , (4.5) = Ω 2 φ i (x) T ν 2 ∇ 2 x ∇ 2 y Q(x, y)dt φ j (y) dxdy. (4.6)
Since, by definition, ∀x ∈ Ω, a(x) = Q(x, x), the whole stochastic ROM is closed by specifying the tensor Q. In practice, this tensor is often heterogeneous -i.e. σB t is non-stationary in space. Thus, its spatially-discrete version is expected to be so large that its estimation or even its storage becomes prohibitive. To overcome this difficulty, a first strategy is to assume a model structure for the covariance as per [START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, part II quasi-geostrophy and efficient ensemble spreading[END_REF][START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF][START_REF] Resseguier | Data-driven versus self-similar parameterizations for stochastic advection by lie transport and location uncertainty[END_REF].

If the small-scale velocity is observed, other techniques can be considered. Indeed, one can build a POD representation of the small-scale velocity σ Ḃ, learned from available snapshots or realizations as in [START_REF] Bauer | Stochastic representation of mesoscale eddy effects in coarse-resolution barotropic models[END_REF][START_REF] Cotter | Numerically modeling stochastic lie transport in fluid dynamics[END_REF][START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF][START_REF] Resseguier | Data-driven versus self-similar parameterizations for stochastic advection by lie transport and location uncertainty[END_REF]. Depending on the desired accuracy for the noise representation, a great number of modes would have to be estimated and the number of coefficients involved in the ROM can quickly grow out of control. Here, we rely on a method specifically devised for this kind of ROM frameworks: the noise structure is again learned from observed residual velocity snapshots, but without passing through the covariance Q. It enables us to directly estimate the correlations of the random ROM's coefficients, as it will be explained later on.

Variance tensor estimation.

Considering the variance tensor a stationary leads to a simple estimator through time averaging as in [START_REF] Resseguier | Stochastic modelling and diffusion modes for proper orthogonal decomposition models and small-scale flow analysis[END_REF]:

a = ∆t v obs v obs T , (4.7)
where

v obs = v obs -w R obs = v obs - n i=0 b obs i φ i , (4.8)
is the observed POD's residual velocity. Since adt = d σB, (σB) T , the estimator (4.7) is consistent thanks to the quadratic covariation definition.

Noise correlation estimation.

Besides the variance tensor, there are n 2 (n + 1) 2 correlations to estimate according to (4.3). For any function ξ in H 2 (Ω) , let us introduce the linear functional:

K jq [ξ] = φ j , -P (ξ • ∇) φ q + δ q0 ν∆ξ , 1 j n, 0 q n. (4.9)
Using this notation, the noise's covariance can be estimated as follows:

Σ α pi,qj = ∆t
∆b obs i ∆t = ∆b obs i ∆t -b obs T (l + f ) •i + b obs T c ••i b obs , (4.14) ∆b obs i ∆t (t k ) = b obs i (t k + ∆t) -b obs i (t k ) ∆t , 0 k N -1. (4.15)
To ensure the noise's covariance matrix to have the desired symmetric non-negative structure, we only keep the symmetric part of the estimated tensor (4.10) and set its possible negative eigenvalues to zero. The part inside the functional K jq in the estimator (4.10) is inspired from [START_REF] Genon-Catalot | Non-parametric estimation of the diffusion coefficient by wavelets methods[END_REF], where products of martingale time increments are projected onto orthogonal functions of L 2 ([0, T ]). For ROMs with a small number of dimensions n, the estimator's computational cost is remarkably low as the observed coefficients of the reduced solution b obs were already computed by the method of snapshots (see section 3.2). Hence, the computational cost for b obs

Σ α pi,qj = lk σ α pi,lk σ α qj,lk , 1 i, j n, 0 p, q n, (4.16)
and O(n 2 ) independent white noises, one can sample realizations of the ROM Gaussian noise terms αdB t . However, σ α has O(n 4 ) coefficients, while the ROM works with only n modes and involves -leaving the noise terms aside -O(n 3 ) coefficients. Therefore, we propose to reduce the noise dimension through a tensorial PCA of Σ α , eventually only keeping the n first eigenvectors. This leads to the following sampling strategy:

αdB t ≈ n k=1 αR k dβ (k)
t , (4.17 It is important to note that the methodology described here is different from a more usual methodology based on PCA decomposition of the residual velocity v keeping n modes and assuming that the corresponding temporal coefficients (b i ) n+1 i 2n are timedecorrelated. The complexity of the final ROM is the same in both methodologies but our ROM maximizes the noise's variance instead of the residual velocity's variance. Thus, our proposed method is better suited in terms of ROM UQ. This assumption is consistent with the fact that the higher-order coefficients of the reduced order solution often tend to have shorter correlation time in fluid dynamics systems. However, in practice, this assumption is not found to be true, and it is a recurrent issue for the datadriven modeling of systems combining fast and slowly evolving components [START_REF] Azencott | Sub-sampling and parametric estimation for multiscale dynamics[END_REF][START_REF] Azencott | Adaptive sub-sampling for parametric estimation of Gaussian diffusions[END_REF][START_REF] Papavasiliou | Maximum likelihood drift estimation for multiscale diffusions[END_REF][START_REF] Pavliotis | Parameter estimation for multiscale diffusions[END_REF].

Consequently, a time down-sampling scheme is proposed to force the noise terms to be as decorrelated as possible.

By assuming that the spatially averaged covariance function has a Gaussian form with a standard deviation equals to the correlation time τ , a simple expression allows us to compute it. For a given unresolved velocity correlation matrix we write:

C v ij = v obs (•, t i ), v obs (•, t j ) = C v ij - n k=1 b obs k (t i )b obs k (t j ), 0 i, j N -1, (4.18)
and its associated stationary covariance function

Cov s (t p ) = 1 N -p N -1-p q=0 C v q,q+p , 0 p N -1. (4.19)
We propose the following correlation time estimation:

τ = 2 Cov 2 s ∆Covs ∆t 2 , (4.20)
using a forward Euler temporal discretization of the stationary covariance:

∆Cov s ∆t (t p ) = Cov s (t p + ∆t) -Cov s (t p ) ∆t , 0 p N -1. (4.21)
These estimations follow over-simplified assumptions known to have a restricted validity.

Experimentally, they still systematically provide the best simulation results when compared to other more complex estimators (whose derivation is outside the scope of this paper).

Moreover, with a white unresolved velocity with a dirac stationary covariance function, we obtain τ = dt. This result is exactly what we would expect from a white unresolved velocity, to prevent an overly aggressive down-sampling.

Before computing the estimations presented in sections 4.2, 4.3 and 4.4, we use this estimated correlation time τ to down-sample both the entire dataset and the observed coefficients of the reduced order solution, leaving us with a time step ∆t ≈ τ . On top of the considerable gain in accuracy, this time down-sampling reduces the amount of data to process, and hence, makes the off-line ROM building process faster.

Having estimated and specified all of the stochastic ROM's parameters in equations (3.3) and (3.4), it is now possible to forecast ensembles of realizations of the ROM through Monte-Carlo simulations. Before presenting numerical results, we discuss some important properties of LU models and of their reduced order versions.

5. Kinetic energy budget. First, the conservative properties of the LU Navier-Stokes representation are recalled, followed by a proof that by combining Galerkin projections and the advecting velocity correction compressibility, an intrinsic energy dissipation appears.

Full-order model budget.

As derived in [START_REF] Resseguier | Geophysical flows under location uncertainty, part I random transport and general models[END_REF][START_REF] Resseguier | Data-driven versus self-similar parameterizations for stochastic advection by lie transport and location uncertainty[END_REF], since the pressure does not influence the energy budget, by neglecting the molecular viscosity, the divergence of w and boundary conditions effects, and by applying Itō's lemma, the expression for the kinetic energy budget writes:

d dt 1 2 w 2 L 2 = -1 2 G * ( G(w)) , w
Loss by diffusion

η 2 HS = Ω 2 d i,j=1 η2 i,j (x, z)dxdz (5.2)
is the squared Hilbert-Schmidt norm of the integral operator η. This enables us to state that the energy is conserved for each realization of the stochastic process, and as a direct consequence of this, the energy intake of the noise and the dissipation by the turbulent diffusion must exactly compensate each other. The latter dissipates the kinetic energy of the mean E{w} 2 L 2 , while the former only releases random energy w -E{w} 2 L 2 into the system. Thus, the time-uncorrelated component of the velocity field drains energy from the mean field to the random component of w.

We can also express this energy transfer with the expectation of equation (5.1):

d dt Ω Var(w) = d dt E w -E{w} 2 L 2 = -d dt E{w} 2 L 2 ∀t ∈ [0, T ]. (5.3)
Besides the physical relevance of energy conservation, variance inflation and its relation to the mean field are also of primary interest for data assimilation or ensemble forecasting issues.

Reduced

d t w R = n i=0 db i φ i = n i=1 (dM R i )(b) φ i = Π φ P(dM)(w R ) ∀t ∈ [0, T ], (5.4)
where Π φ is the projection onto the reduced subspace. Specifically, for any function ζ ∈ L 2 , the projection Π φ is defined as

Π φ [ζ] = n p=1 φ p , ζ φ p .
As previously stated, we can evaluate the variation of kinetic energy with the Itō formula. The following result is proved in Appendix B for every t ∈ [0, T ]:

d dt 1 2 w R 2 L 2 = -1 2 G(w R ) 2 HS
Loss by diffusion

+ 1 2 Π φ (G)(w R ) 2 HS
Energy flux from the noise

-1 2 ∇ • (∇ • a) T , w R 2
Advecting velocity correction compressibility , (5.5) The second term of (5.6) informs about the energy flux due to mode truncation. The ROM's subgrid diffusion extracts energy from some temporal coefficients of the reduced order solution, while the multiplicative noise distributes it to others. These stabilizing and destabilizing effects recreate a large part of the energy fluxes between coefficients of the reduced solution, otherwise lost in deterministic ROMs. Thus, the uncorrelated velocity component drains energy from the coefficients of the reduced solution to give it back to the temporal coefficients of the truncated modes -i.e. modes orthogonal to the reduced space. This energy flux is exactly characterized by the second term of (5.6). But, our stochastic ROM cannot transfer energy from the truncated modes to the reduced solution's coefficients, as the former cannot be specified due to the truncation.

= -1 2 ∇ • (∇ • a) T , w R 2 -1 2 Π ⊥ φ (G)(w R )
Such a dissipation could be perfectly prevented by considering

-1 2 (Π φ [G]) * ( Π φ [G])
represented in the ROM by - Changing the ROM notations from Itō to Stratonovich [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF], the following correction term appears:

1 2 d b T α •p• , B = n q=0 1 2 n i=1 Σ α ip,qi b q dt = n q=0 -1 2 n i=1 Σ α pi,qi b q dt. (5.7)
In Stratonovich form, the ROM also still implicitly includes the diffusion term ensuing from the truncature. Notwithstanding, this additional dissipation in (5.6) is not a problem as it enables the setting of an energy dissipation mechanism representative of a direct energy cascade in the truncated modes' coefficients. Hence, for very turbulent flows described by few modes, the energy leak appears to be necessary as our main concern is to restore energy fluxes between the coefficients of the reduced order system. Indeed, those energy fluxes can be very difficult to model correctly in reduced versions of non-linear systems. As an example, [START_REF] Sapsis | Blending modified Gaussian closure and non-Gaussian reduced subspace methods for turbulent dynamical systems[END_REF] observe that a basis encoding 50% of the energy can lead to a ROM missing more than 98% of energy transfers. These missing energy fluxes create instabilities in some coefficients of the reduced solution (missing negative energy fluxes) while overdamping others (missing positive energy fluxes). Restoring the energy fluxes between coefficients of the reduced system is a significant challenge LU models tackle through the combination of multiplicative noise and turbulent diffusion.

To note, the noise dimensionality reduction of section 4.4 introduces an additional energy leak due to the fact that we are only keeping n eigenvectors over n(n + 1) in the noise covariance matrix's PCA. This introduces another negative energy flux, as the noise variance is reduced while the dissipation is maintained.

6. Ensemble forecasting. As a means to measure the performance of the proposed ROM in UQ scenarios, ensemble simulations were carried out using the data available from DNS simulations of wake flows at Reynolds 100 and 300, the former being quasi two dimensional while the latter, fully three dimensional. Wake flows are well-studied, non-linear, oscillatory flows which are physically produced by a uniform-velocity inflow facing a solid obstacle -here, a cylinder. Vortices are thus created behind the obstacle and periodically detach from it.

The dynamics, and in particular their intrinsic dimensionality, strongly differ between Reynolds 100 and 300. Indeed, as illustrated by figure 1, at Reynolds 100, most of the energy is concentrated in just a few coefficients of the reduced order solution. In contrast, at Reynolds 300, the two first modes are meaningful for a rough approximation, but energy spreads over many degrees of freedom. The top panels of figures 3 and 7 confirm a strong difference in complexity between the two flows. The ROM's results w R will be compared directly to DNS simulations v ref on a time interval outside intervals on which the ROM's coefficients and basis functions have been estimated.

6.1. Baseline ROMs. To better appreciate our stochastic ROM's capabilities, results are compared to two different state-of-the-art algorithms.

d dt b i = Re Re ev b T l •i + b T c ••i b, 1 i n. (6.1)
The eddy viscosity coefficient is typically fitted by least squares using the dataset b obs .

Without that correction, some coefficients of the reduced solution might become unstable because of the missing negative energy fluxes due to mode truncation mentioned in 5.2.

The eddy viscosity term generally manages to stabilize those coefficients, but this does not always make this data-driven method accurate. In order to emulate that unstable behavior, we also simulate a simple DNS POD-Galerkin with no eddy viscosity model, i.e. the ROM (6.1) with Re ev = Re.

Stochastic baseline ROM.

As a deterministic model, the ROM above can hardly be used for dynamics-error UQ nor ensemble-based data assimilation. As a matter of fact, a simple randomization of the initial conditions are known to usually lead to error underestimation (see section 1). So, to keep things simple, we define a second baseline ROM by adding a white noise term to the first baseline ROM. This constitutes a simple adhoc randomization technique for a given dynamical system through a Gaussian additive forcing. Despite its potential lack of physical relevance, such a strategy is very often adopted in data-assimilation applications [START_REF] Chapron | Large-scale flows under location uncertainty: a consistent stochastic framework[END_REF].

db i = Re Re ev b T l •i + b T c ••i b dt + σ ev i• dW t , 1 i n, (6.2)
where σ ev ∈ R n×n is the Cholesky decomposition of the ROM's noise covariance Σ ev = σ ev (σ ev ) T and t → W t ∈ R n is a vector of n independent Brownian motions. Since no prior physical information is available for these baseline models, the noise's covariance is also learned from data as follows: perfectly the reference's phase, and exhibits a slightly damped magnitude. The damping effect can be seen as a consequence of the exchange of energy between the mean and the variance (see equation (5.3)). This interplay is also evidenced by the growth of the confidence interval. The realization we have (randomly) singled-out shows a slowly divergent behavior.

Σ ev ij = ∆t

Velocity forecasts. Having analyzed the model's forecast capabilities in terms

of the ROM's temporal coefficient, we now focus on the analysis of the velocity fields' forecast and compare them to the reference -i.e. the projection of full DNS solution onto the POD basis :

Π φ [v ref ] = Π φ [v ref -v] + v = n i=0 b ref i φ i (with b ref 0 = 1).
Here, our goal will be to assess qualitatively the potential limitations of the ROM. Through equation (3.1), we will compute the velocity fields for the mean, the lowest error realization, and the prediction generated by the eddy viscosity ROM (6.1) of same dimension. At each time step, we define the lowest-error realization as follows:

b min = argmin (b (k) ) k w R -Π φ [v ref ] 2 L 2 = argmin (b (k) ) k b (k) -b ref 2 . (6.6)
We plot in figure 3 the vorticity fields of the ROMs of dimension n = 8. A system of dimension 8 enables us to reproduce fairly well the topology of the velocity fields, and from figure 2 we conclude that our model's best realization is capable of staying in phase with the reference, even 50 seconds into the validation set, while the ensemble realization's mean and the deterministic baseline start losing accuracy. Indeed, the eddies of the latter predictions are late compared to the reference. This misalignment of eddies is a major type of model error in CFD, and it is very problematic for data assimilation because misalignment can hardly be corrected by linear Eulerian algorithms [START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF]. The LU model's cornerstone is precisely the random locations of transported quantities, and in particular, the locations of eddies. This random location idea is encoded in the random flow equation (1.2) and gave rise to the name "Location Uncertainty" (LU). Here, the LU POD-ROM generates an ensemble of simulations. Some of them have eddies in advance compared to the mean and some have late eddies. This variety of eddies' locations allows some ensemble realizations -in particular the best one -to match the reference's eddies' locations.

Global prediction accuracy.

To compare more precisely the performance of the proposed ROM, we systematically plot the RMSE, as well as the ensemble's bias, variance and minimum error. These quantities will also be compared to those obtained by the stochastic baseline ROM (6.2). This enables us to inspect the UQ capabilities of the proposed ROM.

Although the random energy transfers are a necessary feature when the number of modes is not large enough, it increases the temporal coefficients' variance while keeping their biases constant, which yields an increase of the Root Mean Square Error (RMSE) of the whole ensemble:

RMSE = Ê w R -Π φ [v ref ] 2 L 2 = Ê{w R } -Π φ [v ref ] 2 L 2 + Ω Var(w R ). (6.7)
In all the following normalized error plots, the blue line corresponds to the solution computed using a deterministic POD-Galerkin ROMs (with and without eddy viscosity) and the red one, to the RMSE. The green line represents the ensemble bias, whereas the magenta line is the error incurred by the ensemble's solution closest to the reference (6.6), and the gray shade corresponds to 1.96 × the standard-deviation. These quantities are computed for the whole set of generated realizations. The initial condition is common to all and the values are normalized by the square root of the solution's energy averaged over the training set : (6.3). This large additive noise hence quickly fills the reduced subspace overwhelming the physical terms of the reduced dynamics. As for the n = 8 simulation, we observe that the LU POD-ROM differs from the other two simulations by having unstable realizations, thence the sudden peaks in the bias and mean curves. When a realisation diverges, we re-sample it uniformly from the other members of the ensemble.

v 2 L 2 + n i=1 λ i 1/2
Focusing on our stochastic POD-ROM, for any n ∈ {2, 4, 8}, even though the ensemble bias grows as time passes by, the best realization error does not. This is a crucial property for a UQ model to have, as this means that an appropriate filtering technique could eventually retrieve it at each observation, leading to a stable, low error data assimilation system.

6.3. 3D wake flow at Reynolds 300. For a more challenging test, a 3D wake flow behind a circular cylinder at Re = 300 was simulated, constructing the ROM just like in the case of the wake flow at Reynolds 100. The 440 seconds of DNS simulation data (about 88 vortex shedding pseudo-periods) were split into training (80 vortex shedding pseudo-periods) and test sets (8 pseudo-periods). Similarly, the ROM was subjected to the same benchmarks as before, with the equivalent eddy-viscosity UQ-enabled model for comparison. As it is a more complex flow, we expect the different temporal coefficients to be harmonically richer and to capture as much of the small-scale interactions as possible, but yet not to be able to fully reproduce the intricacies of the DNS-simulated flow. 6.3.1. Temporal coefficients forecast. This new scenario being fully 3D and more complex, the advantages of the proposed model are clear. This can be seen from the temporal coefficient forecasts, where the proposed LU stochastic ROM manages to stabilize the system while the deterministic PODs starts diverging after 5 seconds (i.e. before one complete vortex shedding), as evidenced by the first temporal coefficients in figure 5. Intermittency probably fools the learning procedure of the eddy viscosity method and makes it less robust and not adapted to the test set. In contrast, the LU learning procedure -section 4 -shows a greater degree of robustness and leads to more accurate results. It is also interesting to note how the random energy contribution increases with the order of the temporal coefficient in the form of ensemble variance. This is attributed to energy transfers between temporal coefficients, facilitated by the interplay of multiplicative noise and turbulent diffusion. This effect is in fact amplified in the 8-dimension ROM in figure 6, where the temporal coefficients' means get more damped as the amount of random energy in the system increases. This energy dissipation mechanism helps the system to remain bounded, unlike its deterministic counterparts.

Interestingly enough, in the n = 8 dimension ROM, even when the ensemble variance increases with the temporal coefficient order, coefficients 5, 6, 7 and 8 of the random realization (yellow line) also attempt to follow the right reference amplitude and variability, as it can be observed in figure 6.

6.3.2. Velocity forecasts. As in the previous case, we now analyze qualitatively the performance of our proposed method in terms of velocity field predictions.

As can be observed from figure 7, the best prediction remains stable and manages to stay close to the theoretically optimal solution in 8-dimensions on a period of four shedding cycles after the learning period. Conversely, the baseline model starts to diverge, as evidenced by the red zones inside the vortices. This result is quite impressive as, even if the mean velocity field is not as close to the reference as one would like, it still reveals the potential of a method combining our model with a data assimilation technique to retrieve the best realization and provide corrections on the fly. is applied at a rate equal to the residual's correlation time. Finally, in order to restrict the number of coefficients needed to characterize the ROM to the usual O(n 3 ), a technique to reduce the dimensionality of the noise covariance matrix is also proposed.

The conservative properties of the LU closure are discussed, and we demonstrate that, when performing Galerkin projections of the Itō form of LU SPDEs, we end up forfeiting these properties. In the reduced order version, two energy fluxes appear: the first one is attributed to the possible advecting velocity correction divergence, whilst the second one is negative and is directly associated with mode truncation. We argue that this energy loss is understandable, and even desirable, with an interplay between noise and turbulent diffusion to maintain non-linear energy fluxes between coefficients of the reduced solution despite mode truncation.

Numerical comparisons are performed between our stochastic ROM and state-of-the-art deterministic and stochastic ROMs. As test cases, we chose a two-dimensional wake flow at Reynolds 100 with few degrees of freedom, and a more complex three-dimensional wake flow at Reynolds 300 with many more degrees of freedom. Deterministic Navier-Stokes simulations stood for references and all of the ROMs were initialized with their values, and the forecasts, compared to them. For the stochastic ROMs, ensembles of 100 realizations were simulated. The state-of-the-art stochastic ROM's solutions quickly diverged in time, whilst the temporal coefficients of our ROM proved to be neither unstable nor over-damped. The LU POD-ROM solution's biases were even found to be smaller than those of each of the other stochastic ROMs' that we have considered. Moreover, at any given moment, the LU POD-ROM ensembles managed to remain very close to the reference, suggesting that they could be efficient priors for Bayesian inverse problems.

Following this work, we have applied the LU POD-ROM with unknown initial conditions in conjunction with particle filter algorithms to estimate velocity flows in real-time from few local measurements [START_REF] Resseguier | Real-time algorithm to estimate and predict unsteady flows from few measurements using reduced-order models[END_REF].

Finally, when working at very large Reynolds numbers, DNS simulations are no longer an option, but LES (Large Eddy Simulation), DDES (Delayed Detached Eddy Simulation)

and RANS (Reynolds Averaged Navier Stokes) can still provide useful data to build ROMs, at the cost of some approximation error. To tackle this, the closure mechanism's estimation procedure may be adjusted to address the small-scale velocity statistics neglected by LES-like approaches. Another approach that we are currently exploring is to consider the Galerkin projection on a LU version of LES and DDES. Non-polynomial terms of LES and DDES complicate the dimensional reduction, but discrete empirical interpolation methods (DEIM) [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF] may be able to handle it. As illustrated in this paper, our stochastic closure naturally prevents truncation-induced instabilities, and we expect a similar behavior at higher Reynolds. If the proposed in-house stabilization is not sufficient, existing stabilisation and/or data-driven methods [START_REF] Georgaka | A hybrid reduced order method for modelling turbulent heat transfer problems[END_REF][START_REF] Hijazi | Data-driven pod-galerkin reduced order model for turbulent flows[END_REF][START_REF] Torlo | Stabilized weighted reduced basis methods for parametrized advection dominated problems with random inputs[END_REF]] could be added. Furthermore, CFD outputs may also be at a resolution different from the resolution of the measurements. Hence, depending on the relative precision of CFD and measurements, there could be an additional or a missing very-small-scale residual turbulence in the assimilated measures. Accordingly, we will need a simple statistical model for this very-small-scale turbulence and include its effects in the observation model through either an additional "measurement noise" or additional smoothing.

Besides, in realistic applications, Reynolds number, initial and boundary conditions, as well as a plethora of other parameters, are often only approximately known. So, combinations with other state-of-the-art dimensionality reduction methods are most likely necessary. Typically, we could use wRB [START_REF] Chen | A weighted reduced basis method for elliptic partial differential equations with random input data[END_REF][START_REF] Chen | Reduced basis methods for uncertainty quantification[END_REF][START_REF] Torlo | Stabilized weighted reduced basis methods for parametrized advection dominated problems with random inputs[END_REF] for uncertain mean inflow (Reynolds number and angle of attack) and a reduced pressure equation [START_REF] Stabile | POD-Galerkin reduced order methods for CFD using finite volume discretisation: vortex shedding around a circular cylinder[END_REF][START_REF] Stabile | Finite volume POD-Galerkin stabilised reduced order methods for the parametrised incompressible navier-stokes equations[END_REF] for time dependent boundary conditions. Our new implementation of LU terms in the ITHACA-FV library [START_REF] Stabile | POD-Galerkin reduced order methods for CFD using finite volume discretisation: vortex shedding around a circular cylinder[END_REF][START_REF] Stabile | Finite volume POD-Galerkin stabilised reduced order methods for the parametrised incompressible navier-stokes equations[END_REF] will surely aid on these practical aspects.

Nevertheless, note that a wRB method to simulate random velocity effect [START_REF] Torlo | Stabilized weighted reduced basis methods for parametrized advection dominated problems with random inputs[END_REF] does not directly apply to the purpose of the current paper. Indeed, we try to quantify the uncertainty induced by the dimensional reduction of the deterministic Navier-Stokes equations and not to reduce the dimension of an uncertainty quantification problem. The random residual velocity can hardly be considered a wRB parameter because we do not choose its value. There is only one "realization" and this is the one given by the deterministic Navier-Stokes equations. We could specify some statistics for the random residual velocity and simulate the stochastic Navier-Stokes equations (2.1)-(2.2) but this would make our method much more complex, much more CPU demanding, and could possibly lead to over-fitting. Moreover, that random velocity is a residual, and thus, its statistics depend on the reduced dimension n. That would also greatly complicate the wRB algorithm.

Finally, relying purely on data can offer the advantage of directly tackling very high Reynolds number configurations. It would require us to filter noisy measurements, and most likely, to rely on simplified (often 2D) data-driven models [START_REF] Artana | Strong and weak constraint variational assimilations for reduced order fluid flow modeling[END_REF][START_REF] Cammilleri | POD-spectral decomposition for fluid flow analysis and model reduction[END_REF][START_REF] Fick | A stabilized pod model for turbulent flows over a range of reynolds numbers: Optimal parameter sampling and constrained projection[END_REF][START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]. From this latter point of view, modeling under location location offers a great flexibility which might be of great interest. As a matter of fact, the LU formalism enables us to derive random dynamics close to the data but still based on stochastic transport principles [START_REF] Avenel | Stochastic level set dynamics to track closed curves through image data[END_REF].

The noise parameters can be updated online from incoming data and simulations. A full characterization of the noise dynamics through statistical learning procedures would then probably enable us to devise advanced data-driven models that are able to handle even more complex situations.

Appendix A. Estimation formulas. In this appendix, we consider a full probability space (A, F, P ), and a filtration of σ-algebra {F t } t≥0 . For the Hilbert space H 2 (Ω), we understand H 2 = H 2 (0, T ; Ω, F) as the space of strongly measurable, {F t } t≥0 adapted processes u : R d × [0, T ] → R d such that (x → u(x, t)) belong to H 2 (Ω) for every t ∈ [0, T ] and

u 2 H 2 = E T 0 u(•, t) 2 H 2 (Ω) dt < ∞. (A.1)
Similarly, we can define L 2 (0, T ; Ω, F) from L 2 and its norm. We recall below a classical proposition related to quadratic variation process: Proposition: Stochastic integration and quadratic variations. If M is a continuous martingale and X ∈ L 2 (0, T ; Ω, F), then there exists a unique bounded continuous martingale As stated in the introduction, we assume that the operator σ is Hilbert-Schmidt and that its kernel σ : Ω 2 → R d×d is C 2 over the compact bounded set Ω 2 . Therefore, As such, ((x, t) → σ(x)B t ) belong to H 2 . In particular, (x → σ(x)B t ) belong to H 2 (Ω) for every t ∈ [0, T ] and ∂ x i σ, ∂ x i ∂ x j σ are also Hilbert-Schmidt for every 1 i, j d.

Additionally, we assume that the POD modes φ i belong to H 2 (Ω) and have bounded gradients.

After defining boundary conditions, P and K are well defined on H 1 (Ω) and H 2 (Ω) respectively. Then, for every ζ ∈ H 1 (Ω), Pζ L 2 ζ L 2 since P is a projection and then, for every ξ ∈ H 2 (Ω), by Cauchy-Schwarz and triangular inequalities: Then, upon applying the Itō formula to the local kinetic energy we obtain

|K jq (ξ)| φ j L 2 ( ∇φ T q ∞ + ν) ξ H 2 (Ω) , (A.
d 1 2 w R 2 L 2 = Ω d t (w R ) T w R + 1 2 d t (w R ) T , w R ∀t ∈ [0, T ]. (B.3)
To remove the (orthogonal) projection operator Π φ from the first term, we exploit its symmetry, and afterwards, the fact that w R is already in the reduced subspace: From the definition of the projection operator Π φ and using the extended notation (η θ) = Ω η(z) θ T (z)dz, we can express the above quadratic operator with the noise statistics as follows:

Ω Π φ G(φ p ) Π φ G(φ q ) = n i=1 φ i , G(φ p ) φ i , G(φ q ) = n i=1
Σ α pi,qi . (B.12) Finally, by orthogonality, the kinetic energy budget (B.3) simplifies to:

d dt 1 2 w R 2 L 2 = 1 2 G(w R ) 2 HS -1 2 ∇ • (∇ • a) T , w R 2 -1 2 Π φ G(w R ) 2 HS , (B.13) = -1 2 Π ⊥ φ G(w R ) 2 HS -1 2 ∇ • (∇ • a) T , w R 2 , (B.14)
where Π ⊥ φ = I d -Π φ is the projector into the orthogonal complement of the reduced subspace.

2 .

 2 Navier-Stokes model under location uncertainty. Galerkin projections of Navier-Stokes equations do not specifically take into account residual velocity contributions, and thus, nor do they offer a precise quantification of their induced errors. To help remedy this issue, we propose to directly project SPDEs -i.e. the LU Navier-Stokes representationinstead of the classical Navier-Stokes equations. The resulting ROM is expected to describe the same physical system, as we postulate that the solution of the SPDE to be statistically similar to a large-scale representation of the original (deterministic) equation's solution. 2.1. The random physical model. Let us denote by Ω an open bounded subset of R d ,

D

  t involves the usual terms of the deterministic material derivative, on top of three additional new terms: an advecting velocity correction (w * instead of w), a heterogeneous and anisotropic turbulent diffusion, and a multiplicative noise. This last term corresponds to the advection by the unresolved velocity σ Ḃ. Finally, we can recover the classical Navier-Stokes equations by setting the residual velocity to zero -i.e. σ = 0. Let us highlight that for conserved scalars, this operator corresponds to the material derivative -i.e. the derivative along the stochastic flow d (w(X t , t))[START_REF] Resseguier | Geophysical flows under location uncertainty, part I random transport and general models[END_REF]. This stochastic Navier-Stokes model

(4. 3 ) 4 . 1 .

 341 The curse of dimensionality. Computing correlations of Gaussian noises ( αdB t )

λ obs p K

 p jq b obs p ∆b obs i ∆t v obs , 1 i, j n, 0 p, q n, (4.10) where b obs 0 = λ obs 0 = 1 and for 1 i n, b obs i = (φ i , v obs ) ,

  the part inside K jq only involves a projection of the residual velocity along the time dimension. Finally, the n(n + 1) differential operators K jq are computed only on n(n + 1) functions that do not depend on the time variable. We prove the consistency of this estimator in Appendix A by capitalising on the quadratic covariation's definition and the orthogonality of the observed temporal coefficients b obs i . In addition, this can also be extended to non-orthogonal coefficients by solving the linear system engendered by the matrix b p b k pk . 4.4. Noise dimension reduction. Using the Cholesky decomposition σ α of the noise covariance tensor:

  )where αR k k ∈ R (n+1)×n are the matrix forms of the first n eigenvectors (weighted by the corresponding eigenvalues' square roots) and (β (k) ) k are n independent one-dimensional Brownian motions. Since αdB t is a multiplicative noise and the temporal coefficients b i have various amplitudes √ λ i , the covariance matrix Σ α is adequately re-normalized by the amplitudes √ λ i before applying the PCA.

4. 5 .

 5 Time down-sampling rate. Under the LU Navier-Stokes model hypothesis, the unresolved term of the velocity field corresponds to a noise uncorrelated in time.

ξ 2 L 2 =

 22 (ξ, ξ) is the squared norm of L 2 and

  -order model budget. Following (5.1), the full-order LU Navier-Stokes model (2.1)-(2.2) conserves the kinetic energy, up to molecular viscosity and boundary condition effects. With reduced order models, the advecting velocity correction is expected however to create either energy compression or dilation, and the mode truncation, to introduce a small energy leak. The ROM (3.3) does not exactly solve the global LU Navier-Stokes model (2.1)-(2.2) but the Galerkin projection of the divergence-free component of equation (2.1) :

  • stands for the Euclidean norm on R d and Π ⊥ φ = I d -Π φ is the projector onto the orthogonal complement of the space spanned by the reduced basis' functions. The first term of (5.6) corresponds to the fluid compression/dilation created by the velocity correction w * -w = -1 2 (∇ • a) T . By construction, the ROM meets the incompressible condition -i.e. ∇ • w = 0 -but unlike the LU Navier-Stokes model, it does not satisfy the finite-variation part of the mass conservation (2.2): ∇ • w * = 0. The velocity correction divergence is associated to turbulence heterogeneity [6, 16, 67], corresponding to the spatial maxima (positive energy fluxes) and minima (negative energy fluxes) of the turbulence's kinetic energy v 2 .

6. 1 . 1 .

 11 Deterministic baseline ROM. The first state-of-art ROM is a prototype of deterministic ROM widely used in fluid dynamics. After a POD-Galerkin on the classical Navier-Stokes equations, the molecular viscosity coefficient 1/Re is replaced by an eddy viscosity coefficient 1/Re ev 1/Re [15]:

  ev b obs T l •i + b obs T c ••i b obs . (6.5)As for the LU ROM, the implicit assumption of white noise residual∆b obs i ∆t evmay not be valid. Thus, similarly to section 4.5, an optimal down-sampling time step is estimated withn k=1 ∆b obs k ∆t ev (t i ) ∆b obs k ∆t ev (t j ) instead of C v ij, before the noise covariance estimation.

Figure 1 :

 1 Figure 1: Energies of the reduced solution coefficients normalized by the total solution energy temporal mean v 2 L 2 for the 2D flow at Reynolds 100 (left panel) and the 3D flow at Reynolds 300 (right panel).

6. 2 .

 2 2D wake flow at Reynolds 100. Starting with the simpler of the two test cases, 800 seconds of DNS simulation data (about 160 pseudo-periods) were generated for a 2D wake flow behind a circular cylinder at Re = 100. These data were then split into training and test sets, the former (140 pseudo-periods) being utilized for the construction of the ROMs, and the latter (20 pseudo-periods), as a reference to which the ROMs will be compared. Performance will thus be measured through four different metrics: comparison of temporal coefficients in the time domain, comparison of reconstructed velocity fields in the spatial domain, global prediction accuracy as measured by the evolution in time of RMSE, and finally ensemble minimum error. 6.2.1. Temporal coefficients forecast. The 2-dimensional LU ROM of the Reynolds-100 flow is simulated 100 times. In figure 2, the coefficients of the reduced solution b (k) i (where k stands for the k-th simulation) are compared to the coefficients of the reference solution b ref i = (φ i , v ref ) (plotted in black). The ensemble mean (in green) follows almost

Figure 2 :

 2 Figure 2: LU POD-ROM forecast for n = 2 coefficients of the reduced solution for a 2D wake flow at Reynolds 100: ensemble mean (green line), one random realization (yellow line), confidence interval (gray shade). Blue lines correspond to deterministic ROMs: baseline eddy-viscosity POD-ROM (6.1) (light-blue dashed line) and DNS POD-Galerkin (dark-blue solid line). The dash-dot black plots represent the observed references.

Figure 3 :Figure 4 :

 34 Figure 3: Vorticity fields -10 vortex shedding cycles after the learning period -from (from top to bottom) the reference simulation v ref (2D DNS at Reynolds 100: state space dimension of about 10 4 ), its projection onto the POD basis Π φ [v ref ], the LU POD-ROM model's best prediction, the LU POD-ROM model's ensemble mean, and the prediction of the deterministic baseline POD-ROM (6.1) (fitted eddy viscosity) (ROM state spaces are of dimension 8). The vorticity is the velocity curl and is the usual and a convenient way of visualizing 2D flows and their vortices.

Figure 5 :

 5 Figure 5: LU POD-ROM forecast for n = 4 coefficients of the reduced solution of a 3D wake flow at Reynolds 300: ensemble mean (green line), one random realization (yellow line), confidence interval (gray shade). Blue lines correspond to deterministic ROMs: baseline eddy-viscosity POD-ROM (6.1) (light-blue dashed line) and DNS POD-Galerkin (dark-blue solid line). The dash-dot black plots represent the observed references.

Figure 6 : 7 .

 67 Figure 6: LU POD-ROM forecast for n = 8 coefficients of the reduced solution of a 3D wake flow at Reynolds 300: ensemble mean (green line), one random realization (yellow line), confidence interval (gray shade). Blue lines correspond to deterministic ROMs: baseline eddy-viscosity POD-ROM (6.1) (light-blue dashed line) and DNS POD-Galerkin (dark-blue solid line). The dash-dot black plots represent the observed references.

Figure 7 :

 7 Figure 7: Q-criterion iso-surfaces -4 vortex shedding cycles after the learning periodfrom (from left to right and from top to bottom) the reference simulation v ref (3D DNS at Reynolds 300: state space dimension of about 10 7 ), its projection onto the POD basis Π φ [v ref ] , the LU POD-ROM model's ensemble mean, the LU POD-ROM model's best prediction, and the prediction of the deterministic baseline POD-ROM (6.1) (fitted eddy viscosity) (POD-ROM state spaces are of dimension 8). The Q-criterion [39] is a quadratic function of the velocity gradient and is the usual and a convenient way for visualizing 3D flows and their vortices.

Figure 8 :

 8 Figure 8: Normalized error for n = 2, 4 and 8 coefficients of the reduced solution of a 3D wake flow at Reynolds 300 -with the projection of the DNS onto the POD basis Π φ [v ref ] as reference : RMSE (red line), bias (green line), ensemble minimum error (magenta line) and 1.96 × the standard-deviation (shaded gray) for the LU POD-ROM (left panel) and stochastic baseline POD-ROM (right panel) . Blue lines correspond to deterministic ROMs: baseline eddy-viscosity POD-ROM (6.1) (light-blue dashed line) and DNS POD-Galerkin (dark-blue solid line). The black solid line at the top is the error considering only the time mean velocity, i.e. b i = 0, ∀i > 0.

t0

  XdM such that for every continuous martingale N (with zero initial condition)

2 σ 2 + ∂ x i σ 2 + ∂ x i ∂ x j σ 2 Ω 2 T

 22222 < ∞, (A.3) where • stands for the Euclidean norm on R d×d and |Ω 2 | the Lebesgue measure of Ω 2 .

< 0 b p d b i , t 0 ( 0 (b

 000 ∞ . The bound above gives the continuity of K on H 2 (Ω). We also assume that the observed coefficients of the reduced solution b obs i = (φ i , v obs ) are continuous semi-martingales and solutions of the ROM (3.3)-(3.4). From there, the orthogonality of the coefficients of the reduced solution yields, for 1 i, j n and 0 p, q n:T αqj• dB s ) = αki• dB s ) b k , p b k Σ α ki,qj , (A.6) = T λ p Σ α pi,qj . (A.7)

Ω 1 2Π

 1 d t (w R ) T w R = Π φ (dM)(w R ) , w R , (B.4) = (dM)(w R ), Π φ w R , (B.5) = (dM)(w R ), w R , (B.6) = ((F dif + F adv )(w R ) + C(w R , w R ))dt + (GdB t )(w R ), w R , (B.7) = -1 2 G * G(w R ) , w R dt from (F dif (w R ),w R ) using (2.12) + 1 2 ((∇ • a)∇)w R , w R from (F adv (w R ),w R ) + C(w R , w R )dt + (GdB t )(w R ), w R =0 by skew-symmetry of ξ →C(w R ,ξ) and GdBt , (B.8) = 1 2 G(w R ) 2 HS -1 2 ∇ • (∇ • a) T , w R 2 , (B.9)where the second term comes from integration by parts. Besides, the Itō term of the energy budget is straightforward to compute from equation (B.2):Ω d t (w R ) T , w R = 1 2 Π φ (G)(w R ) φ G(φ p ) Π φ G(φ q ) b p b q dt. (B.[START_REF] Brzeźniak | Stochastic partial differential equations and turbulence[END_REF]

Table 1 :

 1 Terms of LU Navier-Stokes equations and their algebraic properties closure mechanism and its reduced order expression. The stochastic transport operator

	Advection by the	GdB t	-(σdB t • ∇)	skew-symmetric	αdB t
	residual velocity			noise
				(if ∇ • σ = 0)
	Molecular viscous	HdB t	1 Re ∆σdB t	additive	θdB t
	dissipation of the			noise
	residual velocity
	Pressure	dP	-∇ (p dt + p σ dB t )	Potential	Included in
	gradient			field	other terms

  G(φ q ) p,q . This correction would have been implicit using the Galerkin projection with Stratonovich calculus. As demonstrated in[START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF][START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF][START_REF] Resseguier | Data-driven versus self-similar parameterizations for stochastic advection by lie transport and location uncertainty[END_REF], using Stratonovich integral, the turbulent dissipation term does not explicitly appear, neither when deriving the ROM from scratch nor when switching notations after the Galerkin projection. Let us stress that the modified advection is still present in the Stratonovich form of the LU model.

	1 2	n i=1 Σ α pi,qi p,q	(see equations (B.11) and (B.12) in Ap-
	pendix B), instead of the reduced version of the full turbulent diffusion operator
	-1 2 φ p , G

*

.

  Reference simulation v ref (2D DNS at Reynolds 100: state space dimension of about 10 4 ) Projection of the reference onto the POD basis Π φ [v ref ]

	(State space
	of dimension 8)
	LU POD-ROM model
	best prediction
	(ROM state space
	of dimension 8)
	LU POD-ROM model
	ensemble mean
	(ROM state space
	of dimension 8)
	Prediction of the deterministic
	baseline POD-ROM (6.1)
	(fitted eddy viscosity)
	(ROM state space
	of dimension 8)

which is coherent with the fact that ∇ • v = 0 from (3.1) because ∇ • v= 0.

Acknowledgments. We warmly thank Pranav Chandramouli for the generation of the three-dimensional wake flow data, Pierre-Louis Lee for his help in the generation of the two-dimensional wake flow data, and Reda Bouaida for his help in the computation of the ROM matrix f . The authors also acknowledge the support of the ERC EU project 856408-STUOD, the ESA DUE GlobCurrent project, the "Laboratoires d'Excellence" CominLabs, Lebesgue and Mer through the SEACS project. Finally, we also thank Darryl D. Holm, Dan Crisan, Wei Pan and Igor Shevchenko for the insightful discussions.

References.

This work was supported by the ERC EU project 856408-STUOD, the ESA DUE GlobCurrent project, the "Laboratoires d'Excellence" CominLabs, Lebesgue and Mer through the SEACS project.

Now, let us note that for 1 j n and 0 q n: αqj• dB s = φ j , -P (σdB s • ∇) φ q + δ q0 ν∆σdB s = d (K jq [σB]) . (A. [START_REF] Berner | Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations[END_REF] Then, with the definition of the quadratic covariation and the increment notation ∆ξ(t k ) = ξ(t k+1 )-ξ(t k ), we obtain the estimator's expression and its consistency for every 1 i, j n and 0 p, q n:

where the continuity of the operator K jp on H 2 (Ω) enabled us to switch the limit in probability and the operator K jp . The martingale flow increments σ∆B t k are approximated by v (•, t k )∆t. In practice, for each 1 i n we can replace ∆b i by

where

Mathematically, this is still correct for very large values of T . Indeed,

0 and b i -b i has finite variations. Thus, b i -b i approaches a finite variation process (for large T ). Numerically, this formulation improves the accuracy of the model as it allows us to remove the smooth-in-time part of ∆b i , thus minimizing estimation error as well.

Appendix B. Energy dissipation.

On top of assumptions of Appendix A, we neglect the viscosity (L = 0) and assume zero Dirichlet boundary conditions for the reduced basis' functions and for x → σ(x)B t for every t ∈ [0, T ]. Note that these fields are divergence-free (since they are learned from a set of incompressible velocity fields). As such, we have:

also Hilbert-Schmidt.

Using the projected Navier-Stokes model (5.4), for every t ∈ [0, T ] we can formally remove the orthogonal projection by moving it into the divergence-free functions space P through integration by parts:


Denis Jabaudon
email: denis.jabaudon@unige.ch

Editorial overview: Developmental neuroscience

The past few years have seen an explosion in new technologies that have brought forward our understanding of the development of the central nervous system at the single-cell level. Singlecell RNA sequencing has shed a crisp light on the molecular diversity of cell types, while advanced 3D imaging techniques such as light-sheet microscopy now allow single axons to be identified and followed in intact cleared brains, bringing a new understanding on the assembly of cells into circuits and networks during development. Our genetic toolbox has been expanded by the development of new molecular tools allowing targeted ablation or silencing of single genes in defined cell types, and incremental progress in the generation and culture of 3D brainlike organoids and assembloids have informed our understanding of normal and abnormal developmental processes in humans. This volume of Current Opinion in Neurobiology is meant to exemplify some of these developmental facets, and contrast and unite them. Investigators at the forefront of their fields provide us with a state-of-the-art view of their topic of research, with at times related topics covered by multiple authors to provide a community point-of-view on complex issues.

How new technologies that can help understanding nervous system development are specifically covered in several articles. Shi and colleagues describe the methods and strategies to make human brain organoids and the role of self-organization (intrinsic) and extrinsic factors (growth factors, morphogens, extracellular matrix) in this process. Organoids can provide important insights into understanding human brain development, neurodevelopmental disorders and brain evolution. However, several challenges still lie ahead like the limited gyration, the low cellular diversity and the lack of vascularization. Le Bail and colleagues also tackle these questions in terms of temporal specification: how diverse are cortical neurons and what are the factors controlling the maintenance of cortical progenitors in vivo and in organoids. They next illustrate what these in vitro models have brought to our knowledge of human brain evolution and the unique features of human cortex development. Nedelec and Martinez-Arias also address these questions but taking the spinal cord as a model. After giving an overview of the molecular mechanisms that control spinal cord patterning and the specification of motor neurons and interneurons, they explain how this knowledge has been used to orient the in vitro differentiation of human pluripotent stem cells into spinal cord neurons, Again, generating the full spectrum of motor neuron diversity remains an important issue. Building locomotor circuits in a dish seems to be within reach which should help understanding the etiology of diseases such as ALS and developing new therapeutics for neuromuscular diseases. Also focusing on human brain organoids, Reiner and colleagues review on how these 3D structures can help us identify and understand human-specific processes of migration during development. emphasizing the value of such approaches in studying the pathophysiological bases of cortical malformations such as Miller-Dieker syndrome or lissencephaly. A final technology-oriented article is the one by Solanelles-Farre and Telley, who review the latest applications to the field of developmental neurobiology of multi-omics technologies such as single-cell RNA sequencing (scRNAseq) and single-cell assay for transposase (scATASseq) provide a compendium of developmental studies in this field to date. Cell profiling techniques, can unravel the transcriptome and epigenome of any cell type and significantly extended the cell repertoires of the brain and spinal cord. Lineage studies, using scRNAseq or genetic barcoding (somatic mutations) have also revealed unknow features of early brain patterning or cortical neuron specification.

In addition to these techniques-oriented reports, a variety of biology-centered reviews highlight multiple facets of brain development, from progenitor biology to neurogenesis, migration, axonal guidance and synaptogenesis.

The complex and emerging question of neuronal lineages in the developing cortex is covered by Lin and colleagues and Llorca and Marin, with a particular focus on how many types of progenitors are required to generate the diverse types of neurons present in the adult. An array of models has been proposed from strictly deterministic (i.e. single neuron types emerge from corresponding single progenitor types) to fully stochastic, whereby neurons emerge somewhat randomly (albeit with a temporal bias) from multicompetent progenitors as corticogenesis unfolds. Fischer and Morin also address the question of progenitor competence with a particular focus on "molecular recorders", which are molecular devices which can in principle count cell divisions or record progressing molecular identities. Finally, Bonnefont and Vanderhaeghen examine the question of temporality in progenitor competence and neuron differentiation, highlighting species-specific developmental paces and critical periods of plasticity. Molecular mechanisms controlling cell fate is discussed by Ortiz-Alvarez and Spassky while Hoye and Silver specifically address translational controls over cell fate, in particular with regard to translation in the basal processes of apical radial glia. Abnormal neuronal migration can alter cortical size, layering and folding leading to various neurological disorders. Castello and Gleeson review malformations of cortical development such as lissencephaly, epilepsy, tuberous sclerosis (TSC) and their current therapies. Biallelic mutations of new culprit genes have been recently identified, in Reelin, several genes encoding component of the extracellular matrix such as laminins and collagens or enzymes (such as TMRC3) controlling protein glycosylation. But in addition, there is mounting evidence suggesting that more focal cortical dysplasia and heterotopia, such as TSC can be attributed to the combination of germline and somatic mutations which are more difficult to mimic in animal models. The relevance of mouse models or human cortical malformations is also discussed by Francis and Capello who explore the underlying molecular mechanisms with a focus on the cytoskeleton and centrosome. Tubulinopathies comprise a large and heterogeneous range of brain disorders associated with mutations in tubulins, associated proteins and motor proteins.

The authors describe the novel mutations identified in kinesins, dyneins and a dynein adaptor or centrosomal components and how, by affecting cytoskeleton dynamics, they perturb neuronal progenitor proliferation and survival in addition to cell migration. They also emphasize the variability in the severity of these malformations across the cortex which remains to be understood. Although neuronal migration is primarily a developmental process, it persists in the rostral migratory stream of rodents, along which neuroblasts migrate from the subventricular zone (SVZ) to the olfactory bulb. Nakajima and colleagues report on the mechanisms that orient the migration of SVZ-derived neuroblasts and the influence of extrinsic factors produced by astrocytes, serotoninergic inputs, the vasculature and choroid plexus. Evidence supporting a conservation of postnatal neuroblast migration in humans, a vivid debate, is also presented. In a last part, the authors discuss recent studies showing that in various models of brain injury, SVZ-derived cells can migrate to the lesions and how the use of biomaterials could facilitate the recruitment of neuroblasts and brain repair.

Increasingly, molecules whose function was thought to be largely elucidated find themselves "repurposed" for other developmental processes, and novel axon guidance molecules emerge in addition to the "canonical" ones.

In their chapter, Guy and Kamiguchi review the unconventional role of lipids in axon guidance taking as examples Lysophosphatidic acid in the barrel cortex, sphingosine-1 phosphate and Lysophosphatidylinositol in the visual system and Lysophosphatidylglucoside in the spinal cord. Although some of their receptors have been identified, their regulation, mechanisms of action and signaling pathways in growth cones are still largely unknown. The current limitations and technical complexity of lipidomics studies is also highlighted by the authors.

Growth factors, such as VEGF, HGF, GDNF or FGF are best known for modulating cell survival and proliferation but in their chapter, Short and colleagues provide a perspective on their role in midline crossing, axon pathfinding in the peripheral nervous system and visual system. They also report on growth factor control of neuronal migration and synaptogenesis.

The pleiotropy of molecules controlling brain development is also highlighted by Faini and colleagues who present the non-canonical roles of Reelin in axonal guidance and by Limoni and Niquille, who review the diverse roles played by semaphorins and their receptors in neurogenesis, neuronal migration and oligodendrocyte development. Last, the chapter by Dorskind and Kolodkin describe recent findings that have provided novel and surprising insights into the mechanism of action of netrin-1 and DCC in axon guidance, the regulation of axon responsiveness to midline-derived cues, and the discovery of novel axon guidance molecules such as teneurins and DPRs.

Marking the transition between pre-natal and post-natal development, Gesuita and Karayannis and Genescu and Garel provide us with their view on cortical layer 1, and on how external inputs act there to sculpt circuit structure and function. The link between development and circuits is also provided by Kennedy and colleagues, who take meso-scale approaches to discuss circuit organization and Artur Kania, who provides his perspective on the development and plasticity of the somatosensory cortex. Circuit assembly and function are further discussed from a multitude of points of view: Hoersting and Schmucker discuss mounting evidence suggesting that neuropsychiatric diseases such as schizophrenia have a neurodevelopmental origin and might be related to abnormal neurite arborization. Understanding the underlying mechanisms requires the development of improved imaging techniques to access neuronal morphology and its heterogeneity at a single cell resolution. Next, they summarize the molecules and pathways that orchestrate the development of axonal and dendritic arborization and neuronal shape. In their review, Blockus and Polleux provide an integrative view on circuit assembly, with a particular focus on the hippocampus where they discuss combinatorial vs. unitary codes underlying synaptic specificity. Finally, Fossati and Charrier discuss recent work on the unexpected role of ionotropic glutamate receptors in synaptogenesis in the cerebellum, hippocampus and cortex, via a promotion of trans-synaptic interactions.

  Together, these reviews address most of the outstanding questions on circuit formation. It shows how new technologies and a multidisciplinary approach should help answering them. It also illustrates the diversity of animal models, neuronal systems and molecules that are studied in this very active research field.understanding how genetic and input-dependent mechanisms interact to control neuronal circuit assembly during cortical development. His laboratory has identified key mechanisms underlying the generation of neuronal diversity during corticogenesis, including showing that electric properties of cortical progenitors are instrumental to generate different types of neurons, identifying sequential molecular steps involved in cortical neuron generation, and revealing that cortical progenitors retain the ability to generate normally previously-born neurons. Denis Jabaudon is FENS KAVLI network of excellence alumnus and his research has earned him several awards including the Cloetta prize and Robert Bing Prize .Alain Chédotal, is research director at the French National Institute for Health and Medical Research (INSERM) and group leader at the Institut de la Vision in Paris. He was an undergrad at the Ecole Normale Supérieure de Lyon, received a PhD in Neuroscience from Pierre & Marie Curie University in Paris and conducted a postdoctoral research at The University of California at Berkeley. Dr Chédotal was elected at the French Academy of Sciences, the Academia Europaea and the European Molecular Biology Organization. He is a member of the European Dana alliance for Brain Initiatives. Dr Chédotal and his colleagues study the cellular and molecular mechanisms which control the development function and evolution of commissural neurons which interconnect the two halves of the nervous system. His team also studies the role of axon guidance molecules in retinal neo-angiogenesis, optic nerve regeneration and cornea innervation. In the past few years they have used novel molecular and imaging techniques (such as tissue clearing an 3D light sheet microscopy) to study axon guidance and embryology and they have started to build the first 3D cellular atlas of the developing human embryo.

	Alain Chédotal
	Institut de la Vision, INSERM, Sorbonne Université, Paris
	Email: alain.chedotal@inserm.fr

and also practices as a clinical neurologist at Geneva University Hospital. His work is funded by the Swiss National Science Foundation and the European Research Council (Advanced grant). Pr. Jabaudon's research focuses on