V Resseguier
email: valentin.resseguier@scalian.com

L Li

• G Jouan

P Dérian

E Mémin

B Chapron

G Jouan

• E Mémin

New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics

Keywords: Uncertainty quantification, ensemble forecast, fluid dynamics, data assimilation, model under location uncertainty

Numerical simulations of industrial and geophysical fluid flows cannot usually solve the exact Navier-Stokes equations. Accordingly, they encompass strong local errors. For some applications -like coupling models and measurements -these errors need to be accurately quantified, and ensemble forecast is a way to achieve this goal. This paper reviews the different approaches that have been proposed in this direction. A particular attention is given to the models under location uncertainty and stochastic advection by Lie transport. Besides, this paper introduces a new energy-budget-based stochastic subgrid scheme and a new way of parameterizing models under location uncertainty. Finally, new ensemble forecast simulations are presented. The skills of that new stochastic parameterization are compared to that of the dynamics under location uncertainty and of randomized-initial-condition methods.

Introduction

Understanding, modeling and tracking high Reynolds flows remain main challenges in current researches. Indeed, beyond economical applications related to weather forecasting and industrial flows analysis, accurate climate projections have become a societal need. The complexity of such systems is mainly due to the non-linear and non-local nature of the evolution laws. These features make large-scale flow structures interact with smaller ones. As such, the large-scale flow components cannot be simulated without the small-scale components. However, the computational expense of solving all the hydrodynamical scales is still today beyond of reach even for turbulent flows of moderated complexity. As an example, today's most accurate oceanic currents numerical simulations use a horizontal mesh resolution of about 1 km [START_REF] Klein | Upper ocean turbulence from high-resolution 3D simulations[END_REF][START_REF] Gula | Gulf stream dynamics along the southeastern us seaboard[END_REF], whereas solving the "real" equations of fluid dynamics, say the Navier-Stokes equations, would require a grid cell of about 1 cm.

The effects of the unresolved so-called turbulent small-scale fluctuations have to be modeled. Turbulent dissipations, advection corrections and backscatterings need to be introduced to mimic the action of the small-scale processes on the large-scale components. They respectively reduce, move and increase energy of large-scale tracers. But, even using the best subgrid model, the true dynamics cannot be fully recovered by a large-scale model. Strong errors remain as only the unresolved dissipation is in general taken into account. The unresolved subgrid components of tracer or velocity remains by definition unknown. In other words, they are uncertain. In this review, uncertainty and stochasticity does not mean noise but rather something that is only known through a probability distribution. Uncertainty can be attached to coherent structures in time and space. To better express this idea, let us introduce a comparison. When looking at a tree, you cannot see all the branches because they are either too small or hidden by leaves. The particular shapes of these branches are uncertain and could be modeled by random processes. However, you know that they respect some features. For instance, each of them is linked to the trunk through one and only one path. These branches are uncertain coherent structures in the same way as unresolved small-scale vortices are. Both must respect appropriate physical laws. Since resolved and subgrid fluid dynamics are coupled, the large-scale resolved dynamics must explicitly take into account this stochastic nature. Understanding, modeling and simulating this randomness or errors is the subject of Uncertainty Quantification (UQ). For sake of concision, this review will not deal with the random models specifically introduced to study extreme events since there is a full literature on this subject (e.g. [START_REF] Franzke | Extremes in dynamic-stochastic systems[END_REF].

The paper is organized as follows. We first motivate the need for UQ in CFD through the data assimilation application. Then, section 3 describes some widely-used CFD UQ tools which are not a priori related to the numerical scale truncation. Section 4 reviews the historical stochastic subgrid parameterizations, which were not originally introduced for UQ purposes. Next, section 5 enumerates the more recent schemes designed for UQ purposes in particular the physically-based ones. In this section, we also propose a new stochastic subgrid parameterization adapted to the CFD dissipation scheme through an imposed energy budget. Section 6 focuses on a specific physically-based family of stochastic subgrid tensors: the models under Location Uncertainty (LU) and the Stochastic Advection by Lie Transport (SALT). After this, section 7 details the metrics used to assess the ensemble forecasts' qualities. Finally, new numerical results are presented for short-term and long-term ensemble forecasts. Several stochastic subgrid parameterizations are compared on a simplified model of geophysical fluid flow dynamics.

Data assimilation, a main motivation

To introduce this review topic, let us first give a brief description about data assimilation and filtering method. Even beyond fluid dynamics applications, coupling numerical model simulations and measurements is of great interest. This is called data assimilation. Some of these methods, derived from variational principle and optimal control theory, are deterministic (Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF]. A functional criterion is optimized in terms of control variables (such as initial conditions) to drive the model as closely as possible to the observations. For instance, in 4D-Var algorithms solutions are functions of time and space but also of the initial conditions. Variational optimization performed using the adjoint tangent dynamical model allows to infer an initial condition with a trajectory that matches at best the observations (with respect to a given distance). This review will instead be placed in the framework of probabilistic data assimilation methods, such as smoothing and above all filtering [START_REF] Doucet | Stochastic theory of compressible turbulent fluid transport[END_REF][START_REF] Doucet | A tutorial on particle filtering and smoothing: Fifteen years later[END_REF][START_REF] Candy | Bayesian signal processing: Classical, modern and particle filtering methods[END_REF][START_REF] Papadakis | Data assimilation with the weighted ensemble Kalman filter[END_REF]. These methods allow combining a random dynamical model together with noisy and partial observations of the system of interest to drive an ensemble forecast -a set of realizations -along time. Hybrid approaches also exist, such as 4DEnVar [START_REF] Buehner | Ensemble-derived stationary and flow-dependent background-error covariances: Evaluation in a quasi-operational nwp setting[END_REF][START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part ii: Observing system simulation experiments with advanced research wrf (arw)[END_REF][START_REF] Yang | Study of variational ensemble methods for image assimilation[END_REF][START_REF] Bocquet | An iterative ensemble kalman smoother[END_REF][START_REF] Yang | Study of variational ensemble methods for image assimilation[END_REF]Mémin, 2017, 2018;[START_REF] Sakov | An iterative ensemble kalman filter in the presence of additive model error[END_REF].

In probabilistic data assimilation algorithms, the randomness of the model is meant to capture the errors of the dynamical system whereas randomness of observations represents the measurement errors. If the model noise is prominent, the estimated filtered variable relies principally on the observations. Conversely, when the observation noise is dominant, the filtered variable trajectory is mainly driven by the model. For this reason, an accurate design of the model errors is crucial in weather and climate communities [START_REF] Allen | Towards objective probabilistic climate forecasting[END_REF]Penland, 2003a,b;[START_REF] Berner | Stochastic parameterization: towards a new view of weather and climate models[END_REF]. Furthermore, due to the huge state-space dimension (∝ 10 7-9 ), the ensemble size is usually very small in comparison (∝ 10 1-2 ). Thus, the randomness of the dynamical model has to be very efficient. Ensemble members, also called particles, have to quickly spread in the phase space. At the same time, each particle should remain a "physically plausible realization" to focus on meaningful regions of the state space.

We will now describe several ways to represent the randomness of the fluid dynamical model. That randomness can have several sources: e.g. wrong physics, unknown parameters or initial conditions. We first detail the approaches which are not a priori related to the unresolved components of the dynamics. Then, we will review random parameterizations which model energy backscatterings and/or errors induced by the subgrid dynamics.

3 UQ not directly related to the coarse-resolution-induced errors 3.1 UQ from parameters, boundary conditions and forcings Some authors inject randomness through the parameters. Indeed, parameters like viscosity, initial and boundaries conditions are often assumed random (e.g. Le [START_REF] Maitre | A stochastic projection method for fluid flow. II. random process[END_REF][START_REF] Sapsis | Dynamical criteria for the evolution of the stochastic dimensionality in flows with uncertainty[END_REF]. The chaotic nature of fluid dynamics increases quickly the eventual errors associated to these parameters. Other authors study the uncertainty arising from forcings. In particular, CO 2 concentration conditions are difficult to specify in climate sciences. [START_REF] Lucarini | Mathematical and physical ideas for climate science[END_REF] approach this problematic with Ruelle response theory.

Random initial conditions

In fluid dynamics, random initial conditions have been widely used for both UQ and predictability studies (e.g. [START_REF] Métais | Statistical predictability of decaying turbulence[END_REF]. For a long time, operational weather forecast centers had relied on random perturbations of initial conditions to spread the ensemble forecasts [START_REF] Ehrendorfer | Predicting the uncertainty of numerical weather forecasts: A review[END_REF]Buizza, 2016b). Different types of perturbations were proposed to that end. For instance, the European center (ECMWF) relied on the so-called singular vectors method (SV) while the American center (NCEP) used the bred vectors (BV). The common idea is to perturb the initial condition along the few directions that will lead to a maximum ensemble variance in the next future. The SV are obtained by maximizing the ensemble variance after a finite-time linearized simulation. The linearization enables to solve the optimization problem with the adjoint equations and a singular value decomposition. The amplitude of the perturbation is then tuned by inflation (see later in this section for a description of the inflation method). In contrast, the bred vectors computation does not need an adjoint code nor a linearization. Two simulations -perturbed and not perturbed -are launched from a previous time. At the current time of interest -the time of the forecast initial condition -the difference between the two simulations is rescaled to give the bred vectors.

In the geophysical data assimilation communities, the initial condition random perturbation is now known to be underdispersive, i.e. it underestimates the errors related to the coarsening of the state variable dynamics [START_REF] Berner | Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations[END_REF][START_REF] Mitchell | Data assimilation in slow-fast systems using homogenized climate models[END_REF][START_REF] Gottwald | The role of additive and multiplicative noise in filtering complex dynamical systems[END_REF][START_REF] Franzke | Stochastic climate theory and modeling[END_REF]. As such the model is overconfident. When an observation -often far from the ensemble -is assimilated, only very few particles -and in the worst case only one -are considered relevant. The others are discarded. This degeneracy is referred to as filter divergence. The weakness of the method can be explained by at least two facts. At the initial time, the random perturbations live in a huge state space. Computational limitation leads to the setup of only a small-size ensemble, thus spanning only a small part of the state-space [START_REF] Mitchell | Data assimilation in slow-fast systems using homogenized climate models[END_REF][START_REF] Gottwald | The role of additive and multiplicative noise in filtering complex dynamical systems[END_REF]. As a consequence, without any adhoc compensation, the ensemble variance is underestimated. Moreover, the initial condition random perturbations are injected at small scales and are hence quickly diffused by the subgrid tensor. From a dynamical system point of view, subgrid tensor makes small-scale components of the solution more stable. Without the fully resolved non-linear mechanisms yielding an energy redistribution, the particles tend to asymptotically align with the most unstable directions [START_REF] Trevisan | Assimilation of standard and targeted observations within the unstable subspace of the observation-analysis-forecast cycle system[END_REF]Trevisan and Palatella, 2011;[START_REF] Mitchell | Data assimilation in slow-fast systems using homogenized climate models[END_REF][START_REF] Sapsis | Attractor local dimensionality, nonlinear energy transfers and finite-time instabilities in unstable dynamical systems with applications to two-dimensional fluid flows[END_REF][START_REF] Gottwald | The role of additive and multiplicative noise in filtering complex dynamical systems[END_REF]. Thus, as time evolves the ensemble spans a smaller and smaller space.

To mitigate the ensemble variance underestimation, a famous compensation method exists. It is called covariance inflation [START_REF] Anderson | A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts[END_REF]. The ensemble covariance is increased (in an additive or multiplicative way) by a carefully-tuned parameter. This scalar can be identified through statistics estimation in the observation space [START_REF] Tandeo | Joint estimation of model and observation error covariance matrices in data assimilation: a review[END_REF]. Indeed, in this space, the innovations are the differences between the predicted and the actual observations. By simple algebra, the variance estimation of that innovation -before the analysis step -can be related to the missing inflation factor, the badly-scaled state-vector covariance, the measurement-error covariance, and the matrix linking state vector and observations. Some operational weather forecasts centers rely on inflation [START_REF] Raynaud | Model error representation in meteo france ensemble nwp system[END_REF]. In any cases, the ensemble covariance is often erroneous and such a compensation may lead to nonphysical behaviors, as exemplified in figure 3.2. Indeed, an ensemble spread underestimation -says a factor 2 -due to for instance an overconfidence in the position of an eddy would lead to a (variance) inflation factor of 4. Thus, in one realization of the ensemble, if the eddy does not overlap with mean eddy, it will stay in the same place but will become 2 times stronger. Therefore, the ensemble mean square error will increase without pushing any realization closer to the truth. A popular inflation variant was introduced by [START_REF] Desroziers | Diagnosis of observation, background and analysis-error statistics in observation space[END_REF]. From the a priori innovation covariance and the crosscovariance between a priori and a posteriori innovations, it is possible to correct both the badly-scaled state-vector covariance and the measurement-error covariance. Nevertheless, similar drawbacks may be expected.

Stochastic backscattering

Other works preferably address the modeling of errors related to wrong dynamics. Indeed, as long as all the scales are not resolved, subgrid dynamics are modeled rather than resolved. This introduces errors which grow in time due to chaotic behaviors. A natural way to address this UQ is to continuously introduce noise in the dynamics. Evolution laws no longer rely on bulk pa-Fig. 1 Inconsistency introduced by covariance inflation for a 2-member ensemble of temperature fields (in Celsius degrees). Here, the optimal inflation factor is about 1.8. The inflation increases the variance without reducing the bias (even though the bias will be modified in the following data assimilation analysis step). Thus, it increases the point-wise Root Mean Square Error (RMSE, see section 7.2 for a definition). Here, inflation makes ensemble members further away from the reference. Moreover, the inflation artificially creates water at 21.1 • and at 14 • from a water at temperature T ∈ [15.5 rameterizations i.e. deterministic models. They become Langevin equations expressed through Stochastic Differential Equations (SDE) or Stochastic Partial Differential Equations (SPDE). Note that the Langevin equations can be non-linear with complicated random forcing and memory terms.

Langevin equations are not new in fluid dynamics. Stochastic subgrid parameterizations have been developed before the use of ensemble forecasts in CFD. Nevertheless, most of these historical stochastic schemes have not been developed for UQ purposes -even though they can often be used to that end.

These first stochastic subgrid parameterizations generally try to mimic the intermittent back-scattering of energy from small scales toward larger scales. Several recent UQ methods are also based or inspired from the stochastic backscattering literature. We briefly describe this literature below.

Stochastic Lagrangian methods

First, let us mention the stochastic Lagrangian models introduced by S. B. Pope and coauthors (see [START_REF] Pope | Lagrangian PDF methods for turbulent flows[END_REF], for a review). The Lagrangian fluid particles represent independent realizations. Each particle moves with a local velocity which is itself the solution of a randomized Navier-Stokes model in a Lagrangian form. In addition to the mean pressure forcing, the authors assume a relaxation mechanism toward the mean velocity and an additive space-time white noise. Its amplitude is proportional to the dissipation D . In order to compensate the induced mean kinetic energy increase, the relaxation term includes a term proportional to D /TKE where TKE denotes the turbulent kinetic energy. Depending on the complexity of the relaxation matrix, the full model is referred to as simplified Langevin model (SLM) or generalized Langevin model (GLM) [START_REF] Pope | A lagrangian two-time probability density function equation for inhomogeneous turbulent flows[END_REF]. In the second case, the relaxation matrix also depends on the mean shear and mean Reynolds stress tensor. The mean terms can be computed by averaging over particles since they are independent. As all coefficients of the SLM and GLM are deterministic, the solution is a Gaussian process. The Langevin models are simulated through Monte Carlo Markov Chain (MCMC). The refined Langevin model [START_REF] Pope | The velocity-dissipation probability density function model for turbulent flows[END_REF]) is a variant which replaces the mean dissipation by a stochastic one in the SLM and GLM. The dissipation is assumed to follow another stochastic equation which makes the dissipation probability law close to a log-normal distribution (i.e. the logarithm of the dissipation is approximately Gaussian). [START_REF] Sawford | Reynolds number effects in lagrangian stochastic models of turbulent dispersion[END_REF] also proposes a Lagrangian model where it is the Lagrangian acceleration which follows a (Gaussian) SLM-like model. Note that considering an acceleration evolution equation is formally equivalent to the inclusion of memory terms in a velocity evolution equation. [START_REF] Berloff | Material transport in oceanic gyres. part II: Hierarchy of stochastic models[END_REF]; [START_REF] Veneziani | Oceanic turbulence and stochastic models from subsurface Lagrangian data for the Northwest Atlantic Ocean[END_REF] also consider Stochastic Lagrangian models. The slightly non-linear resulting evolution models are defined on empirical grounds and Gaussian assumptions.

Eulerian Gaussian backscattering

The most famous Eulerian backscattering method is perhaps the Eddy-Damped Quasi Normal Markovian (EDQNM) model introduced by [START_REF] Orszag | Analytical theories of turbulence[END_REF] and [START_REF] Leith | Atmospheric predictability and two-dimensional turbulence[END_REF]. It closes the large-scale Navier-Stokes equations in the Fourier space by neglecting some phase correlations in the non-linear terms while keeping the energy constant. [START_REF] Chasnov | Simulation of the Kolmogorov inertial subrange using an improved subgrid model[END_REF] uses this framework to set up a forced-dissipative Navier-Stokes model where the Eddy Viscosity is scaledependant and the forcing is Gaussian, homogeneous and isotropic in space and uncorrelated in time. Replacing non-linear interactions by a damping term and a Gaussian forcing is now common practice. This is in particular used to setup evolution laws of subgrid scales where accuracy is of lower concern. The solutions are in this case Gaussian processes. Structural Stability Theory (S3T) is one example in that spirit [START_REF] Farrell | Statistical state dynamics: a new perspective on turbulence in shear flow[END_REF]. The Quasi-Linear (QL) approximation separates the non-linear deterministic dynamics of the mean field and the linearized randomly-forced dynamics of fluctuations. Then, stability analyses are applied to the augmented state-vector of mean and covariance in order to reveal and characterize various turbulent phenomenon. Stochastic superparameterization (SSP) proposes a similar model [START_REF] Grooms | Stochastic superparameterization in quasigeostrophic turbulence[END_REF]. The point approximation separates the large-scale and the small-scale dynamics. The small-scale evolution law is linearized and corrected by the introduction of noise and damping terms. To ensure energy conservation in the stationary regime, noise variances and damping coefficients are related. Then, the second order moments of the solution are known analytically and can feed the subgrid tensor expression of the mean large-scale evolution law. Without involving any theoretical closure, [START_REF] Berloff | Random-forcing model of the mesoscale oceanic eddies[END_REF] considers a Gaussian forcing as well. Yet its noise is inhomogeneous in space and correlated in time.

Well specified inhomogeneity brings phase information making the model more accurate and the forcing more efficient.

Eulerian non-Gaussian backscattering

Phase information can also be encoded by multiplicative noises. Besides, multiplicative noises are the most common non-Gaussian forcing in CFD. [START_REF] Leith | Stochastic backscatter in a subgrid-scale model: Plane shear mixing layer[END_REF] multiplies a white Gaussian noise by a function of the resolved local strain rate. [START_REF] Schumann | Stochastic backscatter of turbulence energy and scalar variance by random subgrid-scale fluxes[END_REF] uses a quadratic function of a homogeneous Gaussian noise. [START_REF] Brankart | Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modeling[END_REF] adds at each time step a multiplicative noise to the active tracers (salinity and temperature) before computing the corresponding density. After this, the density is averaged over realizations. The nonlinearity of the state-equation makes this transitional variability non-negligible and improves the simulation results.

In contrast, [START_REF] Mana | Toward a stochastic parameterization of ocean mesoscale eddies[END_REF] use a non-Gaussian noise process. It is uncorrelated in time and space with a variance that depends on the resolved potential vorticity (PV) gradient.

Except [START_REF] Schumann | Stochastic backscatter of turbulence energy and scalar variance by random subgrid-scale fluxes[END_REF], all these methods defined on empirical grounds have a common characteristic: the factor of the multiplicative noise -or the noise variance -is a function of the gradient of the transported quantity. This comment is also valid for the stochastic Lagrangian models of S. B. Pope and coauthors. This suggests a link between stochastic backscattering and turbulent dissipation, but the justification of this link is generally either unclear or arbitrary.

UQ for errors induced by a coarsening process

Gaussian additive noise

The simplest random dynamics are defined from linear Langevin equations with additive Gaussian noise. This has already been discussed above for stochastic backscattering methods. We may add the linear inverse models [START_REF] Penland | A balance condition for stochastic numerical models with application to the El Nino-southern oscillation[END_REF][START_REF] Penland | The optimal growth of tropical sea surface temperature anomalies[END_REF]. [START_REF] Keating | New methods for estimating ocean eddy heat transport using satellite altimetry[END_REF] also use a linear evolution model with Gaussian noise for a filtering purpose. However, the parameters of the models are themselves Ornstein-Uhlenbeck (OU) processes. An OU process is a gaussian process with an exponential covariance in the stationary regime. Such a process corresponds to an Auto-Regressive (AR) process in discrete time. In [START_REF] Keating | New methods for estimating ocean eddy heat transport using satellite altimetry[END_REF], OU parameters make the evolution model solution not Gaussian and more flexible. This method -well known in the filtering community -is referred to as Stochastic Parameterized Extended Kalman Filter (SPEKF) [START_REF] Gottwald | The role of additive and multiplicative noise in filtering complex dynamical systems[END_REF].

In geophysical fluid dynamics data assimilation, dynamics are often nonlinear. Nevertheless, Gaussian additive random forcings are still widely used. [START_REF] Tandeo | Joint estimation of model and observation error covariance matrices in data assimilation: a review[END_REF] review the additive-Gaussian-noise-covariance estimation procedures. Moment-based, maximum likelihood, and Bayesian estimations are detailed. Several authors also define correlation matrices as solutions of diffusion equations parameterized from observed correlation lengths [START_REF] Pannekoucke | Estimation of the local diffusion tensor and normalization for heterogeneous correlation modelling using a diffusion equation[END_REF]. Others learn a spatial mapping which transforms the heterogeneous and anisotropic spatial covariance into a homogeneous and isotropic one using wavelets (Michel, 2013a,b) or differential geometry tools [START_REF] Pannekoucke | Modelling of local length-scale dynamics and isotropizing deformations[END_REF]. Then, a classical covariance model can be fitted.

Empirical non-Gaussian noises

Among other empirical stochastic models, the Stochastic Perturbed Physics Tendency scheme (SPPT) introduces a correlated multiplicative noise [START_REF] Buizza | Stochastic representation of model uncertainties in the ECMWF ensemble prediction system[END_REF]. All parameterized tendencies are multiplied by the same uniform variable in [0.5, 1.5]. The same random variable is used in a fixed spacetime window. This window sets the noise spatio-temporal correlations. [START_REF] Shutts | A kinetic energy backscatter algorithm for use in ensemble prediction systems[END_REF]; [START_REF] Berner | A spectral stochastic kinetic energy backscatter scheme and its impact on flow-dependent predictability in the ECMWF ensemble prediction system[END_REF][START_REF] Berner | Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations[END_REF] have proposed another method constructed from a spatially homogeneous and isotropic Gaussian noise. Then, they multiply it by the dissipation rate. This method is called Stochastic Kinetic Energy Backscatter (SKEBS). As in [START_REF] Schumann | Stochastic backscatter of turbulence energy and scalar variance by random subgrid-scale fluxes[END_REF] and [START_REF] Brankart | Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modeling[END_REF], the noise is defined from an Ornstein-Uhlenbeck (OU) process with a very small correlation time. As for many stochastic subgrid parameterization, the factor of the multiplicative noise is again a function of the gradient of the transported quantity.

We may also cite a recent variant of SKEBS referred to as Stochastic Convective Backscatter (SCB) scheme [START_REF] Shutts | A stochastic convective backscatter scheme for use in ensemble prediction systems[END_REF]. This scheme is focused on atmospheric flow applications. A Gaussian noise white in space and time is weighted by the vertical variation of the parameterized convective updraught and downdraught mas fluxes. Then, a spatio-temporal smoothing imposes the spatial and temporal correlations.

SPPT and SKEBS methods have been successfully applied in operational weather and climate forecast centers [START_REF] Franzke | Stochastic climate theory and modeling[END_REF]Buizza, 2016b;[START_REF] Pegion | Stochastic parametrization development in the noaa/ncep global forecast system[END_REF][START_REF] Separovic | Model error representation in the canadian ensemble prediction systems[END_REF][START_REF] Mccabe | Representing model error in the met office convection permitting ensemble predicion system[END_REF]Buizza, 2016a). Nevertheless, many drawbacks of these methods have been reported. Above all, conservation laws (e.g. energy and mass conservations) are violated [START_REF] Reynolds | Reports from working group 3 : What are the pros/cons of existing model uncertainty schemes and how should these be measured?[END_REF][START_REF] Leutbechner | Stochastic representations of model uncertainties in the ifs[END_REF]. Precipitation biases have also been observed (Buizza, 2016a). Moreover, the random forcing scales differ from the scales of the error sources that are meant to be quantified [START_REF] Reynolds | Reports from working group 3 : What are the pros/cons of existing model uncertainty schemes and how should these be measured?[END_REF].

Energy-budget-based noises

Energy conservation and redistribution being fundamental aspects of physics, several authors have developed ad hoc schemes to deal with conservation laws.

Existing schemes

Sapsis and Majda (2013c) introduce the Modified Quasilinear Gaussian (MQG) model for -but not restricted to -dimensionally reduced systems (Sapsis and Majda, 2013b,a). This model approximates the third-order moment in the covariance evolution law in order to redistribute the right amounts of energy between modes. Based on stationary regime information, dampings and noises are specified in that way.

Earlier in this paper, we have mentioned the Stochastic superparameterization (SSP) [START_REF] Grooms | Stochastic superparameterization in quasigeostrophic turbulence[END_REF], the EDQNM approximation and other related methods. In these approaches, the random forcing is also specified from an imposed mean energy balance. Nevertheless, these methods are not a priori designed for an UQ purpose.

Similar energy mean budgets have recently been discussed by several authors. [START_REF] Majda | Statistical energy conservation principle for inhomogeneous turbulent dynamical systems[END_REF] refers to this energy mean as the statistical energy. The author derives the evolution law of this energy by adding the evolution equations of the mean energy and the evolution equation of the integrated variance. However, [START_REF] Majda | Statistical energy conservation principle for inhomogeneous turbulent dynamical systems[END_REF] does not specify the random forcing. This is why the latter does not a priori balance the turbulent diffusion. [START_REF] Farrell | Statistical state dynamics: a new perspective on turbulence in shear flow[END_REF] also study the energy mean of stochastic fluid dynamics systems especially under quasi-linear approximations and with an additive Gaussian forcing. Gugole and Franzke (219) enforce the energy conservation by projecting their random forcing into the space of constant energy.

As cited previously, SKEBS parameterization [START_REF] Shutts | A kinetic energy backscatter algorithm for use in ensemble prediction systems[END_REF][START_REF] Berner | A spectral stochastic kinetic energy backscatter scheme and its impact on flow-dependent predictability in the ECMWF ensemble prediction system[END_REF][START_REF] Berner | Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations[END_REF] backscatters a given portion of the global dissipated energy using a flow-dependent correlated noise. [START_REF] Jansen | Parameterizing subgrid-scale eddy effects using energetically consistent backscatter[END_REF] work with a similar idea. Yet, they consider a hyper-viscous diffusion and a noise white in space and time. [START_REF] Dwivedi | Energetically consistent scale adaptive stochastic and deterministic energy backscatter schemes for an atmospheric model[END_REF] also consider hyper-viscosity and neglect spatial correlations but they keep time correlations.

A new random forcing derived from large-scale subgrid dissipation models

We here propose a new energy-budget-based stochastic subgrid model.

In practice, large-scale fluid dynamics models involve dissipation operators called subgrid models. They stabilize the numerical simulations and mimic the action of the unresolved small scales by draining the energy at high wavenumbers. Simple subgrid models often take the form of classical Laplacian operator or of higher-order hyperviscosity operators (typically, some power of a Laplacian). Let us consider the dynamic of a fluid property q. It is assumed to be transported up to a dissipation operator denoted L and a stochastic forcing η -which is a centred process uncorrelated in time:

Dq Dt = L[q] + η, (1)
where D/Dt = ∂ t + v • ∇ stands for the material derivative operator. Through Itô calculus, a conservation of the energy mean -E Ω 1 2 q 2 -would imply:

0 = E Ω d dt ( 1 2 q 2 ) = E Ω (q d dt q + 1 2 d dt < q, q >) = E Ω qL[q] + 1 2 (H[q]) 2 dt , (2)
where

(H[q](x, t)) 2 = E(η 2 (x, t)|q(•, t)), (3)
is the stochastic forcing variance conditioned on field q at time t and < q, q > denotes the quadratic variation term. In order to maintain a desired amount of energy dissipation, we introduce a scaling factor ζ > 0 in the above balance. Specifically, we want to impose:

Ω q -2ζ dt L[q] = Ω (H[q]) 2 . ( 4
)
Now, for a given dissipation operator L, we will construct a noise η which meets the desired balance (4). If L is a negative auto-adjoint operator -which is generally the case -we can define an operator H which is such that

-2ζ dt L = HH , (5)
with H denoting the adjoint of operator H. More precisely, for the Laplacian, bi-Laplacian or higher-order hyperviscosity operators, this operator H reads,

H[q] =              ζα∇q if L[q] = ∇ • (αα T ∇q), ζ∇ • (α∇q) if L[q] = -∇ • α∇ ∇ • (α∇q) , ζα∇∆ p q if L[q] = α 2 ∆ 2p+1 q, ζα∆ p q if L[q] = -α 2 ∆ 2p q, (6)
where ζ = 2ζ dt , p is a positive integer, α a constant and α a -possibly spatially-varying -matrix.

Then, for specific deterministic subgrid tensor L, if the stochastic forcing (conditional) variance integral, Ω (H[q]) 2 , can be set to the integral Ω H[q] 2 , it is possible then to control the energy dissipation and meet the desired balance (4). Indeed, in that case, we have:

Ω (H[q]) 2 = Ω H[q] 2 = Ω qH H[q] = Ω q -2ζ dt L[q] . (7)
Nevertheless, to simulate the stochastic forcing, knowing its global variance is not enough. We also need to choose its spatial structure (e.g. the local variation of the variance, its spatial correlations). In this purpose, we express the noise on a convenient orthonormal basis of L 2 (Ω) (e.g. a wavelet basis)

denoted {e k } k≥0 : η(x, t) = k≥0 λ k (t)ξ k (t)e k (x), (8)
where the ξ k are independent Gaussian white noises of unit variance and:

λ k = Ω H[q] e k . (9)
Then, by definition of H and by the Parseval theorem, the random forcing (conditional) variance integral is given by:

Ω (H[q]) 2 = Ω E(η 2 (x, t)|q(•, t)) = k≥0 |λ k | 2 = Ω H[q] 2 . ( 10
)
Therefore, equation ( 7) is valid and the stochastic forcing η defined by ( 8) and (9) meets our target energy balance (4). Therefore, given a numerical dissipation L together with a chosen orthonormal basis, it is always possible to define a noise (8) that respects the assumed energy balance (4). Note however that this balance constitutes only a global balance. Indeed, locally in space, the random forcing variance is: (11) and moreover,

(H[q]) 2 = k≥0 |λ k | 2 |e k | 2 = k≥0 λ k e k 2 = H[q] 2 ,
H[q] 2 = qH * H[q] = q -2ζ dt L[q] . (12)
The choice of a given basis {e k } k≥0 in this model is of crucial importance, as it influences in particular the variance heterogeneity and the correlation lengths of the random forcing. In order to be close to a local energy balance, basis functions with small supports can be selected. This enables us to restrain the inequality (11) and make it closer to an equality. For instance, Fourier modes have a large support and would lead to a homogeneous random forcing (since |e k | = cst). In this case, the inequality (11) would be (H

[q]) 2 = 1 µ(Ω) Ω H[q] 2 = H[q] 2 (x, t)
, where µ(Ω) is the measure of the spatial domain Ω. In contrast a basis of regularized Dirac functions e k (x) = δ(x -x k ) -defined on the grid points x k -have near zero-measure supports. The inequality (11) would hence become an equality. Nevertheless, such an infinitesimally small basis function support would induce an almost zero random forcing correlation length. We consider this behavior as nonphysical. According to that analysis, a wavelet basis seems a promising trade-off and will be used in the numerical section 8.

Physically-based noises

We now come back to the literature review. In fluid dynamics, due to nonlinearities, the targeted probability density of the solution is highly non-Gaussian and relevant stochastic dynamical models are difficult to derive. In this perspective, an attractive path would be to infer randomness from physics [START_REF] Berner | Stochastic parameterization: towards a new view of weather and climate models[END_REF][START_REF] Craig | Reports from working group 1 : What are the sources of model error and how can we improve the physical basis of model uncertainty representation?[END_REF]). Yet as Navier-Stokes equations are deterministic, this path is not straightforward.

Time-scale separation

Averaging and homogenization

To derive large-scale fluid dynamical model, the time-scale separation assumption is convenient. In the seventies, [START_REF] Hasselmann | Stochastic climate models. part I[END_REF] already relied on it for geophysical fluid dynamics. In his seminal work, the large-scale dynamics were encoded by both mean terms and noise terms. However, eventually only simple multidimensional OU processes were considered.

The time-scale separation assumption is also the foundation of the more rigorous averaging and homogenization theories [START_REF] Kurtz | A limit theorem for perturbed operator semigroups with applications to random evolutions[END_REF][START_REF] Papanicolaou | Asymptotic theory of mixing stochastic ordinary differential equations[END_REF][START_REF] Givon | Extracting macroscopic dynamics: model problems and algorithms[END_REF][START_REF] Mitchell | Data assimilation in slow-fast systems using homogenized climate models[END_REF][START_REF] Gottwald | Homogenization for deterministic maps and multiplicative noise[END_REF][START_REF] Gottwald | The role of additive and multiplicative noise in filtering complex dynamical systems[END_REF][START_REF] Franzke | Stochastic climate theory and modeling[END_REF][START_REF] Gottwald | Stochastic climate theory[END_REF]. As the time-scale separation goes to infinity, the large-scale dynamics will converge according to averaging or to homogenization depending on the structure of the global model. The global dynamics as well as the limit large-scale dynamics can be differential equations or SDEs. In the large-scale equation, terms which are only functions of the small-scale variable often tend to converge to a white-noise-in-time term in the Stratonovich sense with a covariance of the Green-Kubo type. Nevertheless, it is not always true for nonlinear dynamics. Sometimes, the noise has to be understood in the sense of Ito or Marcus stochastic integral. In the last case, the noise is a Levy process [START_REF] Gottwald | Homogenization for deterministic maps and multiplicative noise[END_REF][START_REF] Gottwald | Stochastic climate theory[END_REF].

A successful application of the homogenization theory in geophysics is the MTV scheme [START_REF] Majda | Models for stochastic climate prediction[END_REF][START_REF] Majda | A mathematical framework for stochastic climate models[END_REF][START_REF] Franzke | Low-order stochastic mode reduction for a realistic barotropic model climate[END_REF][START_REF] Majda | An applied mathematics perspective on stochastic modelling for climate[END_REF]. MTV refers to the names of the three main authors: Majda, Timofeyev and Vanden-Eijnden. In practice, the non-linearity of the small-scale equation is empirically replaced by a noise term and a damping term before the homogenization procedure. The homogenized dynamics obtained are cubic with correlated additive and multiplicative (CAM) noises. Even without dealing with Levy processes, this structure is able to produce intermittency and extreme events especially because of the CAM noise. This specific form has also been used to infer data-driven models. [START_REF] Peavoy | Systematic physics constrained parameter estimation of stochastic differential equations[END_REF] proposed an example of such a model which uses energy constrained Bayesian estimators and artificial additional observations through Brownian bridge.

Another method called invariant manifold theory also invokes a time-scale separation. Yet it relies directly on the SDE solution rather than on its probability density function [START_REF] Givon | Extracting macroscopic dynamics: model problems and algorithms[END_REF][START_REF] Gottwald | The role of additive and multiplicative noise in filtering complex dynamical systems[END_REF]. This scheme provided good UQ skills especially because of the multiplicative noise appearing in the limit dynamics. [START_REF] Chekroun | The emergence of fast oscillations in a reduced primitive equation model and its implications for closure theories[END_REF] describes a similar problem for geophysical fluid dynamics on the so-called parameterized manifolds with small Rossby number.

To conclude on these methods, for complicated non-linear dynamical systems, it is still not clear how to perform homogenization and when this is possible. Moreover, the theory does not make the noise covariance explicit enough and it has to be estimated on data. During this step, some Gaussian approximations are usually done when estimating the coefficients of the model. Some homogenization methods like the MTV algorithm may suffer from energy-conservation issues. Nevertheless, workarounds exist [START_REF] Frank | Stochastic homogenization for an energy conserving multiscale toy model of the atmosphere[END_REF][START_REF] Jain | Stochastic mode-reduction in models with conservative fast sub-systems[END_REF]. In addition, the homogenization methods have shown successful results in the context of reduced order models and suggest that geophysical stochastic fluid dynamic models should involve CAM noises.

Edgeworth expansion

In order to alleviate the time scale gap assumption of averaging and homoge-nization procedures, Edgeworth expansions goes to higher order in time scale ratio [START_REF] Wouters | Stochastic model reduction for slow-fast systems with moderate time-scale separation[END_REF].

Velocity time scale separation and skew-symmetric noises

The Kraichnan model is an idealization of passive tracer turbulence [START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF][START_REF] Kraichnan | Anomalous scaling of a randomly advected passive scalar[END_REF][START_REF] Gawȩdzki | Anomalous scaling of the passive scalar[END_REF]. A tracer is forced and advected by two independent spatially-homogeneous time-uncorrelated Gaussian processes. In contrast to the previous section, the time scale gap between the unresolved and the resolved dynamics is directly assumed, without going to a limit. The forcing of the Kraichnan model lives at large spatial scales whereas the energy of the random advecting velocity is distributed over the spatial scales in a self-similar way. Since this velocity is also assumed divergence-free, the advecting term is here a skew-symmetric multiplicative noise of the tracer dynamics. It is a well-studied and well-known model in statistical physics [START_REF] Gawȩdzki | Anomalous scaling of the passive scalar[END_REF][START_REF] Klyatskin | Short-time correlation approximations for diffusing tracers in random velocity fields: A functional approach[END_REF][START_REF] Majda | Models for stochastic climate prediction[END_REF][START_REF] Falkovich | Particles and fields in fluid turbulence[END_REF]. In particular, it is known that the tracer gradient norm has a log-normal point-wise law. Some authors have also added a deterministic advecting velocity to the delta-correlated random term (e.g. [START_REF] Klyatskin | Short-time correlation approximations for diffusing tracers in random velocity fields: A functional approach[END_REF].

More recently, several authors have considered the random transport of fluid dynamics quantities by the combination of two velocity components: a time-correlated one and a time-uncorrelated one. Unlike the Kraichnan model, these stochastic subgrid parameterizations are not restricted to the linear dynamics of passive tracers and can consider large-scale time-correlated flows. In the same time, these new schemes still encompass skew-symmetric multiplicative noises. The mathematical theory relies on stochastic calculus and Ito-Wentzell formula [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF]. [START_REF] Brzeźniak | Stochastic partial differential equations and turbulence[END_REF] first introduces the idea. Then, [START_REF] Mikulevicius | Stochastic Navier-Stokes equations for turbulent flows[END_REF] and [START_REF] Flandoli | The interaction between noise and transport mechanisms in PDEs[END_REF] modify the formula. Their works have focused on pure mathematical aims: existence and uniqueness of SPDE solutions. [START_REF] Neves | Wellposedness for stochastic continuity equations with ladyzhenskaya-prodi-serrin condition[END_REF] also studied the wellposedness of similar SPDE using the derivation of the previous authors. Models under location uncertainty (LU) (Mémin, 2014;Resseguier et al., 2017a,b;Resseguier, 2017;[START_REF] Chandramouli | Deciphering the role of small-scale inhomogeneity on flow structuration: a stochastic approach[END_REF][START_REF] Cai | Motion estimation under location uncertainty for turbulent fluid flows[END_REF][START_REF] Resseguier | Stochastic modelling and diffusion modes for POD models and small-scale flow analysis[END_REF] and Stochastic Advection by Lie Transport (SALT) [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF][START_REF] Crisan | Solution properties of a 3d stochastic euler fluid equation[END_REF][START_REF] Gay-Balmaz | Stochastic geometric models with non-stationary spatial correlations in lagrangian fluid flows[END_REF]Cotter et al., 2018c,a,b;[START_REF] Resseguier | Stochastic modelling and diffusion modes for POD models and small-scale flow analysis[END_REF] are new types of stochastic subgrid tensors based on that stochastic transport. These schemes will be deeply detailed in section 6. [START_REF] Sardeshmukh | Reconciling non-gaussian climate statistics with linear dynamics[END_REF]; [START_REF] Sardeshmukh | Need for caution in interpreting extreme weather statistics[END_REF] also highlight the relevance of skew-symmetric multiplicative noises and its link with an advecting velocity time scale separation assumption. Besides, this separation assumption can be rigorously derived from homogenization technique [START_REF] Cotter | Stochastic partial differential fluid equations as a diffusive limit of deterministic lagrangian multi-time dynamics[END_REF]. With similar assumptions but different mathematical tools (Taylor series instead of stochastic calculus) Dukowicz and Smith (1997) obtain results for the mean tracers which coincide with LU framework. Nevertheless, equations for the mean is not sufficient for UQ purposes.

Memory effects

When there is no time-scale separation, the large-scale system can become non-Markovian. For deterministic dynamics, this is readily shown by the Mori-Zwanzig equation [START_REF] Givon | Extracting macroscopic dynamics: model problems and algorithms[END_REF][START_REF] Gottwald | Stochastic climate theory[END_REF]. Indeed, this explicit expression of the large-scale observables of interest involves a memory term. Accordingly, some authors have proposed non-Markovian stochastic parameterizations (e.g [START_REF] Chekroun | Predicting stochastic systems by noise sampling, and application to the el niño-southern oscillation[END_REF]Kondrashov et al., 2015;[START_REF] Kondrashov | Stochastic modeling of decadal variability in ocean gyres[END_REF][START_REF] Lu | Data-based stochastic model reduction for the kuramotosivashinsky equation[END_REF]. Even though these memory effects are realistic, it is important to note that their simulations imply an increase of the state-space dimension and hence an additional computational cost.

Other approaches based on statistical physics

Using Ruelle response theory, [START_REF] Wouters | Disentangling multi-level systems: averaging, correlations and memory[END_REF] have proposed a systematic way to approximate non-linear dynamical systems. Their main assumption is the weak coupling between two terms of the original evolution equation. In this approximated dynamics, a deterministic first-order correction appears. At second order, there is also a memory term and a stochastic forcing.

In quasi direct interaction closure (DIA), [START_REF] Frederiksen | Subgrid-scale parameterizations of eddy-topographic force, eddy viscosity, and stochastic backscatter for flow over topography[END_REF] proposed to simplify the equations for the mean flow and the fluctuations by assuming that the fluctuation two-time covariance function and the Green function of the diffusion equation (associated with eddy viscosity) are quasi-diagonal in Fourier space. We recall that a diagonal covariance in Fourier space means homogeneity in spatial space.One limitation of the method is the restriction to Fourier-based numerical simulations. [START_REF] Plant | A stochastic parameterization for deep convection based on equilibrium statistics[END_REF] have proposed a physically-based stochastic subgrid parameterization to model mass flux induced by clouds updrafts and downdrafts. Nevertheless, its application is a priori limited to weather forecasts.

Dynamics under Location Uncertainty and Stochastic Advection by Lie Transport

In this section, we focus on a new family of stochastic subgrid parameterizations. These schemes mainly rely on a stochastic transport. Therefore, in this section, we first explain the principle of this random transport and review its main properties. Most of these theoretical results come from Mémin (2014) and Resseguier et al. (2017b). The Stochastic Advection by Lie Transport can be introduced alternatively through stochastic differential geometry (e.g. [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF]. Then, we present the randomized Euler equations of this model family. Finally, we review the existing parameterization choices in this stochastic framework and propose a new one.

Stochastic transport

Informal description

CFD introduces de facto a coarse scale truncation of the system. Those dynamical models emanate from physical deterministic representations whose solutions are assumed to be smooth (i.e. differentiable) in time. Although small-scale fluid flow velocities can be characterized by local and intermittent energy bursts, possibly associated with infinitesimal characteristic timescales, it is generally assumed that these unresolved flow components remain smooth in time. At the model resolution, the resolved (large-scale) flow can thus be considered as a coarse-grained representation of the actual Eulerian flow, with the unresolved flow component rapidly varying in time. From an observer point of view, such sub-grid dynamics can be conveniently modeled by a delta-correlated process. The smooth velocity field, denoted w, represents a large-scale, possibly random, component continuous in time. The unresolved contribution, expressed as σ Ḃ, is then assumed Gaussian, volume preserving (divergence-free) and uncorrelated in time. This contribution can be non-homogeneous and anisotropic in space. Due to the irregularity of the resulting flow, the transport of a conserved quantity, q, by the whole velocity, defined as

q(X t+∆t , t + ∆t) = q(X t , t) (13)
reads in an informal stochastic way as

∂ t q + w • ∇q Corrected advection = ∇ • 1 2 a∇q Diffusion -σ Ḃ • ∇q Random forcing , (14)
with a drift velocity corrected as

w = w -1 2 (∇ • a) T + σ(∇ • σ) T . ( 15
)
Hence, the deterministic (for a fixed realization of velocity field w) evolution equation is replaced by a stochastic equation with respect to σ Ḃ. As a result, the conserved quantity is now advected by an "effective" velocity, w , taking into account the possible spatial variation of the small-scale velocity variance and the possible small-scale velocity divergence. This modified advection is indeed essential to take into account essential physical effects due the small-scale velocity inhomogeneity [START_REF] Chandramouli | Deciphering the role of small-scale inhomogeneity on flow structuration: a stochastic approach[END_REF][START_REF] Pinier | Stochastic flow approach to model the mean velocity profile of wall-bounded flows[END_REF]. The random forcing term in ( 14) relates to the advection by the unresolved velocity σ Ḃ = σ dBt dt . This term continuously backscatters random energy to the system. The diffusion term then accounts for the mixing effect of the small-scale random velocity. This term plays a role similar to the eddy diffusivity models introduced in classical large-scale representations [START_REF] Gent | Isopycnal mixing in ocean circulation models[END_REF][START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF][START_REF] Lilly | On the application of the eddy viscosity concept in the inertial subrange of turbulence[END_REF][START_REF] Smagorinsky | General circulation experiments with the primitive equation: I. the basic experiment[END_REF], in analogy with the molecular diffusion mechanism [START_REF] Boussinesq | Essai sur la théorie des eaux courantes[END_REF]. In particular, for a homogeneous and isotropic small-scale velocity the turbulent diffusion is also homogeneous and isotropic:

a = a 0 I d and ∇ • 1 2 a∇q = 1 2 a 0 ∆q. (16)
The inhomogeneous and anisotropic diffusion coefficient matrix, a, is then defined by the one-point one-time covariance of the unresolved displacement per unit of time:

a = E σdB t (σdB t ) T dt . ( 17
)
This ensures an exact energy balance between the amount of diffusion and the random forcing.

Stochastic flow

To derive more formally the evolution law of a scalar quantity transported by a stochastic flow, the stochastic Lagrangian description of the infinitesimal displacement associated with a particle trajectory X t writes:

dX t = w(X t , t)dt + σ(X t , t)dB t . (18)
In this equation, the second term explicitly figures the flow location uncertainty. Formally, this random field is defined over the fluid domain, Ω ⊂ R d , from a d-dimensional Brownian function B t . Such a function can be interpreted as a white noise process in space and a Brownian process in time.

Formally it is a cylindrical I d -Wiener process (see Da [START_REF] Da Prato | Stochastic Equations in Infinite Dimensions[END_REF][START_REF] Prévôt | A concise course on stochastic partial differential equations[END_REF], for more information on infinite dimensional Wiener process and cylindrical I d -Wiener process). The time derivative of the Brownian function, in a distribution sense, is informally denoted σ Ḃ = σ dBt dt , and is a white noise distribution. The spatial correlations of the flow uncertainty are specified through the diffusion operator σ(., t), defined for any vectorial function, f , through the matrix kernel σ(., ., t):

σ(x, t)f = Ω σ(x, z, t)f (z, t)dz. ( 19
)
This quantity is assumed to have a finite norm. More precisely, the operator σ is assumed to be Hilbert-Schmidt. We also assume that the above expression have periodic or null boundary conditions on the domain frontier. The resulting d-dimensional random field, σ(x, t)dB t , is a centered vectorial Gaussian function, correlated in space and uncorrelated in time with covariance tensor:

Cov(x, y, t, t ) = E (σ(x, t)dB t ) (σ(y, t )dB t ) T , (20)
= Ω σ(x, z, t) σT (y, z, t)dz δ(t -t )dt. ( 21
)
Hereafter, the diagonal of the covariance tensor, a, -also denoted as σσ Twill be referred to as the variance tensor:

a(x, t)δ(t -t )dt = Cov(x, x, t, t ).
Cotter et al. ( 2017) have rigorously shown that the decomposition (18) corresponds to the limit of a deterministic flow when the correlation time of the small-scale velocity goes to zero.

Scalar advection

For a fluid flow defined by equation ( 18), the material derivative (expressed in Ito form) of a quantity q writes:

D t q = (q(X t+dt , t + dt) -q(X t , t)) | X t =x = (22) d t q = q(x,t+dt)-q(x,t) Time increment + (w dt + σdB t ) • ∇q Advection -∇ • 1 2 a∇q Turbulent diffusion dt + tr ((σ T ∇) H T ) Coupling turbulence-forcing dt,
where the quantity q is forced as follows

D t q = F dt + H T dB t . ( 23
)
The time increment d t q can be interpreted as the analog of the partial time derivative ∂ t q in deterministic partial differential equations. The above expression -derived by Resseguier et al. (2017a) -is a reformulation of the Ito-Wentzell formula (theorem 3.3.1 page 91 [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF][START_REF] Chow | Stochastic partial differential equations[END_REF].

Stratonovich notations:

Equations ( 22) and ( 23) rely on the so-called Ito notations. They are convenient to derive ensemble mean -since terms in dB t are always centered -and for the numerical implementations. From these notations, we also get more insights on the involved physical processes (e.g. the turbulent diffusion). Nevertheless, it is also possible to write the very same equations in Stratonovich notations (i.e. with "•dB t " instead of "dB t "):

D S t q = d t q = q(x,t+dt/2) -q(x,t-dt/2) Centered time increment + (w S dt + σ • dB t ) • ∇q, (24)
where

D S t q = q X t+ dt 2 , t + dt 2 -q X t- dt 2 , t -dt 2 | X t =x , ( 25
) = D t q -1 2 tr ((σ T ∇) H T ) dt, ( 26
) = F -1 2 d dt H T , B t -1 2 tr ((σ T ∇) H T ) dt + H T • dB t , (27)
is the Stratonovich material derivative and

w S = w -1 2 d i=1 ∂ i σσ T i• = w * -1 2 σ (∇ • σ) T , (28)
is the Stratonovich drift of the flow (18). The Stratonovich material derivative, D S t q, coincides with the Ito material derivative, D t q, for H = 0, i.e. for time correlated forcing in the transport equation ( 23) -which is generally the case. The Stratonovich drift, w S , coincides with the corrected drift, w * , for divergence-free unresolved velocity -, which is the most common case. Equation ( 24) is the Ito-Wentzell formula for Stratonovich notations [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF][START_REF] Chow | Stochastic partial differential equations[END_REF]. The reader not familiar with stochastic partial differential equations can also refer to section 10.1.2 of Resseguier ( 2017) for more details. Equations ( 27) and ( 26) are derived in Appendix B by application of the Stratonovich-Ito-notation-change formula (Kunita, 1997, theorem 3.2.5 page 60).

Energy Energy conservation:

To ensure a stochastic isochoric flow, incompressibility constraints on the modified drift, ∇ • w = 0 and on the small-scale velocity, ∇ • σ = 0, are required. One can show (Resseguier et al., 2017a) that those two constraints enable to establish a strong energy conservation property for any realizations of a tracer (i.e. with zero forcing in ( 23)):

d dt Ω 1 2 q 2 = 0. ( 29
)
The noise energy intake is exactly compensated by the diffusion term. In particular, the mean energy, E Ω 1 2 q 2 , is also conserved. So, the LU-SALT also meets the balance (4) of section 5.3.

2 (with ζ = 1, L[q] = ∇ • 1 2 a∇q and (H[q]) 2 = H[q] 2 = 1
dt ∇q T a∇q) even though the SALT-LU random forcing is not of the form (8) in general.

The energy conservation (29) also implies that the decrease of the energy of the mean always goes with a variance increase:

d dt Ω 1 2 Var(q) = - d dt Ω 1 2 (E(q)) 2 , ( 30
)
This process is useful for UQ. While the tracer interacts with the unresolved scales, the tracer is continuously randomized.

Local energy in spatial space: Since D t is a derivative (considering Lagrangian coordinates), it is straightforward that:

D t 1 2 q 2 = qD t q = 0. ( 31
)
Hence, all the discussions about the tracer q are also valid for the local energy

1 2 q 2 .
Local energy in Fourier space (spectrum):

Here, q(k, t) = Ω q(x, t)e -ik•x dx stands for the Fourier transform of q. In order to get some insights on where the uncertainty goes in Fourier space, we can derive the time evolution of the tracer spectrum Γ (k, t) = E|q(k, t)| 2 . Indeed, if the small-scale velocity is homogeneous and isotropic in space, we can show (see appendix A) that:

∂ t Γ (k, t) = 2 (2π) d E w • ∇q (k, t)q(k, t) Usual spectrum variation + R d F σ (k, k , t)dk Noise-induced spectrum flux , ( 32
)
where {z} is the real part of z and the spectrum flux -induced by the small-scale velocity -from the wave-vector kk toward the wave-vector k is:

F σ (k, k , t) = k 2 (2π) d Γ σ (k , t) sin 2 ((k, k ))Γ (k -k , t) - 1 d Γ (k, t) , (33)
where

Γ σ (k , t) = 1 dt E σdB t (k , t) 2 . ( 34
)
The first term of the budget (32) encompasses the usual large-scale non-linear effects. The second term is the tracer spectrum increase which is due to the small-scale velocity. This energy flux does not depends on the tracer phases, but only on the tracer spectrum. Furthermore, the tracer mode q(k -k , t) gives energy to q(k, t), by the intermediate of σdB t if and only if:

Γ σ (k , t) = 0 and |q(k, t)| |q(k -k , t)| < √ d sin((k, k )) . (35)
In this case, the energy transfer from |q(k -k , t)| 2 toward |q(k, t)| 2 is mostly a random process, since it is due to the positive difference between the noise energy intake -which is a random process -and the turbulent dissipation.

Moreover, if we assume that the (bidirectional or tridirectional) spectrum is decreasing with the wave number k and that the dimension d = 2, a backscattering of energy from the small-scale Fourier mode q(k -k , t) to the large-scale Fourier mode

|q(k, t)| > |q(k -k , t)| implies: cos π 2 -(k, k > 1 √ 2 = cos π 4 (36) i.e. k, k modulo π = -k, (k -k ) modulo π ∈ π 4 , 3π 4 . ( 37
)
Note that this results holds even if the statistics of the small-scale velocity depends on the tracer q, as long as the small-scale velocity is homogeneous and isotropic.

Moments of passive tracers:

Consider now that the expectation corresponds to a conditional expectation given the effective drift. This applies to passive scalar transport for which the drift does not depend on the tracer. In equation ( 22), terms in dB t have zero-mean, and the mean passive scalar evolution can be immediately derived taking the conditional expectation of the stochastic transport:

∂ t E(q) + w • ∇E(q) Advection = ∇ • 1 2 a∇E(q) Diffusion . ( 38
)
Since w * is divergent-free, it has no influence on the energy budget. The mean field energy decreases with time due to diffusion. As for the variance, its evolution equation, derived in Resseguier et al. (2017a), reads:

∂ t V ar(q) + w • ∇V ar(q) Advection = ∇ • 1 2 a∇V ar(q) Diffusion + (∇E(q)) T a∇E(q) Variance intake . (39)
This is also an advection-diffusion equation, with an additional source term.

Integrating this equation on the whole domain, with the divergent-free condition, and considering the divergence form of the first right-hand term, we

obtain d dt Ω V ar(q) = Ω (∇E(q)) T a∇E(q) 0. ( 40
)
It shows that the stochastic transport of a passive scalar creates variance. The dissipation that occurs in the mean-field energy equation is exactly compensated by a variance increase. This mechanism is very relevant for ensemblebased simulations. The uncertainty modeling directly incorporates a large-scale dissipating sub-grid tensor, and further encompasses a variance increase mechanism to balance the total energy dissipation. Such a mechanism is absent in ensemble-based data assimilation development [START_REF] Berner | Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations[END_REF][START_REF] Gottwald | The role of additive and multiplicative noise in filtering complex dynamical systems[END_REF][START_REF] Snyder | Performance bounds for particle filters using the optimal proposal[END_REF].

Extensive properties : Reynolds-transport theorem

Similar to the deterministic case, the stochastic Reynolds transport theorem shall describe the time differential of a scalar function, q(x, t), integrated over a material volume, V(t), transported by the random flow (18): d

V(t) q = V(t) D t q + ∇ • (w dt + σdB t ) q + d t 0 D t q, t 0 ∇ • σdB t .
(41) This expression, rigorously derived in Resseguier et al. (2017a), was first introduced in a slightly different version by Mémin (2014). In most cases, the unresolved velocity component, σ Ḃ, is divergence-free and, the source of variations of the extensive property V(t) q is time-differentiable, i.e. with a differential of the form d V(t) q = Fdt. In such a case, for an arbitrary volume, the transport theorem takes the form D t q = f dt, and the material derivative can be replaced by the equation ( 22) to provide an intrinsic expression of this stochastic transport theorem.

Euler models

We will now express the incompressible Euler models existing in this family of random schemes. The Euler model under Location Uncertainty (LU) differs from Stochastic Advection by Lie transport (SALT) Euler model. For sake of simplicity and in contrast to the rest of the paper, we present them only in their Stratonovich forms.

Additionally, in this review, we will not considered the compressible case, since it brings other complexities. Accordingly, the density ρ is assumed to be constant. Using the mass conservation, this implies the following divergencefree constraints (Mémin, 2014;Resseguier et al., 2017a):

∇ • w * = 0, (42) ∇ • σ = 0. ( 43
)
After presenting SALT and LU Euler models, we will present a deterministic LES-like model refers to as pseudo-stochastic model. It has also been referred to as Euler model under Location Uncertainty (Mémin, 2014;[START_REF] Resseguier | Geophysical flows under location uncertainty, part I: Random transport and general models[END_REF][START_REF] Harouna | Stochastic representation of the reynolds transport theorem: revisiting large-scale modeling[END_REF]Resseguier et al., 2017a,c;Chandramouli et al., 2018;Cintolesi and Mémin, 2019a,b). Therefore, we briefly explain it here for disambiguation.

Euler LU

The Euler LU assumes stochastic transport of the Ito drift w up to some forcings:

D t w = D S t w Due to the transport of ρw = g -1 ρ ∇p dt. ( 44
)
The forcings can also encompass other terms (e.g. viscous terms, time-uncorrelated pressure forcing). By the properties of the stochastic transport, the Navier-Stokes LU conserves the (Ito) kinetic energy, Ω 1 2 w 2 (chapter 8 of Resseguier, 2017).

SALT Euler

In contrast, the SALT Euler is derived from differential geometry arguments with a Kelvin circulation preservation constraint and reads:

D S t w * Due to the transport of ρw * + d j=1 ∇(σ • dB t ) j w * j Additional term = g -1 ρ ∇p dt. ( 45
)
The reader may refer to [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF] for a full presentation of the derivation. The SALT Euler does not conserve kinetic energy (chapter 10 of Resseguier, 2017) but conserves the (Stratonovich) helicity, H = Ω w * • ∇ × w * [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF][START_REF] Resseguier | Stochastic modelling and diffusion modes for POD models and small-scale flow analysis[END_REF]. The advantage of this model is that it leads to a classical form of the vorticity equation. However, this at the price of an additional term in the momentum equation and a loss of energy conservation. [START_REF] Crisan | Solution properties of a 3d stochastic euler fluid equation[END_REF] have proofed the local well-posedness of the 3D incompressible SALT Euler model. In this review, we will not consider the compressible case, since it brings other complexities.

Pseudo-stochastic Euler LU

The pseudo-stochastic Euler model under location uncertainty writes:

∂ t w + (w • ∇) w -∇ • 1 2 a∇w T T = g -1 ρ ∇p. ( 46
)
Compared to the (stochastic) Euler model under location uncertainty ( 44), the pseudo-stochastic Euler model under location uncertainty (46) does not involve the stochastic forcing (σdB t • ∇) w. This is because it is assumed that the drift w is a smooth function of time (more precisely a finite variation process) or, equivalently, that an uncorrelated pressure forcing compensates the stochastic forcing -(σdB t • ∇) w -, which ensues from the unique decomposition in terms of martingales and finite variation terms. Accordingly, this pseudo-stochastic model is deterministic and hence easier to simulate for LES-like applications. The pseudo-stochastic model still retains the action of the unresolved small-scale components in terms of dissipation and modified advection, but does not include the backscaterring effect brought by the random component. As a consequence, for UQ purposes, the (stochastic) Euler model under location uncertainty (44) may be more suitable (Resseguier et al., 2017a). Note that this pseudo-stochastic model has been successively used to define efficient reduced order models (Resseguier et al., 2017c), in which the subgrid dissipation is directly defined from the neglected modes. This representation has the advantage to provide also new diagnosis enabling to quantify local energy dissipation as well as the effect of the turbulence inhomogeneity on the large-scale flow [START_REF] Pinier | Stochastic flow approach to model the mean velocity profile of wall-bounded flows[END_REF]Resseguier et al., 2017c). This ability is a strong asset of the LU formalism for reduced order modeling.

Parameterization of LU-SALT models

LU-SALT models are parameterized by the unresolved velocity spatial covariance (20). It can be an anisotropic, heterogeneous, time-dependant or even be a function of the large-scale quantities (e.g. q). Choosing this covariance is equivalent to choosing the linear operator σ or choosing its kernel σ. In this section, we review the different parameterizations that have been considered so far. At the end of the section, we also propose a new one. It enables in particular time-dependant heterogeneous unresolved velocity statistics without relying on off-line high-resolution simulation outputs.

Kraichnan model

Similarities between LU-SALT dynamics and the Kraichnan model (Kraichnan, 1968, 1994;[START_REF] Gawȩdzki | Anomalous scaling of the passive scalar[END_REF][START_REF] Majda | Models for stochastic climate prediction[END_REF] were already highlighted by Mémin (2014). [START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF] focuses on the passive scalar transport by a homogeneous (Cov(x, y, t, t ) = Cov(xy, t, t )), isotropic (Cov(x, y, t, t ) = Cov( xy , t, t )) and divergence-free velocity field uncorrelated in time. With LU notations, it implies:

w = 0, (47) ∇ • σ = 0, ( 48
) σdB t = σ * dB t , (49)
where * denotes a convolution and σ is a particular (isotropic) filter. The homogeneity condition (49) makes the parameterization easier in Fourier space:

E σdB t (k 1 , t) σdB t T (k 2 , t ) = Γ σ (k 1 )δ(k 1 -k 2 )δ(t -t )dt, ( 50
)
where Γ σ can be taken to mimic the self-similarities observed in real turbulence flow [START_REF] Gawȩdzki | Anomalous scaling of the passive scalar[END_REF]:

Γ σ (k) = A 2 P (k) 1 + k 2 κ 2 m -2α , ( 51
)
with A is a constant, 1/κ m specifies the velocity correlation length and

P (k) = I d - kk T k 2 , (52)
is the projection -expressed in Fourier space -onto the space of solenoidal functions. The Fourier transform of the associated filter σ (appearing in (49)) is:

σ (k) = A P (k) 1 + k 2 κ 2 m -α . ( 53
)
Finally, the Fourier transform of the small-scale velocity writes:

σ Ḃ (k) = A √ ∆t P (k) 1 + k 2 κ 2 m -α dBt √ ∆t (k), (54)
where dB t is the spatial Fourier transform of dB t , with dB t √ ∆t, a vector of d independent discrete white noise processes of unit variance in space and time. To sample the small-scale velocity, we first sample dB t √ ∆t, to get dB t √ ∆t, and finally σ Ḃ (k) with the above equation. Note that the same unresolved simulation is possible without relying on Fourier transform. Indeed, one can use the fact that

M (k) = 1 + k 2 κ 2 m -α
is the Fourier transform of a Matérn covariance M of degree ν = α -1 and of range ρ = 2(α -1)/κ m . Therefore, the unresolved velocity is:

σdB t = AP {M * dB t } , (55)
with

P = I d -∆ -1 ∇∇ T .
(56)

The homogeneous stationary model

Assuming zero large-scale velocity is obviously not possible in CFD. However, combined with a resolved large-scale velocity, a homogeneous and stationary small-scale velocity can be a relatively good approximation for some flows.

Note that even a homogeneous small-scale vecolity leads to a heterogeneous non-Gaussian random forcing, due to the its multiplicative structure. Resseguier et al. (2017b) show numerical simulations of a LU version of a 2-dimensional geophysical model referred to as the Surface Quasi-Geostrophic (SQG) model [START_REF] Blumen | Uniform potential vorticity flow: part I. theory of wave interactions and two-dimensional turbulence[END_REF][START_REF] Held | Surface quasi-geostrophic dynamics[END_REF][START_REF] Lapeyre | Surface quasi-geostrophy[END_REF]. For this purpose, a modification of ( 53) is considered. Instead of involving the operator Resseguier et al. (2017b) work on the streamfunction ψ σ dB t . As such no additional divergence-free constraint is needed. The streamfunction kernel, ψσ is defined such as:

P = I d -∆ -1 ∇∇ T ,
σ(x)dB t = ∇ ⊥ ψ σ (x)dB t = ∇ ⊥ ψσ dB t (x). ( 57
)
with

∇ ⊥ = (-∂ y , ∂ x ) T .
Accordingly, a single cylindrical Wiener process, B t , is sufficient to sample the Gaussian process. This is specific to two-dimensional domains. In 3D, a vector of 3 independent I d -cylindrical Wiener processes, and a curl must be considered to simulate an isotropic small-scale velocity. Thus, the small-scale velocity can be conveniently specified from its omnidirectional spectrum:

k → 1 µ(Ω) E [0,2π] dθ k k σ Ḃ (k) 2 = 2π ∆t k 3 ψσ ( k ) 2 , (58)
where µ(Ω) is the surface of the spatial domain Ω, θ k is the angle of the wavevector k and ∆t the simulation time-step. Consistent with SQG turbulence, the omni-directional spectrum slope, denoted s, is fixed to -5/3. For 2D Euler equations, the slope would be set to -3. The unresolved velocity should be energetic only where the dynamics cannot be properly resolved. Consequently, the authors apply to the spectrum a smooth band-pass filter, f BP , which has non-zero values between two wavenumbers κ m and κ M . The parameter κ m is inversely related to the spatial correlation length of the unresolved component. κ M is set to the theoretical resolution, π/∆x, and κ m to the estimated effective resolution (κ m = κ M /2 in Resseguier et al. (2017b)). The small scales' energy is specified by the diffusion coefficient a and the simulation time step:

E σ Ḃ σ Ḃ T = 1 ∆t a = 1 ∆t a 0 0 0 a 0 . ( 59
)
The diagonal structure of the variance tensor is due both to incompressiblity and isotropy. The scalar variance tensor, a 0 , is similar to an eddy viscosity coefficient. So, a typical value of eddy viscosity used in practice is a good proxy to setup this parameter. Otherwise, this parameter can be tuned. The time step depends itself, through the CFL conditions, on both the spatial resolution and the maximum magnitude of the resolved velocity. Finally, similarly to the Kraichnan model, the Fourier transform of the small-scale velocity writes:

σ Ḃ (k) = A √ ∆t ik ⊥ f BP ( k ) k -α dBt √ ∆t (k) with s = 3 -2α = - 5 3 , ( 60
)
where 2017b) have numerically demonstrated the good UQ skills of this method for SQG flows.

A is a constant to ensure E σ Ḃ 2 = tr(a)/(∆t) = 2a 0 /(

The homogeneous non-stationary model

A main drawback of the previous choice of σ is the needed tuning of the parameters. [START_REF] Resseguier | Stochastic modelling and diffusion modes for POD models and small-scale flow analysis[END_REF] propose a more general homogeneous parameterization where no tuning is needed. A new σ will be defined at each time step from the resolved velocity kinetic energy spectrum.

In the previous parameterization, the absolute diffusivity (i.e. KE times correlation time [START_REF] Falkovich | Particles and fields in fluid turbulence[END_REF]Penland, 2003b;[START_REF] Klyatskin | Stochastic equations through the eye of the physicist: Basic concepts, exact results and asymptotic approximations[END_REF][START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF][START_REF] Keating | Diagnosing lateral mixing in the upper ocean with virtual tracers: Spatial and temporal resolution dependence[END_REF])) of the unresolved velocity is twice the variance tensor trace T (a) whereas the unresolved kinetic energy is T (a)/dt. Clearly, this kinetic energy has no physical meaning. Indeed, it depends on the simulation time step and one should have the possibility to choose the time step as close as possible to zero. Thus, the unresolved velocity amplitude is specified through an absolute diffusivity rather than through a KE. In the mathematics literature of homogenization, Kubo-type formulas may be seen as what physicists call absolute diffusivities. More generally, since the variance of a time-continuous white noise is infinite, it is more relevant to deal with absolute diffusivity rather than kinetic energy in order to describe the statistics of the time-uncorrelated velocity. Thus, keeping a spectral approach, we define -for any spatio-temporal field -an Absolute Diffusivity Spectral Density (ADSD) denoted A(κ) at the wave-number κ. We will rely on this ADSD rather than on the KE spectrum, E(κ) . Since the absolute diffusivity is the variance multiplied by the correlation time, it is naturally to defined the ADSD as follows:

A(κ) = E(κ)τ (κ) where τ (k) = 1/κ v κ = 1/κ κE(κ) , (61)
is the eddy turnover time at the scale 1/κ and v κ is at characteristic velocity a this scale. Accordingly, we have:

A(κ) = κ -3/2 E 1/2 (κ). ( 62
)
If in addition we assume a KE self-similar distribution,

E(κ) = C 2 κ -s , (63)
we obtain:

A(κ) = Cκ -r , (64)
where r = (s + 3)/2.

We aim at defining the unresolved velocity ADSD from the large-scale velocity. For this purpose, we will assume the self similar model ( 64) is valid at all spatial scales. At each time step, we compute the ADSD of the largescale velocity, A w , from formula (62). Then, we fit the coefficients C and r of equation ( 64). Let us denote with C w and r w these coefficients. Note that they are time-dependent because they depend on w which is. More precisely, we estimate the coefficients C w and r w in a wavenumber interval which approximately represents a inertial range of fully-resolved scales (i.e. before the spectrum roll-off).

From there, we can define the operator σ in such a way that the total velocity -resolved plus unresolved -meets (64) at small spatial scales with the same coefficients (C w , r w ). Since the two velocity components are not correlated, the total ADSD is the sum of the ADSD of each velocity component. Therefore, σ is chosen such as the unresolved ADSD compensates the resolved ADSD roll off -introduced by the deterministic subgrid parameterizationat small scales. Specifically, the unresolved ADSD is set to:

A σ Ḃ (κ) = max(0, C w κ -rw -A w (κ)) f 2 BP (κ). ( 65
)
As previously, f 2 BP is a band-pass filter between κ m and κ M . In practice, we set κ M to the theoretical resolution:

κ M = π ∆x , (66)
and κ m to the effective resolution which is estimated as follows:

κ m = ln(0.95) ln(0.1) 1/p κ M , ( 67
)
where p is the order of the Laplacian used as deterministic subgrid tensor (i.e. D t b = -ν(-∆) p b dt). The justification of the above formula is left in the Appendix C. Compared to the work of Resseguier et al. (2017b), the value of κ m is less critical. Indeed, equation ( 65) implies a weaker unresolved ADSD at larger scales where the resolved ADSD, A w is stronger. This softens the threshold effect introduced by the band-pass filter f BP .

In practice, we set an upper-bound for the estimation of r w . Without this upper bound, a concentration of energy at relatively large wave-numbersscales smaller than κ m -in the resolved fields can become unstable. Indeed, this localized energy concentration would decrease the r w estimation, and hence increase the unresolved ADSD A σ Ḃ through (65) at large wave-numberslarger than κ m . This implies a larger noise intake, which can induce a larger concentration of energy at relatively large wave-numbers in the resolved fields, resulting in a positive feedback loop. To prevent these unphysical instabilities, the slope r w is bounded.

In order to link the unresolved ADSD to the kernel σ which defines the unresolved velocity, we note that

A σ Ḃ (κ) = E σ Ḃ (κ)dt = 1 µ(Ω) [0,2π] dθ k κ σ(t, κ) 2 = 2πκ 3 ψσ (t, κ) 2 . ( 68
)
From formulas ( 65)-( 68), we can finally express the unresolved velocity as follow:

σ Ḃ(t, k) = 1 √ ∆t ik ⊥ ψσ (t, k ) dB t √ ∆t (k), (69)
= 1 √ ∆t ik ⊥ max(0, C w k -rw -A w ( k )) 2π k 3 f BP ( k ) dB t √ ∆t (k). ( 70
)
Again the simulated unresolved velocity ADSD is physically relevant while the KE spectrum is not. Indeed, the simulated unresolved velocity ADSD is expected to match the true (time-correlated) unresolved velocity ADSD, whereas the KE spectra of the simulated and true unresolved velocities differ. Indeed, the true unresolved velocity correlation time spectral distribution τ (κ) is not restricted to the time step dt.

Resseguier et al. ( 2019) have successfully applied this parameterization for UQ purpose in a SQG flow.

Heterogeneous modulation of homogeneous models

The SALT-LU stochastic parameterization randomly folds tracer isolines. This process is often desirable. For instance, it can trigger physically-relevant instabilities, such as the filament instabilities highlighted by Resseguier et al. (2017b). After these instabilities have been randomly triggered, eddies are formed by non-linear processes. However, a homogeneous small-scale velocity may also perturb the tracer isolines which should remain still (i.e. which remain still in high-resolution deterministic simulations), e.g. sharp, straight, coherent fronts. A typical application of this problem in more realistic flow simulations is the simulation of jets like the Gulf Stream and regions of the Antarctic Circumpolar Current. These real-world jets are associated with diffusivity suppression [START_REF] Ferrari | Suppression of eddy diffusivity across jets in the Southern Ocean[END_REF], and this effect is not present in our formulation so far. If we seek to preferentially perturb some tracer gradients, a heterogeneous small-scale velocity is required. Note that the heterogeneity discussed here needs to be non-stationary and thus cannot be represented by the stationary EOF presented later in this paper. Besides, in a small ensemble of realizations, relevant heterogeneity of the small scales may make the spreading more accurate for UQ and enable comparable ensemble forecast accuracy with fewer members. We here propose a possible heterogeneous version of the previous method.

In order to obtain a heterogeneous model of the unresolved velocity, we need a heterogeneous version of the ADSD (64). Since the wave-number κ cannot depends of the position x, the constant C w and/or the spectrum slope r should do. A spatially varying spectrum slope is probably difficult to estimate. Hence, we restrict ourselves to a spatially-varying constant C w and a spatiallyhomogeneous spectrum slope. The constant may also varies with the time and the wave number. According to the Kolmogorov theory (e.g. [START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF] and ( 63):

C w (x, t, κ) = cst. p F (x, κ), (71)
where F is the energy flux through the spatial scales and p = 1/3 for a SQG flow. More specifically, the energy flux describes the energy moving from scales larges than 1/κ toward smaller scales and can be computed as follows:

F (x, t, κ) = q < κ ((w • ∇)q) < κ , ( 72
)
where q is the transported (up to possible source terms) quantity, g < κ is the low-pass filtered version of g (setting to zero the Fourier modes of g which have frequencies larger than κ) and • stands for the spatial average [START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF]. For a SQG flow, q corresponds to the buoyancy normalized by the stratification: q = b/N . The energy flux is essentially a third-order moment. It is very important because it describes the cascade of the flow by non-linear energy transfers [START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF].

If the energy flux through scale is understood locally in space (as indeed Smagorinsky (1963) also assumes), the formula (71) provides a natural parameterization of the unresolved velocity heterogeneities. We simply modulates the unresolved ADSD (65) by the heterogeneous ratio p F / p F -averaged over the resolved inertial range wave-numbers:

(73) σ Ḃ(t, k) = 1 √ ∆t ik ⊥ max(0, C w k -rw -A w ( k )) 2π k 3 f BP ( k ) dB t √ ∆t (k), σ Ḃ(t, x) = P p F (t, x) p F (t)
Heterogeneous modulation σ Ḃ(t, x)

Homogeneous velocity ,

where P = I d -∆ -1 ∇∇ T is the projector onto the space of free-divergence functions. This parameterization is physically meaningful since, locally in space, a stronger direct cascade at large scales (larger F and thus larger C w ) suggests that the unresolved velocity (large) should maintains this cascade by folding smaller-scale tracer structures. Furthermore, considering that the energy flux is a third-order structure makes this parameterization relevant to differentiate between strait fronts and curved structures (e.g. eddies). Indeed, at least three points are needed to define a curvature and differentiate between these structures.

In order to keep a divergence-free velocity and the ensuing properties (e.g. energy conservation), the modulated velocity is projected onto the space of free-divergence functions, using the operator P. Because of that we do not consider the advection correction w * -w = -1 2 (∇ • a) T of the LU formalism. Indeed, here the variance tensor has the simple form a = 1 d tr(a)I d . As such, the advection correction is a gradient field and is hence removed by the projection onto the space of free-divergence functions. [START_REF] Resseguier | Stochastic modelling and diffusion modes for POD models and small-scale flow analysis[END_REF] have numerically shown that this parameterization indeed keeps the sharp straight fronts while perturbing the eddies and the meanders in SQG flows.

The heterogeneous spatially-uncorrelated model

Yang and Mémin (2018) have considered a spatially heterogeneous subgrid velocity. The variance tensor a is estimated on-line by data assimilation. To simplify, they neglect the spatial correlation of that velocity.

Off-line learning of the EOF from velocity

It is always possible to write the unresolved velocity as follows (see e.g. Mémin, 2014):

σ(x, t) Ḃ = n KL k=1 ξ k (x, t) Ẇk , (75)
with W k are independent one-dimensional Brownian motions and n KL possibly infinite. The (ξ k ) k are so-called (weighted) Empirical Orthogonal Functions (EOFs). If we denotes by λ k the L 2 norm of the EOFs, their normalized ver-

sions ξk = ξ k /λ k k
are the eigenvectors of the auto-adjoint operator defined by the kernel σ(x, t)σ T (y, t). The (λ k ) k are the corresponding eigenvalues.

That kernel corresponds -up to the factor 1 dt -to the one-time two-point covariance of the unresolved velocity. In the case of stationary EOFs, the type of decomposition ( 75) is supported by the mathematical theory of σσ T -Wiener processes [START_REF] Da Prato | Stochastic Equations in Infinite Dimensions[END_REF][START_REF] Prévôt | A concise course on stochastic partial differential equations[END_REF].

First, we will consider the off-line learning of the EOFs (ξ k ) k from a velocity field set. Specifically, we will assume that we have a history of "small-scale" velocity fields

(v (x ij , m)) 1 m n O at n O different times t m (or n O different "realizations") sampled on a spatial grid of M points (x ij ) ij = ((i∆x, j∆y)) ij .
Typically, the velocity snapshots are high-resolution simulation outputs or high-pass-filtered versions of it. We also need to assume that the EOFs are stationary in time. For the sake of clarity, we restrict the methods presentation to d = 2. Notwithstanding, the generalizations to three-dimensional flow are straightforward.

Candidate for the increment realizations

In order to estimate the EOFs, (ξ i ) i , we assume that we can observe increments

tm+∆T tm σdB t = n KL k=1 ξ k (x ij )(W k (t m + ∆T ) -W k (t m )). ( 76
)
at several times t m . Here, ∆T is the time step of a low-resolution simulation. We will interpret the following snapshots:

∆X m ij = v (x ij , m)∆T, (77)
as a realizations of the flow increments (76). Since Brownian time increments (76) are independent, by the above interpretation, we implicitly assume that v (x ij , m) are independent for different values of m. In this case, it is equivalent to consider them as as a set of independent "realizations" or as a set of snapshots at several times t m .

Preprocessing

The increments are supposed to be centered and divergence free (∇ • ξ i = 0). Therefore, after computing the residual flow increments ∆X

m ij m
, they are centered:

∆X m ij = ∆X m ij -Ê{ ∆X m ij }, (78)
with the estimator

Ê{ ∆X m ij } = 1 n O n O m=1 ∆X m ij , (79)
and projected onto the space of divergence-free functions:

∆X m ij = P{∆X m ij } with P = I d -∆ -1 ∇∇ T . ( 80
)
Note that the projection can also be applied directly on the final EOFs ξ kat the end of the estimation process -rather than on each increments ∆X m ij . Besides, if ∇ • v = 0 then the centred flow increments are already divergencefree and the step (80) can be skipped.

Covariance and EOF

Then, we can define the EOFs from an estimate of the spatial covariance of the residual flow increments (averaging over the time index m) :

γ σ (x ij , x pq ) ≈ 1 (n O -1)∆T n O m=1 ∆X m ij (∆X m pq ) T . (81)
Besides, we have:

γ σ (x ij , x pq ) = n KL k=1 ξ k (x ij ) ξ k T (x pq ). (82)
In order to properly define the EOFs, we must add the following orthogonal constraint:

pg ξ i (x pq ) • ξ j (x pq ) = 0 if i = j. ( 83
)
According to equations ( 82) and ( 83), the EOFs can be obtained from a diagonalization of the covariance (γ σ (x ij , x pq )) (ij),(pq) .

Dual problem for the EOF estimation

Since the number of grid points M is often larger than the number of increments, n O , the covariance (81) can be a very large matrix (M d × M d).

It is easier to consider the dual problem to estimate the EOFs. As such, the temporal covariance is first computed:

γσ (t p , t q ) = 1 M ij 1 ∆T ∆X p ij • ∆X q ij . (84)
After diagonalizing this n O × n O matrix, the EOFs are obtained by projecting the eigenvectors onto the normalized increments 1

√ ∆T ∆X k ij .
Finally, after estimating the (ξ k ) k off-line, the ensemble forecast can be generated on-line with the formula (75). This techniques has recently been applied for the stochastic simulation of turbulent channel flow. These stochastic simulation show how the inhomogeneity of the small-scales velocity components allows structuring the large-scale component in terms of streaks [START_REF] Chandramouli | Deciphering the role of small-scale inhomogeneity on flow structuration: a stochastic approach[END_REF].

Off-line learning of EOF from flows

Now, we detail a procedure proposed by Cotter et al. (2018c,a,b) for the estimation of the EOFs, (ξ k ) k , involved in the unresolved velocity definition (75).

Here again, the EOFs are assumed stationary in time. The main difference with the off-line learning of EOF from velocity 6.3.6 is the candidates which are used for the increments ∆X m ij . The data-driven method of Cotter et al. (2018c,a,b) relies on Lagrangian paths defined at two "resolutions". The first paragraph of this section defines these two types of Lagrangian paths. Then, we explain how to obtain the candidates for the flow increment realizations from these Lagrangian paths.

Preliminary definitions

We introduce two types of velocity field:

a high-resolution velocity v on a fine spatial mesh-grid, a low-resolution velocity v on a coarse spatial mesh-grid. This velocity field is a spatially-low-pass-filtered version of v.

Then, two types of Lagrangian path are defined:

a "high-resolution flow", X ij (t 0 , t), defined by the high-resolution velocity v:

dX ij dt (t) = v (t, X ij (t 0 , t)) and X ij (t 0 , t 0 ) = x ij . ( 85
)
a "low-resolution flow", X ij (t 0 , t), defined by the low-resolution velocity v:

dX ij dt (t 0 , t) = v t, X ij (t 0 , t) and X ij (t 0 , t 0 ) = x ij . (86)

Candidate for the increment realizations

In order to estimate the EOFs, we will interpret the following residual flow increments as realizations of the unresolved stochastic flow increments (76):

∆X m ij = X ij (t m , t m + ∆T ) -X ij (t m , t m + ∆T ). ( 87
)
Then, the following steps are the same than in 6.3.6 (preprocessing and covariance diagonalization).

This method (Cotter et al., 2018c,a) may seem more natural than the one of 6.3.6, since the increments (87) are really flow increments. Nevertheless, the time step ∆t being small, ∆X m ij is close to (v -v) (t m , x ij ) ∆T . Therefore, the main difference between the two methods is probably in the definition of the "residual velocities": v on one hand and (v -v) on the other hand.

The method 6.3.7 has been successfully tested on a 2D Euler flow (Cotter et al., 2018c), a quasi-geostrophic flow (Cotter et al., 2018a) and a SQG flow [START_REF] Resseguier | Stochastic modelling and diffusion modes for POD models and small-scale flow analysis[END_REF]. [START_REF] Resseguier | Stochastic modelling and diffusion modes for POD models and small-scale flow analysis[END_REF] also compare the data-driven method 6.3.7 and the homogeneous non-data driven method 6.3.3 for forced homogeneous SQG turbulence. Similarly good UQ results are obtained.

Off-line learning of the stochastic advection operator

The accuracy of the unresolved velocity Karhunen-Loeve decomposition ( 75) is influenced by the number of kept modes EOFs, n KL . In fluid mechanics, the state-space dimension is often huge. Therefore, that kind of spectral decomposition of the covariance often involves a number of modes much smaller than the state-space dimension. Hence, the decomposition (75) strongly reduces the sampling computational cost since there is no need to manipulate the full spatial covariance.

Nevertheless, for specific applications like dimensionally-reduced models, the state space can be reduced to a small well-chosen vector space (e.g. [START_REF] Holmes | Turbulence, coherence structures, dynamical systems and symetry[END_REF]. By definition, the resolved velocity component lies on this subspace, whereas the unresolved component does not (Mémin, 2014;[START_REF] Resseguier | Geophysical flows under location uncertainty, part I: Random transport and general models[END_REF]Resseguier et al., , 2017c)). If the state-space dimension -says n r -is small enough, it can be more interesting to directly work with the statistics of the random operator:

f → σdB t • ∇f, (88)
expressed in the reduced basis, rather than dealing with the unresolved velocity statistics. Indeed, the reduced representation of the stochastic operator is a (conditionally-) Gaussian matrix of size n r × n r . It is centered and time-uncorrelated. Therefore, the probability law of this matrix is exactly and entirely determined by its covariance, which have n 2 r × n 2 r coefficients. In comparison, the approximate unresolved velocity Karhunen-Loeve decomposition would lead to stochastic operator

f → n KL k=1 (ξ k • ∇f ) Ẇk , (89)
characterized by n 2 r × n KL coefficients. For n KL < n 2 r , this last representation is necessary an approximation (since the covariance of the exact representation of (89) in the reduced subspace is of full rank, i.e. of rank n 2 r ). For small n r , one can rely on the random operator reduced representation instead. Rigorous and low-cost estimators exist for this representation covariance matrix coefficients and provide good UQ skills (Resseguier, 2017, chapter 8).

Nevertheless, for large state-space dimension n r , this n 2 r × n 2 r covariance matrix of the random operator is too large to handle. Hence, additional assumptions on the unresolved small-scale velocity (e.g. homogeneity, parametric model, transport) or a Karhunen-Loeve decomposition become necessary. 6.3.9 Transport of the unresolved velocity Methods 6.3.6 and 6.3.7 learn the unresolved small-scale structure from highresolution simulations. However, such simulations are not always available. Moreover, methods 6.3.6 and 6.3.7 assume heterogeneity without time dependence. While some turbulence heterogeneities could be fixed in time (e.g. heterogeneity due to boundaries conditions), we can expect that a large part of it moves with the flow.

Accordingly, Gay-Balmaz and Holm (2018) assume a "transport equation" -in the sense of differential geometry -for the EOFs:

D S t ξ k = (ξ k • ∇)(w * dt + σ • dB t ). ( 90
)
As such, the unresolved velocity is both heterogeneous and non-stationary in time, without requiring any learning. However, as far as we know, this method has not been tested in numerical simulations yet.

On-line learning of EOF from velocity

We propose here a new approach -based on the EOF decomposition (75) -but where the EOFs are time-dependent. They are estimated on-line from a coarse simulation. For this purpose, we propose to generate pseudo-observations of the small scales directly from the resolved velocity. Then, we will compute an EOF representation as in 6.3.6. Finally, we will adapt the noise variance to smaller scales using a turbulence power-law scaling (originally proposed by [START_REF] Harouna | Stochastic representation of the reynolds transport theorem: revisiting large-scale modeling[END_REF].

Pseudo-observations

The approach proposed in this section defines n O pseudo-observations v at each simulation grid point before computing the singular value decomposition (SVD).

For a given time t and a given resolved low-resolution velocity component, w, we build pseudo-observations by sliding a n w × n w (n w odd) window over the spatial grid. We denote by L = n w l the spatial scale of the window, where l is the smallest scale of the simulation. At every grid point x ij , we list the n 2 w velocity values contained in the window centered at that point -with appropriate boundary conditions (replication, periodicity, etc. ) when looking at a point on the border:

I(x ij , t) = w(x pq , t) |p -i| n-1 2 , |q -j| n-1 2 . ( 91
)
Then, for each m ∈ {1, ..., n O }, for each point x ij independently, we set the value of the pseudo-observation v (x ij , t, m) by randomly choosing a value in the set I(x ij , t).

After this, we proceed as in 6.3.6 to obtain the EOFs, but averaging over the pseudo-realization index m instead of averaging over the time steps t m .

Rescaling

From the SVD, we obtain a set of EOFs ξ (L) k and hence a model for the unresolved velocity σ (L) (x, t) Ḃ:

σ (L) Ḃ = n KL k=1 ξ (L) k Ẇk . ( 92
)
The pseudo-realizations v = v L have been generated at a spatial scale L = nl. These fluctuations correspond to a virtual observation scale L (n × n window) and must be scaled down to the "simulation scale" l. Therefore the unresolved velocity variance tensor, a, is rescaled by a coefficient proposed by [START_REF] Harouna | Stochastic representation of the reynolds transport theorem: revisiting large-scale modeling[END_REF], and which in 2D reads:

a l = l L 2/3 a L , (93)
where a L and a l are the variance tensors at the scales L and l respectively. This scaling relies on the Kolmogorov-Richardson cascade assumption with the velocity fluctuations at scale as u ∝ 1/3 1/3 . The unresolved velocity, σ Ḃ, can be finally simulated -at the "simulation scale" l -as:

σ Ḃ = n -1/3 σ (L) Ḃ. ( 94
)
Let us note that in 3D the scaling must be adapted as the power exponent is 4/3 [START_REF] Harouna | Stochastic representation of the reynolds transport theorem: revisiting large-scale modeling[END_REF]. The on-line SVD method 6.3.10 will be numerically tested in section 8. There, the LU parameterizations 6.3.2, 6.3.3 and 6.3.10 will be compared to the new scheme 5.3.2 and to ensembles generated by initial condition perturbations.

Metrics for UQ

In this section, we detail some important metrics to quantify ensemble forecasts prediction skills.

This section is a short summary about existing verification ensemble tools. More detailed review could be found in [START_REF] Hamill | Verification of eta-rsm short-range ensemble forecasts[END_REF][START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF][START_REF] Scheuerer | Variogram-based proper scoring rules for probabilistic forecasts of multivariate quantities[END_REF][START_REF] Thorarinsdottir | Verification: assessment of calibration and Trevisan A, Palatella L (2011) On the Kalman filter error covariance collapse into the unstable subspace[END_REF].

Talagrand diagram

A Talagrand diagram (or rank histogram) is a technique used to check the reliability of an ensemble forecast or a set of quantile. The idea is considers the availability of N observation q o 1 , ..., q o N in a set of rank r. The sets of rank are built from rank statistic which is a sorted ensemble {(q (i1) 1 , ..., q (in e ) 1

), ..., (q

(i1) N , ..., q (in e ) N
)} of n e members [START_REF] Anderson | A method for producing and evaluating probabilistic forecasts from ensemble model integrations[END_REF][START_REF] Talagrand | Evaluation of probabilistic prediction systems[END_REF][START_REF] Hamill | Verification of eta-rsm short-range ensemble forecasts[END_REF]:

∀(j, k) ∈ {1, ..., n e } × {1, ..., N }, r j = 1 N N k=1 P (q (ij-1) k ≤ q o k < q (ij ) k ), ( 95
)
where P is the estimated probability of ranking an observation between two sorted ensemble members. In the context of (95), N represents the spatial points or times steps of the observation.

A calibrated ensemble should result in a flat histogram. However, a flat histogram does not guarantee a calibrated ensemble [START_REF] Hamill | Interpretation of rank histograms for verifying ensemble forecasts[END_REF]. A flat histogram mostly indicates that the ensemble and observation are sampled from a common distribution. A U-shaped rank histogram corresponds to an underdispersion or to a conditional bias [START_REF] Hamill | Interpretation of rank histograms for verifying ensemble forecasts[END_REF] of the ensemble members. A dome-shaped shows an overdispersion of the ensemble. A nonsymmetric histogram is the footprint of bias.

Mean absolute and mean squared error

The mean absolute error (MAE) and mean squared error (MSE) are the stateof-the-art verification and selection tools specially used in regression problems as decision function. In an ensemble verification problem, the MAE is defined by the following formulation:

∀(j, k) ∈ {1, ..., n e } × {1, ..., N }, MAE = 1 N n e N k=1 ne j=1 |q o k -q (j) k |, (96)
and the MSE by this one:

∀(j, k) ∈ {1, ..., n e } × {1, ..., N }, MSE = 1 N n e N k=1 ne j=1 (q o k -q (j) k ) 2 . ( 97
)
The goal behind these two equations is to calculate the total error between a sample of observation and ensemble of size N . So, the empirical mean of all error between an observation q o k and each member q

(j)
k of the ensemble of size n e is calculated. Both metrics are proper scores. A proper score is negatively oriented, such that a lower score indicated a better ensemble forecast. A proper score converging to zero means that the observation cannot be distinguished from the ensemble members. There are several proper scoring rules.

Continuous ranked proper score

Another proper score used for ensemble model evaluation is the continuous ranked proper score (CRPS). This kind of score is characterized by the socalled predictive distribution F -the distribution represented by the ensemble forecast {q (1) , ..., q (ne) } realizations and the observation q o . The CRPS is defined by the following equation [START_REF] Matheson | Scoring rules for continuous probability distributions[END_REF][START_REF] Hersbach | Decomposition of the continuous ranked probability score for ensemble prediction systems[END_REF]:

CRPS(F, q o ) = R (F (B) -H(B -q o )) 2 dB, ( 98
)
where H is the heaviside function (equals to one if q o ≤ B and zero otherwise), F is the cumulative distribution function (CDF) of the random process q (represented by the ensemble) at one spatial point x and time t. The CRPS is the distance between the heaviside function and the distribution of the random variable q. This heaviside function represents the inequality between an ensemble member and the observation. If F has a finite first-order moment, there is another representation of the CRPS [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF]:

CRPS(F, q o ) = E F |q -q o | -1 2 E F |q -q |, (99)
where q and q are independent realizations of F . If a priori assumption are made on the behaviour of the random variable q. It is possible to deduce an analytic expression of the CRPS thanks to the equation ( 98) and ( 99).

In the case where F is issued from a normal distribution N (µ, σ 2 ), the CRPS formula is [START_REF] Gneiting | Calibrated probabilistic forecasting using ensemble model output statistics and minimum crps estimation[END_REF]:

CRPS(F, q o ) = σ ω(2Φ(ω) -1) + 2φ(ω) - 1 √ π , (100)
where ω = (q o -µ) σ , φ is probability density function (PDF) of the standard normal distribution and Φ is the CDF of the standard normal distribution.

Also following [START_REF] Grimit | The continuous ranked probability score for circular variables and its application to mesoscale forecast ensemble verification[END_REF], for an ensemble forecast {q (1) , ..., q (ne) }, a natural fair estimator of the CRPS is given by:

CRPS (F, q o ) = 1 n e ne i=1 |q (i) -q o | - 1 2n 2 e ne i=1 ne j=1 |q (i) -q (j) |. (101)
The expression (101) allows to compute the CRPS without making a priori assumptions on our random variable. However, this expression has a O(n 2 e ) computational complexity.

An algebraically equivalent representation of the CRPS based on the generalized quantile function is proposed by [START_REF] Laio | Verification tools for probabilistic forecasts of continuous hydrological variables[END_REF])

CRPS ( F , q o ) = 2 n 2 e ne j=1 (q (ij ) -q o )(n e H(q (ij ) -q o ) -j + 1 2 ) , (102)
where {q (i1) , ..., q (in e ) } is the sorted ensemble forecast. The computational complexity of the expression ( 102) is O(n e log(n e )). The interested reader can refer to [START_REF] Jordan | Facets of forecast evaluation[END_REF] for mathematical details on the discrete form of the CRPS.

Energy score

In a CFD problem, observations and ensembles evolve in space and time. The CRPS is a pointwise score, hence cannot summarize the whole uncertainty of the ensemble in one value in a multivariate case. Therefore, to perform the evaluation of our ensemble in a multivariate case, the CRPS will not be enough. Of course, the means of the CRPS can be use. However, we would miss some information. Instead, the energy score can be considered. There, our observation q o and our random variable q is multivariate. The corresponding ensemble {(q (1) 1 , ..., q

(1) N ), ..., (q (ne) 1 , .., .q (ne) N )} contains n e realizations. Each of them represents the field of interest in N time steps or spatial points. F is now the multivariate CDF associated to q. We need a generalized form of the CRPS to work in this case. Following [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF][START_REF] Scheuerer | Variogram-based proper scoring rules for probabilistic forecasts of multivariate quantities[END_REF], a possible generalization of the CRPS is the energy score (ES):

ES(F, q o ) = E F q -q o -1 2 E F q -q , ( 103
)
where . is the euclidean norm. So, an estimator of the ES would be:

ES ( F , q o ) = 1 n e ne i=1 q (i) -q o - 1 2n 2 e ne i=1 ne j=1
q (i) -q (j) (104)

Variogram score -p

The energy score can provide an evaluation of the uncertainty of a multivariate ensemble. Nevertheless, this score is often not sufficiently sensitive to misspecification of the correlations structures of the multivariate observation and the multivariate ensemble. Those correlations cannot be study by univariate proper score unlike the means and variance. Thus, it does not evaluate the capacity of the ensemble model to reproduce the correlations structures in time and/or space of the observation.

One tool used to study these correlations is the variogram of order p (also called structure function). This tool consists in studying the pairwise differences of the component of the multivariate ensemble forecast and observation by considering the statistics of our multivariate observation. Following [START_REF] Scheuerer | Variogram-based proper scoring rules for probabilistic forecasts of multivariate quantities[END_REF], the p-variogram score (VS-p) is an alternative solution which overcomes this drawback. We express it using the same multivariate random variable q and observation q o :

V S p (F, q o ) = N i=1 N j=1 1 w ij (|q o i -q o j | p -E F |q i -q j | p ) 2 , ( 105
)
where the w ij are non negative weights that allows one to emphasize or downweight pairs of component combinations based on subjective expert decisions. p is the order of the variogram score. N represents the number of spatial points or time steps. F is the multivariate cumulative distribution function of q. According to [START_REF] Scheuerer | Variogram-based proper scoring rules for probabilistic forecasts of multivariate quantities[END_REF], a typical choice of weight to discriminate time or space dependencies is the inverse distance weight:

w ij := 1 √ |i-j|
and p := 0.5. The VS-p score measures the dissimilarity between approximations of the variograms of order p of the observations and forecast over all pairs of the components of our variable of interest.

Numerical results : Ensemble forecasts verification for SQG dynamics

In this section, we will present new numerical UQ results. As a test case, a simple geophysical fluid dynamics -the stochastic surface quasi-geostrophic (SQG) model -will be considered. After presenting it, several stochastic subgrid parameterizations will be compared for short-term and for long-term ensemble forecasts. We will focus on SALT-LU schemes, the new random forcing derived from numerical dissipation of section 5.3.2 and random perturbation of initial condition.

Test flows

Deterministic and stochastic surface quasi-geostrophic (SQG) dynamics

From here, the field q = b will denote the buoyancy. It is proportional to the density anomaly

ρ b = -g ρ ρ 0 with ρ(x, y, z, t) = ρ 0 1 - N 2 g z + ρ (x, y, z, t), ( 106
)
where ρ is the density and N the stratification. In the ocean, the density anomaly is small compared to the total density (Boussinesq approximation) and the flow is approximately isochoric. The conservation of salinity and temperature, with a linearized equation of state, provides the transport of buoyancy [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF]. Then, considering rapid rotation, strong stratification and uniform potential vorticity leads to the so-called SQG model [START_REF] Blumen | Uniform potential vorticity flow: part I. theory of wave interactions and two-dimensional turbulence[END_REF][START_REF] Held | Surface quasi-geostrophic dynamics[END_REF][START_REF] Lapeyre | Surface quasi-geostrophy[END_REF].

A random version of this model (denoted SQG M U ) can be derived from the location uncertainty principle (Resseguier et al., 2017b). It will keep the same structure except that the buoyancy is now transported in the stochastic sense ( 14). The horizontal velocity u = w is related to the buoyancy b in Fourier space through the usual SQG relation:

û = ik ⊥ b N k , ( 107
)
where k is the horizontal wave-vector. The unresolved velocity, σ Ḃ, is also horizontal. Consequently, the variance tensor, a, is a 2 × 2 matrix and the SQG M U model is two-dimensional.

Our SQG simulations

High-resolution deterministic SQG simulations of test flows will provide references to which we will compare random simulations performed at a lower resolution. For this purpose, the high-resolution ( 5122 ) simulation outputs will be projected onto the space of low-resolution (128 2 ) fields, i.e. it will be addequately filtered and subsampled to a low-resolution. We will refer to this projected fields as "observations" and denote them b o . For all simulations, a standard hyperviscosity (HV) scheme has been introduced [START_REF] Held | Surface quasi-geostrophic dynamics[END_REF]:

D t b = α∆ 4 b dt, (108)
with a positive coefficient α proportional to M -8

x where M x denotes the grid size (i.e. 128 or 512). The domain size is a square box L x ×L y = 1000 km×1000 km and the boundary conditions are doubly periodic.

Several simulations of these models have been performed based on two types of initial fields and on various paraterisations of uncertainty. A first type of initial buoyancy field is shown in figure 2 (a). As in Resseguier et al. (2017b), this field consists of a spatially smooth buoyancy field with two warm elliptical anticyclones (positive buoyancy) and two cold elliptical cyclones (negative buoyancy) given by:

b 0 (x) = F x - 250 km 250 km + F x - 750 km 250 km -F x - 250 km 750 km -F x - 750 km 750 km , (109)
with

F (x) = B 0 exp - 1 2 x 2 σ 2 x + y 2 σ 2 y and σ x = 67 km, σ y = 133 km. ( 110
)
The amplitude of the buoyancy and the stratification are set to: B 0 = 10 -3 m.s -2 and N = 3f 0 . The Coriolis frequency f 0 is fixed to 1.028 × 10 -4 s -1 , which corresponds to a latitude of 45 • .

Another type of initial buoyancy field is shown in figure 2(b), which is a homogeneous Gaussian random field generated from a spectrum of buoyancy with a prescribed slope equals to -5/3. This slope corresponds to the power law of a developped SQG turbulence. We will now forecast the initial buoyancy fields of figure 2 with serveral ensembles of simulations and assess the performance of each of these ensembles. Short term and long term forecasts will be treated separately because they are very different in nature.

Short-term ensemble forecast

As in the study Resseguier et al. (2017b), we first focus on the first-month forecast of the smooth initial field (figure 2(a)). Figure 3 shows the reference high-resolution simulation for this first month. After two weeks, filament instabilities [START_REF] Lapeyre | Surface quasi-geostrophy[END_REF] create developed turbulence. Here, we also study the free-decaying SQG turbulence flow initialed by the rough field (figure 2(b)). The free-decaying turbulence can be seen in figure 4. A part of the initial energy is dissipated by the deterministic subgrid tensor. But, a part of this initial energy cascades to the larger scales by creating larger vortices from the merging of small vortices. Each ensemble cointains 30 realizations.

Several low-resolution ensembles have been forecast. Two ensembles rely on initial conditions random perturbations, three on dynamics under location uncertainty and five on the dissipation-adapted noise (see section 5.3.2). Specifically, the small-scale velocity in the dynamics under location uncertainty and the random perturbation of initial condition in SQG P IC have been both generated using three approches: a homogeneous stationary model (denoted "Spectral"), a homogeneous non-stationary model (denoted "ADSD") and a heterogeneous non-stationary model (denoted "SVDpseudo"). For the model under location uncertainty, it corresponds to parameterizations 6.3.2, 6.3.3 and 6.3.10 respectively, of section 6.3. For the "Spectral" method, the small-scale energy is specified by the diffusion coefficient a 2 = 9 m 2 .s -1 . We also need to choose the spectrum slope and the minimum and maximum wavelengths of the unresolved velocity. For more, information on the values of these parameters, the reader can refer to Resseguier et al. (2017b). In contrast, for the improved "ADSD" method, we do not need to choose any parameter. The parameterization must adapt itself using large-scale informations at each time step. Despite the difference of flow and resolution, the very same code -without any modification -is used here and in [START_REF] Resseguier | Stochastic modelling and diffusion modes for POD models and small-scale flow analysis[END_REF] to define the "ADSD" small-scale velocity. For the ensembles with perturbed initial conditions, similar sampling methods are adopted. But the obtained "small-scale" random fields are used only at the initial time. As an example, figure 5 shows the initial perturbation spectrum used in the "Spectral" method for the first initial condition (109). Then, it is added to the initial condition b(x, t = 0). Besides, the dissipation-adapted noise 5. | f ( κ)

| 2 κ( r a d . m -1 )

D e t e r mi n i st i c an d r an d om p ar t s of t h e i n i t i al c on d i t i on

Fig. 5 Kinetic energy spectrum (m 2 .s -4 /(rad.m -1 )), at the initial time, of the mean buoyancy, in blue, spectrum of its random perturbation, in red, and slope -5 3 in black. The initial perturbation is restricted to a narrow spectral band. This random initial condition has been used to simulate an ensemble with the deterministic SQG model.

variance contains a larger part of the numerical dissipation), more and more small-scale structures are presented in physical field.

Once the ensembles have been produced by the previous models, we try to measure the quality of ensemble forecasts by some easy-to-implement criterions. We have N e ensemble members for each forecast with b (i) denotes the i-th ensemble member of the forecast. A first necessary condition for reliability is that the mean squared bias (MSB) of the ensemble (i.e. the MSE of the ensemble mean) is close to the mean intra-ensemble variance (MEV), up to an ensemble size-dependant scaling factor :

M SB ≈ N e + 1 N e M EV, (111)
where

M SB = 1 M M j=1 b o -E b 2 (t, x j ), (112)
M EV = 1 M M j=1 V ar (b)(t, x j ), (113)
M the number of grid points x j , E {f } and V ar (f ) the empirical mean and the empirical variance of f , computed from the ensemble, respectively. At the same time, a classical error metric -the ensemble mean square error (MSE) is also considered : Fig. 6 Comparison of one realisation of the buoyancy field after 17 days of advection, according to various energy scaling factor in "WavHypervis".

MSE(t) = E   1 M M j=1 (b o -b) 2 (t, x j )   . (114)
In figures 7(a) and 7(b), we compare the criterions ( 111) and ( 114) produced by our models. For the 4 vortices field 2(a), it seems that for all models produce slight errors during the first 10 days. Their forecast skills tend to a stationary state after about 17 days' advections. In this case, WavHypervis and SVDpseudo under location uncertainty have a better ensemble spread, yet conversely the homogeneous models -Spectral and ADSD have lower MSE. As illustrated in figures 7(c) and 7(d), in the WavHypervis method, when the factor ζ increases, both the spreading and the ensemble errors increase. In figures 7(e) and 7(f ), we compare the results with the spectral initial field 2(b), in which the ensembles models tend to a stationary state much more faster. In this case, WavHypervis and SVDpseudo under location uncertainty have still better ensemble spread. The forecasts by Spectral and ADSD consist of less errors in ensemble.

Another intuitive estimation of ensembles dispersions have been performed by Talagrand diagram. As shown in table 1, under the 4 vortices background, the small perturbations method applied on the initial condition is too underdispersive. Instead, Spectral and ADSD under location uncertainty formed a ∪-shape obtaining underdispersion and so higher but still underestimated variance. SVDpseudo under location uncertainty estimates high uncertainty. Instead, WavHypervis with 50% energy allocated from the numerical dissipation provides a slight overdispersion at last. This is probably due to a premature bifurcation in this ensemble (see section 8.3 for more details about the bifurcation). From table 2, with the spectral initial field, WavHypervis, and SVDpseudo under initial perturbations provide almost perfect ensemble spread. The ensemble skills converge toward a stationary state much faster for this flow.

As explained above, distinct models yield distinct spreading based criterions ( 111 type provide opposite conclusions. An objective analysis need hence others metrics like proper score (see section 7).

To begin, we will compare the numerical results with the proper score CRPS (see the section 7.3). The CRPS is defined as a pointwise score. Applied on ensemble of spatio-temporal fields, maps of CRPS can be represented at each fixed time step as in the figure 8. As expected on all maps of the figure 8, the normalized CRPS is relatively high on small scale structures and low on the center of each vortices. Indeed, turbulent structures are the hardest things to reproduced. Moreover, the figure 8 suggests that the model Spectral in SQG P IC is the worst model since it has the largest CRPS values. Model WavHypervis in SQG M U shows smaller CRPS. Thus reveals that this model makes fewer local errors at this advection time. The CRPS has only been computed at one fixed time. So, again, it is difficult to conclude on which model is better.

To push further the analysis, we maps at different times. The figure 9 reveals a part of the spatial evolution of those CRPS maps for the model Spectral in SQG M U . With the random gaussian field as initial condition, we can see on the figure 9 that the CRPS structures are first small and spread over the spatial domain. Then, these small scale CRPS structures merge and create larger structures of larger intensity. The merging is due to the inverse energy cascade of SQG turbulence [START_REF] Lapeyre | Surface quasi-geostrophy[END_REF] which merges the badly resolved small-scale turbulence structures. On the top of that, the new large vortices have chaotic trajectories. The difficulty for the models to track these trajectories yield large CRPS values in the centers of those vortices. Our analysis describes the variability of the CRPS in time and space. But, the difficulty to select a model remains to summarize the CRPS information.

The figure 10 summarizes model performances by representing each CRPS by one box plot. As such, we can directly see the variability of CRPS values over one CRPS map. For the four vortices initial condition, the model Spectral in SQG P IC has the most spread CRPS. Moreover, at 30 days of advection, it has the largest CRPS mean. Therefore, the figure 10 confirms that the model Spectral in SQG P IC has the worst uncertainties quantification skill. The models Spectral in SQG M U and ADSD in SQG M U obtain the lowest CRPS mean and spread at 20 and 30 days of advection. However, it is difficult to choose between these two models. Their CRPS distributions are too close. On the left panel of the figure 10 spectral background, it is easier to select a model. Indeed, at 10 days of advection, the model Spectral in SQG P IC gets the lowest CRPS mean and at 20 days the model ADSD in SQG M U has the lowest mean. Nevertheless, it is still difficult to conclude about the most efficient model. Boxplots of CRPS give statistical description about this proper score. But, it does not analyze the multivariate structure of the error between the ensemble and the reference.

The figure 11 plots the normalized energy scores (see 7.4) of advection times. There, the analysed multivariate structure is the spatial structure of the random fields. The figure 11 shows that the model Spectral in SQG P IC has the largest energy score for both initial conditions. So, this model provides the worst uncertainties quantification with multivariate spatially ensemble.

The model SVDpseudo in SQG P IC has a low energy score after two weeks of advection for the spectral background and after three weeks for the fourvortice initial condition. With the spectral background, after two weeks largescale structures has formed and the model SVDpseudo in SQG P IC makes less error with the reference. Nevertheless, during the first week of advection, the model SVDpseudo in SQG P IC has a really large energy score small-scale structures are present and are not well reproduced by SVDpseudo in SQG P IC . So, in this case, the models SQG P IC are less efficient than the models SQG M U and WavHypervis.

For the four-vortice initial condition, the model SVDpseudo in SQG P IC has a larger energy score at 16 to 20 days of advection. During this period of advection, some small-scale structures are produced and well resolved by the simulations. After this period, the energy score of SVDpseudo in SQG P IC decreases probably due to the bifurcation phenomenon described in section 8.3.

The figure 11 suggests that models SVDpseudo in SQG M U and WavHypervis make the least error. Indeed, during the most important times steps and for both initial condition, these models show the lowest energy scores. In particular, the energy score of the WavHypervis method reveals that it is the most efficient model according to this metric.

We will now consider proper scores and variograms where the temporal structures are considered instead of the spatial ones. The figure 12 reflects that the model Spectral in SQG P IC has the largest normalized energy score and variogram-0.5 score (see section 7.5). On the edges of the centers of the vortices, this model encompasses the worst temporal structures. The models SVDpseudo in SQG M U and WavHypervis seem to perform better. Nonetheless, in the model Spectral in SQG M U , minimum values are located in the center of the vortices which indicates that the temporal evolution of the vortices is well reproduced by this model.

Conclusion of the short-term forecasts analysis

For this numerical study, ensemble forecast skills have been assessed through a set of verification tools. The table 3 summarizes validation scores estimated for each UQ model. Talagrand histogram is an evaluation of the calibration behavior of an ensemble, while proper scores focus on dispersion and errors between univariate or multivariate ensembles and references. Each verification score has his own specificity, it is essential to be aware of their properties to select the most suitable measure (see section 7). Also, to avoid misselection of ensemble methods, the assessment of a number of scoring rules is advised.

In the table 3, SQG P IC methods obtained the lowest overall verification scores. In contrast, ADSD in SQG M U , SVDpseudo in SQG M U and WavHypervis models present the best performances. Therefore, we recommend one of this method for UQ tasks, and we strongly advice to avoid relying only on initial conditions randomization.

The model ADSD in SQG M U has the best overall MSE, meaning that it introduces less errors than other UQ methods. The model SVDpseudo in SQG M U and WavHypervis get lower MSE performances. Yet, theirs CRPS, ES and VS are the best among all the methods. This illustrates the accurate univariate and multivariate UQ produced by SVDpseudo in SQG M U and WavHypervis algorithms. WavHypervis model reveals a weaker MSB score analysis than SVDpseudo in SQG M U induced by a lack in the ensemble spread. For this reason, if one UQ method has to be chosen among ADSD in SQG M U , SVDpseudo in SQG M U and WavHypervis models, it would probably be SVDpseudo in SQG M U . Nevertheless, the higher UQ skills of this method also come with a slightly higher CPU time.

This numerical study have characterized UQ skills of short-term ensemble forecasts. To complete the analysis, we now focus on long-term forecast.

Long-term ensemble forecast

After the first month of advection, the flow initiated with (109) (figure 3) breaks its symmetries. This leads to a chaotic behavior and a complete loss of (deterministic) predictability. In this section, we will focus on this regime. As demonstrated here, the bifurcation associated with the symmetry breaking is efficiently tracked using the model under location uncertainty. Using few realizations, the probability density functions of each subsequent scenarios is well characterized, whereas methods based on random initial conditions do not converge. The identification of several scenarios, done at each time step in a reduced subspace, is obtained by Principal Component Analysis (PCA), also termed Empirical Orthogonal Functions method (EOF).

Finally, a diagnosis of bifurcation is performed and discussed from both the model under location uncertainty and a method based on randomized initial conditions.

Chaotic test flow and resolution issues

The boundaries conditions of the simulation are doubly periodic and, for the initial condition (109), there is a meridional line of symmetry at x = 500 km. Therefore, the zonal period of the initial condition is L x /2 = 500 km. This periodicity is relatively stable and holds during the first month. Nevertheless, the SQG dynamics, is subject to an inverse cascade of energy (Capet et al., 2008), and vortices of the same sign tend to merge. When this merging occurs, this affects the global shape of the flow. In particular, the periodicity that remains in the first month eventually disappears. This symmetry is hence metastable rather than stable. The symmetry breaking corresponds to a transition from one "state" to another. By "state", we mean a relatively "compact" and connected subspace of the state space. Warm vortices can merge at x = 0 or at x = 500km. In the following, we will refer to the first case as "scenario 1", and to the second case as "scenario 2". Because of the periodic boundary conditions, these two possible transitions are likely to occur. In a deterministic numerical simulation, the appearance of one transition or the other is determined by an infinitesimal asymmetry in the initial condition or possibly by a numerical error. This is a bifurcation. The bifurcation related to the merging of cold vortices is similar. With those two simultaneous bifurcations, there are thus 2 × 2 = 4 likely transitions.

To trigger a particular transition, we introduce two infinitesimal modifications in the initial condition. For the sake of simplicity, we focus on the bifurcation associated with the warm vortices. To do so, the merging of cold vortices in x = 500 km will be forced, by adding an infinitesimal small-scale cold eddy in (x, y) = (480, 750) (in km). The barycentre of northern structures becomes slightly closer to x = 500 km. This gives rise to the desired transition, as shown in Figures 13 and14. To trigger the bifurcation associated with the southern warm vortices, an infinitesimal small-scale eddy has been added in (x, y) = (20, 250) (in km). If the eddy is warm, the southern barycentre is moved closer to x = 0, and the two warm vortices merge near x = 0 (scenario 1), as shown in Figure 13. If this eddy is cold, the southern barycentre is moved closer to x = 500 km, and the two warm vortices merge near x = 500 km (scenario 2), as shown in Figure 14

where the large-scale field, b 0 , and the two-dimensional Gaussian function, F , remain defined by ( 109) and ( 110). The factor s w is set to 1 (respectively -1) if one wants to force the scenario 1 (respectively the scenario 2). The size of the large vortices of b 0 is of the order of the Rossby radius L d , whereas the small-scale eddies spread only over few kilometers.

The high-resolution simulations corresponding to scenario 1 and 2 are displayed in Figures 13 and14, respectively. At this resolution, the evolution toward scenario 1 or scenario 2 is determined by the value of the parameter s w (i.e.+/ -1). The associated variations are hardly taken into account by the low-resolved SQG model, as shown in Figure 15. The two SQG simulations correspond to 70 days of advection at high and low resolution, respectively. Both simulations have been initialized in the same way (s w = 1). The lowresolution field differs from the high-resolution field, as it apparently followed the wrong transition. As understood, whatever the random or deterministic nature of the tracer evolution law, a dissipation or a filtering at small scales is necessary to remove aliasing effects. In present simulations, the dissipation is created by an hyperviscosity scheme. At low resolution, the initial perturbation is rapidly diffused (few days). When the symmetry breaking occurs, after 40 days of advection, this initial perturbation has been completely forgotten. Moreover, another infinitesimal asymmetry triggers the other likely transition. According to the expression of b 0 (equations ( 109) and ( 110)), the large scale of the initial condition is not exactly zonally periodic with period 500 km. The southern part is slightly warmer in the middle. Indeed, the value of b 0 on (x, y) = (500, 250) (in km) is about 1.8 × 10 -5 m.s -2 , and on (x, y) = (0, 250) (in km) is about 8.8×10 -6 m.s -2 . The initial barycentre of the southern structures is thus closer to x = 500 km. This asymmetry has a very weak amplitude but a large spatial length scale which prevents its diffusion. This explains the merging in the wrong location.

In the next section, we will show that low-resolution simulations of the dynamics under location uncertainty can retrieve the right scenario.

Stochastic analysis

Unlike the deterministic SQG model, the SQG dynamics under location uncertainty, with the exact same initialization s w = 1, yields several likely transitions. Here, for simplicity, only the variant Spectral of SQG M U is considered (neither ADSD nor SVDpseudo). The calibration Spectral of SQG M U is detailed in section 6.3.2. In Figure 16, we show two realizations of the SQG M U dynamics. One of those realizations corresponds to the reference scenario (scenario 1), the other does not. The model encodes several likely transitions, and Here, a cold and a warm very small eddies are added respectively in the top and the bottom of the initial condition. These eddies are highlighted by respectively a blue square and a red square. This small difference in the initial condition does not modify the flow until the onset of the symmetry breaking, the 40 th day. Since the flow is chaotic, the small perturbation at t = 0 determines how the symmetry breaking occurs a month and a half later.

thus several potential scenarios. Indeed, the random forcing provides various small-scale perturbations that may trigger these transitions. As this triggering is random, the large-scale changes are also random. In other words, there is a backscattering of uncertainty toward the large scales, as illustrated in Figure 17. We decomposed the mean omni-directional spectrum, i.e. the mean energy at a given scale, Ê{Γ b }, into the spectrum of the mean tracer, Γ Ê{b} , (blue (shaded grey):

Ê{Γ b } = Ê 1 µ(Ω) k 2 | b| 2 dθ k , (116)
= 1 µ(Ω) k 2 | Ê{ b}| 2
Energy of the mean

+ Ê | b -Ê{ b}| 2 Variance dθ k , ( 117
) = Γ Ê{b} + Ê Γ b-Ê{b} , (118)
where µ(Ω) denotes the area of the domain Ω, Ê{f } the empirical mean of f , computed from the ensemble, and f the Fourier transform of f . The more thick the shaded grey area is, the more variance is contained at this scale. At t = 30 days, the variance remains at small scales but this small-scale uncertainty activates the symmetry breaking. This results in a variance backscattering with a thickening of the spectrum of the random component at large scales (visualized at t = 70). Since large scales influence strongly the small scales the small-scale variance is also enforced.

One can then wonder about the number of realizations following the right scenario. This necessitates the analysis of the 4-dimensional spatio-temporal random field b (i) (x j , y k , t l ) ijkl . The superscript (i) designates the i-th realization of the ensemble. To reduce the associated dimension, a Principal Component Analysis (PCA) -also termed Empirical Orthogonal Function (EOF) representation -is performed over the realizations, at a fixed time t. Within this analysis, unlike usual EOF representation, the time axis is replaced by the realization index. At a given time, it helps to represent the whole ensemble (200 realizations of the random field) by the ensemble mean field and few other

4 pdf of 1 st EOF coefficient pdf of EOF coefficients -2 0 2 1 st EOF coefficient ×10 -4 -3 -2 -1 0 1 2 3 2 nd EOF coefficient ×10 -4
After the transition

pdf of EOF coefficients -2 -1 0 1 2 1 st EOF coefficient ×10 -4 -2 -1 0 1 2 2 nd EOF coefficient ×10 -4
Fig. 18 Energy of the EOFs (m 2 .s -4 ) (left), probability density function of the first EOF coefficient (middle) and joint probability density function for the two first EOF coefficients (right) after 30 days (top) and 70 days (bottom). The ensemble is simulated according to the SQG model under moderate uncertainty at resolution 128 2 . At t = 30 days, the energy of the mean field, denoted as a 0-th order EOF (index+1 = 1), is much higher than the variance field, described by the other EOFs. On the contrary, at t = 70 days (after the symmetry breaking), the energy of the mean and variance have the same order of magnitude. The variance is mainly explained by the first EOF. Therefore, at the first order, the random component of the tracer can be approximated by this first EOF. The probability density function of the first EOF coefficient and the joint probability density function for the two first EOF coefficients are unimodal at t = 30 days and bimodal at t = 70 days. The breaking symmetry has created two likely scenarios, which are very different from one another. The scenario A corresponds to negative values of the first EOF coefficient (probability of 47%) whereas the scenario B corresponds to positive value (probability of 53%). The red line separates the probability density function between the two scenarios.

B. The red line at zero separates the probability density function between the two scenarios. In Figure 19, the same probability density function along time is plotted. The bifurcation is clearly visible. Also shown, the mean buoyancy (m.s -2 ) of the two likely scenarios are represented after 70 days of advection.

The two mean fields, Ê(b|1) and Ê(b|2), are calculated as:

Ê(b|1) = Ê(b) + Ê(c 1 |1)Ψ 1 = Ê(b) +    1 # i|c (i) 1 < 0 c (i) 1 <0 c (i) 1    Ψ 1 , (124) Ê(b|2) = Ê(b) + Ê(c 1 |2)Ψ 1 = Ê(b) +    1 # i|c (i) 1 0 c (i) 1 0 c (i) 1    Ψ 1 , (125)
where # stands for the cardinality of a set. The scenario A is quite close to the scenario 1, which is the reference with this initial condition, whereas scenario B is close to scenario 2. The stochastic model has enabled the ensemble to track both scenarios and to describe them statistically. Let us point out that the shape of the isotropic small-scale velocity expression has been loosely fixed by an a priori form of the spectrum. Some learning procedures of the noise topology from past data could lead to express more informative heterogeneous random fields, and to statistically favor the most likely transition. Moreover, since the two scenarios are very different, the introduction of few observations, through an assimilation procedure, could very easily help to select the right scenario.

For sake of comparison, we also show results obtained using the deterministic SQG model, with the "Spectral" perturbation of the initial conditions. This method was already used in the previous section for short-term forecast. We first briefly recall how these perturbations are generated. Initial smallscale buoyancy perturbations are assumed Gaussian and sampled from a (-5

3 ) spectrum, as shown in Figure 5. These perturbations should not change the large-scale flow before the predictability time (about one month). Accordingly, we require those perturbations to be of small amplitudes and restrict them to small scales. The same analysis is performed, including the EOF decomposition and the distinction between two likely scenarios. Figure 20 gathers the results. After 70 days of advection, the scenario A is hardly visible in the probability density function of the first EOF. The ensemble estimates a probability of only 39% for the reference scenario. The SQG M U ensemble estimated a probability of 47%. Moreover, the probability density function is very noisy. This suggests that the ensemble may not be converged, i.e. the empirical statistics of the ensemble will change if the ensemble size grows. The scenario A completely disappears in the joint probability density function for the two first EOF coefficients. The probability density function of the first EOF along time, before and after the symmetry breaking, exhibits very narrow branches associated with high probabilities, compared to Figure 19. It indicates that randomized initial conditions may lead to underdispersive ensemble. Furthermore, in Figure 20, trajectories of some realizations are still visible after the bifurcation. This confirms that the probability density function did not converge. Indeed, to estimate this density, we use the well-known Parzen-Rosenblatt estimator [START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF][START_REF] Parzen | On estimation of a probability density function and mode[END_REF]: each realization is associated with a kernel and the estimator is the sum of those kernels. Here, some realizations or set of few realizations are isolated and create spikes in the estimator. More realizations would be needed to have almost continuously distributed realizations. In other words, the ensemble is not converged. This drawback could be expected for at least two reasons. First, the structure of the initial noise contains little physical information, while the dimension of the state space is huge. Without phase information, covering all the possibilities requires a very large number of realizations. Furthermore, the subgrid tensor diffuses the small-scales components of the tracer where the ensemble variability is encoded. This is a known feature of ensemble forecasts: ensemble members tend to align with most unstable directions of the dynamics (Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF][START_REF] Ng | The role of model dynamics in ensemble Kalman filter performance for chaotic systems[END_REF][START_REF] Gottwald | The role of additive and multiplicative noise in filtering complex dynamical systems[END_REF]. Since small scales are stabilized by the subgrid tensor, the ensemble shrinks to span a smaller large-scale unstable subspace [START_REF] Sapsis | Attractor local dimensionality, nonlinear energy transfers and finite-time instabilities in unstable dynamical systems with applications to two-dimensional fluid flows[END_REF]. On the contrary, the stochastic model associates phase and intermittency with the noise and continuously injects it into the dynamics. The phase information or inhomogeneity as well as the non-Gaussianity come from the multiplicative structure. Even though the uncorrelated velocity is only prescribed by a spectrum, the tracer gradients have phase and dynamically constraint the regions of application of the noise. This process makes the stochastic forcing much more efficient. Hence, a smaller number of realizations are needed. In Figure 20, the convergence of the probability density function of scenario B (positive values of the first EOF coefficient) seems slightly better than the density of scenario A. Unfortunately, the bottom Figures shows that the scenario B is not the one followed by the high-resolution simulation. Let us note that the reference is deterministic. Accordingly, the reference probability density is a dirac measure. Indeed, the deterministic reference initial condition is assumed to be known and is used in all large-scale simulations. Moreover, the reference dynamic defined by the high-resolution SQG model is deterministic as the real ocean dynamics is. The bad description of the scenario A tends to confirm that the SQG model with randomized initial conditions fails to describe the bifurcation.

Another argument is the distance between the reference and each ensemble. In Figure 21 reports the error corresponding to the realization closest to the reference, i.e. with minimal error, for each ensemble. In geophysical data assimilation, large confidence is often given to observations. Thus, the maximum a posteriori estimator is almost equal to the minimal error realization. Figure 21 shows that the ensemble from the stochastic model is closer to the reference than the ensemble with randomized initializations. It suggests that the stochastic method should lead to a better maximum a posteriori estimator.

Conclusion of the long-term forecasts analysis

In this study, long-time forecasts of two different SQG models have been compared. The first one corresponds to the classical SQG equations with a random initial condition. The second one is the SALT-LU SQG, derived from a Same diagnostic based on a deterministic model with random initial conditions stochastic expression of the transport equations (see section 6). Both models are compared to a high-resolution simulation reference.

The chosen high-resolution reference is subject to a bifurcation after 40 days of advection. An infinitesimal modification of the initial condition determines the global shape of the flow two months later. Depending on the value of this initial modification, two different scenarios are isolated. For the same initial condition, the deterministic high-resolution and the low-resolution simulations do not follow the same transition. Indeed, the sub-grid tensor associated with the coarser resolution diffuses the crucial initial perturbation before the transition. This makes the deterministic forecast useless. This result questions the classical definition of predictability and associated error which only rely on initial perturbations [START_REF] Lorenz | The predictability of a flow which possesses many scales of motion[END_REF]. Note that before the symmetry breaking, the large-scale errors induced by slight modifications of the initial conditions are negligible in front of large-scale model errors. To recover, the true scenario at low resolution, a possible solution could be to randomize the initial condition. As shown, this solution would require a large number of realizations. At the opposite, the proposed stochastic model tracked both scenarios. Its efficiency is mainly explained by the continuous injection of multiplicative noise. This structure dynamically constrains both the phase and the intermittency of the noise. The model achieved to predict the likelihood and the point-wise tracer probability density in each case.

To identify and separate the scenarios, a simple threshold is sufficient regarding the high energy distributed along the first PCA axis. For more realistic flow, the number of scenarios is likely to be larger. Hence, this classification may be inadequate and more advanced clustering methods (e.g. k-means and expectation-maximization algorithms) are probably required. In this paper, this ensemble aggregation was used as a diagnosis. But, it could also improve filtering methods. Indeed, for geophysical flows, the state-space dimension being large, the ensemble size is generally too small to encode all possibilities. So, when an observation is assimilated, only few realizations of the ensemble are close enough, and all the others are considered useless. This often leads to filter degeneracy. Considering distance to likely scenarios rather than distance to realizations may help preventing this deficiency.

While the numerical analysis of this section 8 is applied to a randomized version of the SQG model with toy initial conditions, the multiplicative structure of the noise and the balance with diffusion hold for any fluid dynamics models under location uncertainty. This suggests that similar conclusions could be expected in more complex problems. Hence, this SALT-LU method opens for new ensemble forecasts methods, for both short-term forecasts and climate projections where uncertainty quantification is a main issue [START_REF] Allen | Towards objective probabilistic climate forecasting[END_REF].

Conclusion

We have reviewed existing methods to generate ensemble forecasts quantifying modes errors related to coarse resolution in computational fluid dynamics. The accuracy of this quantification being a main issue in data assimilation, the paper began by recalling the principles of that simulation-measurement coupling framework. We have categorized UQ methods based on their relations to coarse-graining. Though they have not been proposed in UQ context, we have also covered stochastic backscattering models, since they have inspired several UQ methods. Many UQ methods deal with random parameters, forcing, boundary or initial conditions. The latter type of methods has been widely misused to quantify resolution-induced error. While covariance inflation mitigates this issue, most of the data assimilation community now reckons that randomizing initial conditions underestimates resolution-induced errors. To tackle this issue, many authors introduces Gaussian and non-Gaussian noises in the dynamical equations. Most of these stochastic models are empirical, but more and more are based on energy budgets or derived from physical principles. Many of the latter assume a time scale separation in the dynamic.

After presenting a new UQ method (WavHypervis) adaptable to most deterministic subgrid dissipation, we have focused on a recent family of stochastic subgrid models: the dynamics under location uncertainty (LU) and the stochastic advection by Lie transport (SALT). These frameworks rely on a time scale separation of the velocity field, and a stochastic Navier-Stokes model. In the latter, the large-scale velocity component is transported -up to some forcings -by the small-scale and the large-scale random velocity components. Rigorously derived from stochastic calculus theory, this stochastic transport can be decomposed into a skew-symmetric multiplicative noise and an eddy-viscosity-like diffusion term. Being a transport, it naturally ensures the conservation of many physical invariants. LU and SALT models differ in the interpretations of "transport" (classical fluid dynamics VS geometric interpretation) and of "large-scale velocity". The first interpretation difference implies in particular that the dynamics under location uncertainty conserves kinetic energy whereas stochastic advection by Lie transport conserves helicity and circulation. The second interpretation difference leads to a modification of the large-scale advection in LU models. After an extended theoretical description, we have reviewed existing parameterization choices -i.e. subgrid velocity statistics modeling -for the SALT-LU framework.

Then, after briefly presenting the state-of-art metrics to assess UQ skills, new numerical results have been presented. We have compared the WavHypervis method, the dynamics under location uncertainty and the randomization of initial conditions. Test cases were short-term and long-term predictions of free decaying turbulence. The chosen dynamics was a two-dimensional meteorologic and oceanic model called the Surface Quasi-Geostrophic dynamics. For this flow, SALT and LU models mostly coincide (only the large-scale advection correction differs). Our results confirm that initial conditions randomization is not adapted to resolution-induced UQ, and that the dynamics under location uncertainty accurately spreads ensemble members along time. LU and WavHypervis short-term forecasts show very good calibrations and UQ scores. In particular, the recent LU parameterizations ADSD and SVDpseudo obtain the best UQ scores. At long term, even after a bifurcation of the reference fluid dynamics, the model under location uncertainty accurately and smoothly describes the likely scenarios.

SALT and LU dynamics are now mature frameworks. Yet, many questions remain.

First, it is still unclear whether SALT or LU is more appropriate for UQ purpose. Many numerical studies with appropriate UQ metrics and/or with data assimilation procedures would probably be necessary to try to answer this question.

Another possible research focus is the improvement of SALT-LU parameterizations. ADSD method is a promising one since it is parameter-free, fast, accurate and could in theory adapt itself to any turbulent dynamics. Nevertheless, this method is currently defined in the Fourier space. Developing a physical-domain-based ADSD implementation would be very useful. Scaling estimations of relative diffusion [START_REF] Keating | Diagnosing lateral mixing in the upper ocean with virtual tracers: Spatial and temporal resolution dependence[END_REF] or velocity structure function [START_REF] Gawȩdzki | Anomalous scaling of the passive scalar[END_REF] and Matérn covariances [START_REF] Williams | Gaussian processes for machine learning[END_REF][START_REF] Lim | Generalized whittle-matérn random field as a model of correlated fluctuations[END_REF][START_REF] Lilly | Fractional brownian motion, the Matérn process, and stochastic modeling of turbulent dispersion[END_REF][START_REF] Resseguier | Stochastic modelling and diffusion modes for POD models and small-scale flow analysis[END_REF] could probably help in this task. An adaptation of this method to complex boundary conditions would also be needed. The Gaussian unresolved velocity moments conditionally to the unresolved velocity boundary conditions may give a path toward this direction. New spatial anisotropic and heterogeneous SALT-LU parameterizations -like SVDpseudo -would also be helpful. Current works based on data-driven methods (related e.g. machine learning, Koopman operator, Girsanov-based maximum likelihood estimations) are on going.

LU and SALT long-term purpose concerns data assimilation. Premiminary studies have been published towards that direction. Using the pseudostochastic Navier-Stokes model under location uncertainty (see section 6.2.3) and a dynamics error model, [START_REF] Yang | High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4denvar[END_REF] have assimilated high-A Spectral energy flux in models under location uncertainty

In the homogeneous case, σdBt = σ * dBt and a = a 0 I d does not depend on x. Therefore, a -possibly active -tracer q is solution of: 0 = Dtq = dtq + (wdt + σdBt) • ∇q -a 0 2 ∆qdt.

For sake of simplicity, we do not explicitly show the time dependence in q and its Fourier transform, in this appendix. In Fourrier space, we get an infinite dimensional SDE (not an SPDE):

dt q(k) = - 1 (2π) d R d ŵ(k )dt + σ(k ) dBt(k ) • i(k -k )q(k -k )dk - a 0 2 k 2 2 q(k)dt. ( 127
)
To find a PDE associated to the Fourier modes of the tracer (Fokker-Planck equation), we need an SDE in R. The previous one is in C. Either we consider the closed coupled equations on the real part, (q(k)), and imaginary part, (q(k)), either we decompose trough modulus and phase. The equation on the modulus is not closed but still very instructive. From now on, the dimension of the spatial space, d, will be 2. However, we expect that the following results are still correct for d = 3. Since σdBt is divergence free, we set σ = ∇ ⊥ ψσ. We can notice that: The classical Ito formula in C, equation (127) and the above result lead to: dt|q(k)| 2 = 2 q(k)dq(k) + d < q(k), q(k) > ,

d Bt(x 1 ), Bt(x 2 ) = δ(x 1 -x 2 )dt, (128)
= -a 0 k 2 2 q(k)dt + 2 (2π

) d R d ŵ(k ) • i(k -k )q(k -k )q(k)dk =F 1 dt + 2 (2π) d R d σ(k ) dBt(k ) • i(k -k )q(k -k )q(k)dk =dF 2 (132) + (2π) d (2π) 2d R d (k -k ) T σ(k )q(k -k ) 2 dk dt. ( 133
)
dF 2 does not modify the mean spectral energy budget since E(dF 2 ) = 0, even if it modifies the law of |q(k)| 2 . Note that the mean spectral energy budget is the study of the evolution of E|q(k)| 2 whereas the spectral energy budget of the mean is the much more trivial study of |E(q(k))| 2 . We can explicit the Laplacian term using the following expression of the variance tensor, based on the homogeneity assumption:

a ij = a ij (0) = R d σi (z)σ j (z)dz = R d (∇ ⊥ ψσ(z)) i (∇ ⊥ ψσ(z)) j dz, ( 134
) = 1 (2π) d R d (∇ ⊥ ψσ(k )) i (∇ ⊥ ψσ(k )) j dk , ( 135
) = 1 (2π) d R d (k ⊥ ) i (k ⊥ ) j ψσ( k ) 2 even function if i = j, odd function if i = j dk , (136)
= 1 d 1 (2π) d R d k ψσ( k ) 2 dk =a 0 δ ij . (137)
Injecting this into the spectral energy budget yields:

dt|q(k)| 2 = F 1 dt + dF 2 - 1 (2π) d d R d k k ψσ(k )q(k) 2 dk dt + 1 (2π) d R d k • k ⊥ ψσ(k )q(k -k ) 2 dk dt, ( 138
) = F 1 dt + dF 2 + 1 (2π) d R d k k ψσ(k ) 2 sin 2 ((k, k ))|q(k -k )| 2 - 1 d |q(k)| 2 ) dk dt. ( 139
)
Taking the expectation of the previous equation highlights the terms which contribute to the mean spectral energy budget.

∂tE|q(k)| 2 = E(F 1 ) + 1 (2π) d E R d k k ψσ(k ) 2 sin 2 ((k, k ))|q(k -k )| 2 - 1 d |q(k)| 2 ) dk . (140)

B Stratonovich material derivative

Here, we express the Stratonovich material derivative:

D S t q(x, t) = q X t+ dt 2 , t + dt 2 -q X t- dt 2 , t -dt 2 | X t =x , (141)
as a function of the Ito material derivative Dtq(x, t) = (q(X t+dt , t + dt) -q(Xt, t)) | X t =x (142)

and its forcing

Dtq = F dt + H T dBt. ( 143
)
To do derive these results, we will apply the Stratonovich-Ito-notation-change formula:

r • ds = rds + 1 2 d < r, s >, (144)
which can be found in [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF], theorem 3.2.5 page 60. First, let us rewrite equation ( 143) in Lagrangian coordinates: d (q(Xt, t)) = q(X t+dt , t + dt) -q(Xt, t) = F (Xt, t)dt + H T (Xt, t)dBt.

(145)

Then, we apply the formula (144) to obtain th corresponding Stratonovich equation: After this, we can reapply formula (144) in the Eulerian grid to get: D S t q(x, t) = F -1 2 tr σ T ∇ H T (x, t)dt + H T (x, t)dBt,

d (q(Xt, t)) = q X t+ dt 2 , t + dt 2 -q X t- dt 2 , t -dt 2 , (146)
= Dtq(x, t) -1 2 tr σ T ∇ H T (x, t)dt.

(152) by identification with the balance (142).

C Effective resolution and inertial range

Let us assume the simulated evolution law is Dtq = -ν(-∆) p q dt. The deterministic subgrid model -ν(-∆) p q acts, in a finite time t, as a low-pass filter. In Fourier space, this filter is:

F ( k ) = exp -νt k 2p . ( 153
)
If the hyperviscosity ν is well chosen, we may expect that at the Shanon resolution π/∆x = κ M , only 10% of the energy is left by the filter, i.e.

F (κ M ) = 1/10.

(154)

A ratio smaller than 10% may lead to an over-damped simulation. Moreover, the precise value of this ratio does not influence much our final estimate. We may define the effective resolution as the scale κ = κm where the deterministic subgrid model influence is negligible. There, we may expect the filter to be equal to 95%, i.e.: F (κm) = 95/100.

(155)

The ratio κm/κ M can then be derived from formulas ( 153), ( 154) and ( 155).

Fig. 2

 2 Fig. 2 Two initial buoyancy (m.s -2 ) fields : (a) -A smooth field with four given vortices (cyclones in blue and anti-cyclones in red); (b) -A homogeneous Gaussian field generated from a -5/3-spectrum.

Fig. 3

 3 Fig. 3 Buoyancy (m.s -2 ) at t = 10, 15, 16, 18, 20, 30 days of advection for the usual SQG model at resolution 512 2 , based on the smooth initial field -figure 2(a).

Fig. 4

 4 Fig. 4 Buoyancy (m.s -2 ) at t = 3, 6, 9, 12, 15, 18 days of advection for the usual SQG model at resolution 512 2 , based on the random initial field -figure 2(b).

  3.2 has been implementedin SQG M U M AT LAB code. For this method -denoted "WavHypervis" -we use a wavelet basis (e k ) k∈Z . We have forecast 5 ensembles with 5 different values of the scaling factor ζ. For instance, figure6shows us that under the same low-resolution, when the scaling factor ζ increases (i.e. when the noise

Fig. 7

 7 Fig. 7 Comparisons of ensemble forecasts for different stochastic models based on the two tested inital conditions : (a) -The mean square error produced by different models under the background of 4 vortices; (b) -The corresponding difference between mean square bias and mean ensemble variance; (c), (d) -Using the empirical noise (5.3.2) with wavelet basis and various energy scaling factors, under the background of 4 vortices; (e), (f ) -Results under the spectral background. Notice that all these results are normalized by the amplitude -B 0of the referent initial buoyancy fields.

  ), (114) and Talagrand diagram. Nonetheless, these two criterions Talagrand diagram on some forecast days under the background of 4 vortices : (a) -MU Spectral; (b) -MU ADSD; (c) -WavHypervis 50%; (d) -MU SVDpseudo; (e) -PIC Spectral; (f ) -PIC SVDpseudo.

Fig. 8

 8 Fig. 8 CRPS of each model after 20 days of advection with an initial condition based on a smooth field with four given vortices

Fig. 9

 9 Fig. 9 Normalized CRPS of the model mu spectrum during 5,10 and 20 days of advection with an initial condition based on an initial random gaussian field.

Fig. 10

 10 Fig. 10 Box plot of the normalized CRPS computed on the ensemble at each spatial points of each model at different time step of advection for two different initial condition (left: Random gaussian field, right: given vortices).

Fig. 11

 11 Fig. 11 Normalized energy score for each model and spatial multivariate ensemble on all time steps and both initial condition (left: Random gaussian field, right: Four given vortices).

Fig. 12

 12 Fig. 12 Normalized energy score and variogram-0.5 score of the temporal multivariate ensemble of each models at each spatial points with an initial condition based on the four vortices (left column: normalized energy score of model, right: normalized variogram-0.5 score of each model).

  . The exact expression of the initial condition is the following: b(x, t = 0) = b 0 (x) + 0.3s w F

Fig. 13

 13 Fig.13Buoyancy (m.s -2 ) at t = 0, 30, 40, 50, 55, 60, 70 and 80 days of advection for the SQG model at resolution 512 2 . Here, a cold and a warm very small eddies are added respectively in the top and the bottom of the initial condition. These eddies are highlighted by respectively a blue square and a red square. This small difference in the initial condition does not modify the flow until the onset of the symmetry breaking, the 40 th day. Since the flow is chaotic, the small perturbation at t = 0 determines how the symmetry breaking occurs a month and a half later.

Fig. 14 Fig. 15

 1415 Fig.14Buoyancy (m.s -2 ) at t = 0, 30, 40, 50, 55, 60, 70 and 80 days of advection for the usual SQG model at resolution 512 2 . Here, two very small cold eddies are added in the top and the bottom of the initial condition. They are highlighted by two blue squares. This small difference in the initial condition does not modify the flow until the onset of the symmetry breaking, the 40 th day. Since the flow is chaotic, the small perturbation at t = 0 determines how the symmetry breaking occurs a month and a half later.

Fig. 16 Fig. 17

 1617 Fig.16Two realizations of buoyancy (m.s -2 ) after 70 days of advection for the SQG under location uncertainty at resolution 128 2 . Event though the small-scale perturbation in the initial condition is dissipated in few days, the small-scale component of the dynamics under location uncertainty randomly triggers the symmetry breaking. Therefore, some realizations follow the right transitions (left) and some do not (right).

Fig. 19

 19 Fig. 19 Probability density function of the first EOF coefficient along time (top), buoyancy (m.s -2 ) after 70 days of advection for the mean of the two likely scenarios A and B of the SQG MU model at resolution 128 2 (respectively middle left and middle right) and the reference scenarios 1 and 2 at high resolution 512 2 (respectively bottom left and bottom right). The bifurcation is clearly visible on the top plot. The two likely scenarios differ from the sign of the first EOF coefficient. They are almost associated with an equal probability of occurrence: a probability of 47% for scenario A and a probability of 53% for scenario B. The scenario A (respectively B) is similar to the scenario 1 (respectively 2).

Fig. 20

 20 Fig. 20 Energy of the EOFs (m 2 .s -4 ) (top left), probability density function of the first EOF coefficient (top middle), joint probability density function for the two first EOF coefficients (top right) after 70 days, probability density function of the first EOF coefficient along time (middle), buoyancy (m.s -2 ) after 70 days of advection for the means of the two likely scenarios A and B (bottom left and bottom right respectively). The ensemble is simulated according to the usual deterministic SQG model with random initial conditions at resolution 128 2 . The joint probability density function for the two first EOF coefficients suggests only one likely scenario. The central Figure confirms that this model fails to correctly describe the bifurcation. Indeed, the probability density function appears to be not converged in this case. The bottom Figures show that the worst resolved scenario (the scenario A) is the one similar to the true reference scenario (the scenario 1).

  y)e -ik 1 •x e +ik 2 •y dxdy dt (129) = (2π) d δ(k 1 -k 2 )dt.(130)

=

  F (Xt, t)dt -d H T (Xt, t), Bt + H T (Xt, t) • dBt, , Bt -1 2 tr σ T ∇ H T (Xt, t)dt + H T (Xt, t) • dBt,(148)since, by application of the Ito-Wentzell-formula,d H T (Xt, t) = (dH) T (Xt, t) + (σdBt • ∇) H T + ( others terms in dt). (149)Finally, rewriting everything in the Eulerian grid Xt = x gives the result: , Bt -1 2 tr σ T ∇ H T (x, t)dt + H T (x, t) • dBt.(150)

  ∆t) (see (59) above), dB t is the spatial Fourier transform of dB t , with dB t √ ∆t, a discrete scalar white noise process of unit variance in space and time. To sample the small-scale velocity, we first sample dB

t √ ∆t, to get dB t √ ∆t, and finally σ Ḃ (k) with the above equation. A MATLAB code simulating the model under location uncertainty with this parameterization for a SQG flow is available online (http://vressegu. github.io/sqgmu). Resseguier et al. (

Table 3

 3 Model performance by score, MSE: Mean squared error (pointwise), MSB: Mean squared bias (pointwise), CRPS: Continuous ranked proper score (pointwise), ES: Energy score (spatial and time).

	Model	MSE	MSB	CRPS	ES	VS
	SQG M U SP EC	+	+	+	+	+
	SQG M U ADSD	++	+	+	+	+
	SQG M U SV D	-	++	++	++	++
	W avHypervis 50	-	+	++	++	++
	SQG P IC SP EC	-	--	--	--	-
	SQG P IC SV D	--	-	+	+	+

Acknowledgments

The authors acknowledge the support of the ESA DUE GlobCurrent project, the "Laboratoires d'Excellence" CominLabs, Lebesgue and Mer through the SEACS project. The authors also acknowledge Darryl D. Holm, Francois Gay-Balmaz, Dan Crisan, Wei Pan and Igor Shevchenko for helpful discussions.

EOF spatial modes. In the present case, the buoyancy is approximated as:

≈ Ê(b)(x, t)

where N e is the size of the ensemble, N EOF N e is the number of EOF modes chosen to described the whole ensemble and the (Ψ n ) 1 n N EOF denote the EOF spatial modes. Those spatial fields are orthogonal:

with M the number of grid points. The mean energies -or variances -of EOF coefficients correspond to the eigenvalues of the two-points correlation matrix; they are ordered in decreasing order and represent the energies associated with each spatial mode. To describe the ensemble with respect to a maximal variance point of view, only the EOF coefficients c n concentrating the largest part of the buoyancy mean energy are kept. This energy, which differs from the energy of the mean, reads:

Energy of the mean

Mean energy of the n-th EOF .

(123)

In Figure 18, left part, the energy associated with the different spatial modes is displayed. At t = 30 days, the energy of the mean field, denoted as a 0-th order EOF (index+1 = 1), is much larger than the variance field, described by the other EOFs. On the contrary, at t = 70 days, after the breaking symmetry, the energy of the mean and the variance have the same order of magnitude. The variance is mainly explained by the fist EOF. Thus, at the first order, the randomness of the tracer is approximately encoded by this first EOF coefficient, c 1 . Its probability density function and the joint probability density function for the two first EOF coefficients presented in Figure 18 are unimodal at t = 30 days and bimodal at t = 70 days. Note that the tracer is clearly non-Gaussian. The symmetry breaking has created two likely scenarios in the ensemble. The scenario corresponding to negative values of the first EOF coefficient (probability of 47%) is called scenario A and the scenario corresponding to positive values (probability of 53%) is called scenario