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

Valentin Resseguier, Long Li, Gabriel Jouan, Pierre Dérian, Etienne Mémin, Bertrand Chapron

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

## To cite this version:

Valentin Resseguier, Long Li, Gabriel Jouan, Pierre Dérian, Etienne Mémin, et al. New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics. Archives of Computational Methods in Engineering, 2021, 28 (1), pp.215-261. 10.1007/s11831-020-09437-x. hal-02558016v2

HAL Id: hal-02558016

https://inria.hal.science/hal-02558016v2

Submitted on 15 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

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

Valentin Resseguier, Long Li, Gabriel Jouan, Pierre Dérian, Etienne Mémin, Bertrand Chapron

► To cite this version:

Valentin Resseguier, Long Li, Gabriel Jouan, Pierre Dérian, Etienne Mémin, et al. New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics. Archives of Computational Methods in Engineering, 2021, 28 (1), pp.215-261. 10.1007/s11831-020-09437-x. hal-02558016v2

HAL Id: hal-02558016

https://inria.hal.science/hal-02558016v2

Submitted on 15 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
New trends in ensemble forecast strategy:
uncertainty quantification for coarse-grid
computational fluid dynamics

Received: date / Accepted: date

**Abstract** Numerical simulations of industrial and geophysical fluid flows can-
not usually solve the exact Navier-Stokes equations. Accordingly, they encom-
pass 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 ap-
proaches 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 lo-
cation 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.

**Keywords** Uncertainty quantification · ensemble forecast · fluid dynamics · data assimilation · model under location uncertainty

# 1 Introduction

Understanding, modeling and tracking high Reynolds flows remain main chal-
lenges 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.
# New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics

Received: date / Accepted: date

**Abstract** Numerical simulations of industrial and geophysical fluid flows can-
not usually solve the exact Navier-Stokes equations. Accordingly, they encom-
pass 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 ap-
proaches 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 lo-
cation 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.

**Keywords** Uncertainty quantification · ensemble forecast · fluid dynamics · data assimilation · model under location uncertainty

# 1 Introduction

Understanding, modeling and tracking high Reynolds flows remain main chal-
lenges 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 (Klein et al., 2008; Gula et al., 2015), 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 backscat- terings 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 re- main as only the unresolved dissipation is in general taken into account. The unresolved subgrid components of tracer or velocity remains by definition un- known. 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 specif- ically introduced to study extreme events since there is a full literature on this subject (e.g. Franzke, 2017).

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 parameteriza- tions, 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 sim- plified model of geophysical fluid flow dynamics.
As an example, today's most accurate oceanic currents numerical simulations use a horizontal mesh resolution of about 1 km (Klein et al., 2008; Gula et al., 2015), 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 backscat- terings 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 re- main as only the unresolved dissipation is in general taken into account. The unresolved subgrid components of tracer or velocity remains by definition un- known. 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 specif- ically introduced to study extreme events since there is a full literature on this subject (e.g. Franzke, 2017).

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 parameteriza- tions, 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 sim- plified model of geophysical fluid flow dynamics.
## 2 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 Dimet and Talagrand, 1986). 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 (Doucet et al., 2001; Doucet and Johansen, 2009; Candy, 2011; Papadakis et al., 2010). 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 (Buehner, 2005; Liu et al., 2009; Yang, 2014; Bocquet and Sakov, 2014; Yang and Mémin, 2017, 2018; Sakov et al., 2017).

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 (Allen and Stainforth, 2002; Penland, 2003a,b; Berner et al., 2015). Furthermore, due to the huge state-space dimension ($\propto 10^{7-9}$), the ensemble size is usually very small in comparison ($\propto 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.
## 2 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 Dimet and Talagrand, 1986). 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 (Doucet et al., 2001; Doucet and Johansen, 2009; Candy, 2011; Papadakis et al., 2010). 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 (Buehner, 2005; Liu et al., 2009; Yang, 2014; Bocquet and Sakov, 2014; Yang and Mémin, 2017, 2018; Sakov et al., 2017).

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 (Allen and Stainforth, 2002; Penland, 2003a,b; Berner et al., 2015). Furthermore, due to the huge state-space dimension ($\propto 10^{7-9}$), the ensemble size is usually very small in comparison ($\propto 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 Maitre et al., 2002; Sapsis and Lermusiaux, 2012). 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₂ concentration conditions are difficult to specify in climate sciences. Lu-carini et al. (2014) approach this problematic with Ruelle response theory.

## 3.2 Random initial conditions

In fluid dynamics, random initial conditions have been widely used for both UQ and predictability studies (e.g. Métais and Lesieur, 1986). For a long time, operational weather forecast centers had relied on random perturbations of initial conditions to spread the ensemble forecasts (Ehrendorfer, 1997; 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 (Berner et al., 2011; Mitchell and Gottwald, 2012; Gottwald and Harlim, 2013; Franzke et al., 2015). 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 (Mitchell and Gottwald, 2012; Gottwald and Harlim, 2013). As a consequence, without any adhoc compensation, the ensemble variance is underestimated.
# 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 Maitre et al., 2002; Sapsis and Lermusiaux, 2012). 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₂ concentration conditions are difficult to specify in climate sciences. Lu-carini et al. (2014) approach this problematic with Ruelle response theory.

## 3.2 Random initial conditions

In fluid dynamics, random initial conditions have been widely used for both UQ and predictability studies (e.g. Métais and Lesieur, 1986). For a long time, operational weather forecast centers had relied on random perturbations of initial conditions to spread the ensemble forecasts (Ehrendorfer, 1997; 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 (Berner et al., 2011; Mitchell and Gottwald, 2012; Gottwald and Harlim, 2013; Franzke et al., 2015). 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 (Mitchell and Gottwald, 2012; Gottwald and Harlim, 2013). 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 (Trevisan and Uboldi, 2004; Trevisan and Palatella, 2011; Mitchell and Gottwald, 2012; Sapsis, 2013; Gottwald and Harlim, 2013). 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 (Anderson and Anderson, 1999). 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 (Tandeo et al., 2018). 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 (Raynaud et al., 2016). 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 Desroziers et al. (2005). From the *a priori* innovation covariance and the cross-covariance 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.

# 4 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 parameterizations *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
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 (Trevisan and Uboldi, 2004; Trevisan and Palatella, 2011; Mitchell and Gottwald, 2012; Sapsis, 2013; Gottwald and Harlim, 2013). 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 (Anderson and Anderson, 1999). 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 (Tandeo et al., 2018). 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 (Raynaud et al., 2016). 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 Desroziers et al. (2005). From the *a priori* innovation covariance and the cross-covariance 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.

# 4 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 parameterizations *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
<figure><img src="image_3.png" /><figcaption>**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 \in [15.5^\circ, 19.5^\circ]$.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**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 \in [15.5^\circ, 19.5^\circ]$.</figcaption></figure>
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.

## 4.1 Stochastic Lagrangian methods

First, let us mention the stochastic Lagrangian models introduced by S. B. Pope and coauthors (see Pope, 1994, 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 $\epsilon_D$. In order to compensate the induced mean kinetic energy increase, the relaxation term includes a term proportional to $\epsilon_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) (Pope, 1983). 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 (Pope and Chen, 1990) 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). Sawford (1991) 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.

Berloff and McWilliams (2002); Veneziani et al. (2004) also consider Stochastic Lagrangian models. The slightly non-linear resulting evolution models are defined on empirical grounds and Gaussian assumptions.

## 4.2 Eulerian Gaussian backscattering

The most famous Eulerian backscattering method is perhaps the Eddy-Damped Quasi Normal Markovian (EDQNM) model introduced by Orszag (1970) and Leith (1971). It closes the large-scale Navier-Stokes equations in the Fourier space by neglecting some phase correlations in the non-linear terms while
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.

## 4.1 Stochastic Lagrangian methods

First, let us mention the stochastic Lagrangian models introduced by S. B. Pope and coauthors (see Pope, 1994, 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 $\epsilon_D$. In order to compensate the induced mean kinetic energy increase, the relaxation term includes a term proportional to $\epsilon_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) (Pope, 1983). 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 (Pope and Chen, 1990) 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). Sawford (1991) 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.

Berloff and McWilliams (2002); Veneziani et al. (2004) also consider Stochastic Lagrangian models. The slightly non-linear resulting evolution models are defined on empirical grounds and Gaussian assumptions.

## 4.2 Eulerian Gaussian backscattering

The most famous Eulerian backscattering method is perhaps the Eddy-Damped Quasi Normal Markovian (EDQNM) model introduced by Orszag (1970) and Leith (1971). 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. Chasnov (1991) uses this framework to set up
a forced-dissipative Navier-Stokes model where the Eddy Viscosity is scale-
dependant 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 (Farrell and Ioannou, 2014). 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 co-
variance in order to reveal and characterize various turbulent phenomenon.
Stochastic superparameterization (SSP) proposes a similar model (Grooms
and Majda, 2014). 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, Berloff (2005) 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.

4.3 Eulerian non-Gaussian backscattering

Phase information can also be encoded by multiplicative noises. Besides, mul-
tiplicative noises are the most common non-Gaussian forcing in CFD. Leith
(1990) multiplies a white Gaussian noise by a function of the resolved local
strain rate. Schumann (1995) uses a quadratic function of a homogeneous
Gaussian noise. Brankart (2013) adds at each time step a multiplicative noise
to the active tracers (salinity and temperature) before computing the corre-
sponding density. After this, the density is averaged over realizations. The non-
linearity of the state-equation makes this transitional variability non-negligible
and improves the simulation results.

In contrast, Mana and Zanna (2014) 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 Schumann (1995), 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 turbu-
lent dissipation, but the justification of this link is generally either unclear or
keeping the energy constant. Chasnov (1991) uses this framework to set up
a forced-dissipative Navier-Stokes model where the Eddy Viscosity is scale-
dependant 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 (Farrell and Ioannou, 2014). 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 co-
variance in order to reveal and characterize various turbulent phenomenon.
Stochastic superparameterization (SSP) proposes a similar model (Grooms
and Majda, 2014). 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, Berloff (2005) 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.

4.3 Eulerian non-Gaussian backscattering

Phase information can also be encoded by multiplicative noises. Besides, mul-
tiplicative noises are the most common non-Gaussian forcing in CFD. Leith
(1990) multiplies a white Gaussian noise by a function of the resolved local
strain rate. Schumann (1995) uses a quadratic function of a homogeneous
Gaussian noise. Brankart (2013) adds at each time step a multiplicative noise
to the active tracers (salinity and temperature) before computing the corre-
sponding density. After this, the density is averaged over realizations. The non-
linearity of the state-equation makes this transitional variability non-negligible
and improves the simulation results.

In contrast, Mana and Zanna (2014) 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 Schumann (1995), 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 turbu-
lent dissipation, but the justification of this link is generally either unclear or
arbitrary.

# 5 UQ for errors induced by a coarsening process

## 5.1 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 (Penland and Matrosova, 1994; Penland and Sardeshmukh, 1995). Keating et al. (2012) 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 Keating et al. (2012), 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) (Gottwald and Harlim, 2013).

In geophysical fluid dynamics data assimilation, dynamics are often non-linear. Nevertheless, Gaussian additive random forcings are still widely used. Tandeo et al. (2018) 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 (Pannekoucke and Massart, 2008). 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 (Pannekoucke et al., 2014). Then, a classical covariance model can be fitted.

## 5.2 Empirical non-Gaussian noises

Among other empirical stochastic models, the Stochastic Perturbed Physics Tendency scheme (SPPT) introduces a correlated multiplicative noise (Buizza et al., 1999). All parameterized tendencies are multiplied by the same uniform variable in [0.5, 1.5]. The same random variable is used in a fixed space-time window. This window sets the noise spatio-temporal correlations. Shutts (2005); Berner et al. (2009, 2011) 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 Schumann (1995) and Brankart (2013), 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
arbitrary.

# 5 UQ for errors induced by a coarsening process

## 5.1 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 (Penland and Matrosova, 1994; Penland and Sardeshmukh, 1995). Keating et al. (2012) 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 Keating et al. (2012), 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) (Gottwald and Harlim, 2013).

In geophysical fluid dynamics data assimilation, dynamics are often non-linear. Nevertheless, Gaussian additive random forcings are still widely used. Tandeo et al. (2018) 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 (Pannekoucke and Massart, 2008). 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 (Pannekoucke et al., 2014). Then, a classical covariance model can be fitted.

## 5.2 Empirical non-Gaussian noises

Among other empirical stochastic models, the Stochastic Perturbed Physics Tendency scheme (SPPT) introduces a correlated multiplicative noise (Buizza et al., 1999). All parameterized tendencies are multiplied by the same uniform variable in [0.5, 1.5]. The same random variable is used in a fixed space-time window. This window sets the noise spatio-temporal correlations. Shutts (2005); Berner et al. (2009, 2011) 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 Schumann (1995) and Brankart (2013), 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 (Shutts, 2015). 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 (Franzke et al., 2015; Buizza, 2016b; Pegion et al., 2016; Separovic et al., 2016; McCabe et al., 2016; Buizza, 2016a). Nevertheless, many drawbacks of these methods have been reported. Above all, conservation laws (e.g. energy and mass conservations) are violated (Reynolds et al., 2016; Leutbechner et al., 2016). 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 (Reynolds et al., 2016).

## 5.3 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.

### 5.3.1 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) (Grooms and Majda, 2014), 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. Majda (2015) 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, Majda (2015) does not specify the random forcing. This is why the latter does not a priori balance the turbulent diffusion. Farrell and Ioannou (2014) also study the energy mean of stochastic fluid dynamics systems especially under quasi-linear approximations and with an additive Gaussian force.
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 (Shutts, 2015). 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 (Franzke et al., 2015; Buizza, 2016b; Pegion et al., 2016; Separovic et al., 2016; McCabe et al., 2016; Buizza, 2016a). Nevertheless, many drawbacks of these methods have been reported. Above all, conservation laws (e.g. energy and mass conservations) are violated (Reynolds et al., 2016; Leutbechner et al., 2016). 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 (Reynolds et al., 2016).

## 5.3 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.

### 5.3.1 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) (Grooms and Majda, 2014), 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. Majda (2015) 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, Majda (2015) does not specify the random forcing. This is why the latter does not a priori balance the turbulent diffusion. Farrell and Ioannou (2014) also study the energy mean of stochastic fluid dynamics systems especially under quasi-linear approximations and with an additive Gaussian force.
ing. Gugole and Franzke (219) enforce the energy conservation by projecting their random forcing into the space of constant energy.

As cited previously, SKEBS parameterization (Shutts, 2005; Berner et al., 2009, 2011) backscatters a given portion of the global dissipated energy using a flow-dependent correlated noise. Jansen and Held (2014) work with a similar idea. Yet, they consider a hyper-viscous diffusion and a noise white in space and time. Dwivedi et al. (2019) also consider hyper-viscosity and neglect spatial correlations but they keep time correlations.

### 5.3.2 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 $\mathcal{L}$ and a stochastic forcing $\eta$ – which is a centred process uncorrelated in time:

$$ \frac{Dq}{Dt} = \mathcal{L}[q] + \eta, \qquad (1) $$

where $D/Dt = \partial_t + v \cdot \nabla$ stands for the material derivative operator.

Through Itô calculus, a conservation of the energy mean – $\mathbb{E} \int_{\Omega} \frac{1}{2} q^2$ – would imply:

$$ 0 = \mathbb{E} \int_{\Omega} \frac{d}{dt} (\frac{1}{2} q^2) = \mathbb{E} \int_{\Omega} (q \frac{d}{dt} q + \frac{1}{2} \frac{d}{dt} \langle q, q \rangle) = \mathbb{E} \int_{\Omega} (q \mathcal{L}[q] + \frac{1}{2} (\mathcal{H}[q])^2) dt, \qquad (2) $$

where

$$ (\mathcal{H}[q](x, t))^2 \triangleq \mathbb{E}(\eta^2(x, t) | q(\bullet, t)), \qquad (3) $$

is the stochastic forcing variance conditioned on field $q$ at time $t$ and $\langle q, q \rangle$ denotes the quadratic variation term. In order to maintain a desired amount of energy dissipation, we introduce a scaling factor $\zeta > 0$ in the above balance. Specifically, we want to impose:

$$ \int_{\Omega} q \left( -\frac{2\zeta}{dt} \mathcal{L}[q] \right) = \int_{\Omega} (\mathcal{H}[q])^2. \qquad (4) $$

Now, for a given dissipation operator $\mathcal{L}$, we will construct a noise $\eta$ which meets the desired balance (4). If $\mathcal{L}$ is a negative auto-adjoint operator – which is generally the case – we can define an operator $H$ which is such that

$$ -\frac{2\zeta}{dt} \mathcal{L} = HH^*, \qquad (5) $$
ing. Gugole and Franzke (219) enforce the energy conservation by projecting their random forcing into the space of constant energy.

As cited previously, SKEBS parameterization (Shutts, 2005; Berner et al., 2009, 2011) backscatters a given portion of the global dissipated energy using a flow-dependent correlated noise. Jansen and Held (2014) work with a similar idea. Yet, they consider a hyper-viscous diffusion and a noise white in space and time. Dwivedi et al. (2019) also consider hyper-viscosity and neglect spatial correlations but they keep time correlations.

### 5.3.2 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 $\mathcal{L}$ and a stochastic forcing $\eta$ – which is a centred process uncorrelated in time:

$$ \frac{Dq}{Dt} = \mathcal{L}[q] + \eta, \qquad (1) $$

where $D/Dt = \partial_t + v \cdot \nabla$ stands for the material derivative operator.

Through Itô calculus, a conservation of the energy mean – $\mathbb{E} \int_{\Omega} \frac{1}{2} q^2$ – would imply:

$$ 0 = \mathbb{E} \int_{\Omega} \frac{d}{dt} (\frac{1}{2} q^2) = \mathbb{E} \int_{\Omega} (q \frac{d}{dt} q + \frac{1}{2} \frac{d}{dt} \langle q, q \rangle) = \mathbb{E} \int_{\Omega} (q \mathcal{L}[q] + \frac{1}{2} (\mathcal{H}[q])^2) dt, \qquad (2) $$

where

$$ (\mathcal{H}[q](\mathbf{x}, t))^2 \triangleq \mathbb{E}(\eta^2(\mathbf{x}, t) | q(\bullet, t)), \qquad (3) $$

is the stochastic forcing variance conditioned on field $q$ at time $t$ and $\langle q, q \rangle$ denotes the quadratic variation term. In order to maintain a desired amount of energy dissipation, we introduce a scaling factor $\zeta > 0$ in the above balance. Specifically, we want to impose:

$$ \int_{\Omega} q \left(-\frac{2\zeta}{dt} \mathcal{L}[q]\right) = \int_{\Omega} (\mathcal{H}[q])^2. \qquad (4) $$

Now, for a given dissipation operator $\mathcal{L}$, we will construct a noise $\eta$ which meets the desired balance (4). If $\mathcal{L}$ is a negative auto-adjoint operator – which is generally the case – we can define an operator $H$ which is such that

$$ -\frac{2\zeta}{dt} \mathcal{L} = HH^*, \qquad (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] \triangleq \left\{
\begin{array}{ll}
\tilde{\zeta} \alpha \nabla q & \text{if } \mathcal{L}[q] = \nabla \cdot (\alpha \alpha^T \nabla q), \\
\tilde{\zeta} \nabla \cdot (\alpha \nabla q) & \text{if } \mathcal{L}[q] = -\nabla \cdot (\alpha \nabla (\nabla \cdot (\alpha \nabla q))), \\
\tilde{\zeta} \alpha \nabla \Delta^p q & \text{if } \mathcal{L}[q] = \alpha^2 \Delta^{2p+1} q, \\
\tilde{\zeta} \alpha \Delta^p q & \text{if } \mathcal{L}[q] = -\alpha^2 \Delta^{2p} q,
\end{array}
\right.
\quad (6)
$$

where $\tilde{\zeta} = \sqrt{\frac{2c}{dt}}$, $p$ is a positive integer, $\alpha$ a constant and $\alpha$ a – possibly spatially-varying – matrix.

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

$$
\int_{\Omega} (\mathcal{H}[q])^2 = \int_{\Omega} \|\mathbf{H}[q]\|^2 = \int_{\Omega} q \mathbf{H}^* \mathbf{H}[q] = \int_{\Omega} q \left( -\frac{2c}{dt} \mathcal{L}[q] \right). \quad (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²(Ω) (e.g. a wavelet basis)
denoted {e_k}_{k≥0}:

$$
\eta(\mathbf{x}, t) = \sum_{k \ge 0} \lambda_k(t) \xi_k(t) e_k(\mathbf{x}), \quad (8)
$$

where the $\xi_k$ are independent Gaussian white noises of unit variance and:

$$
\lambda_k = \int_{\Omega} \|H[q]\| e_k. \quad (9)
$$

Then, by definition of $\mathcal{H}$ and by the Parseval theorem, the random forcing (conditional) variance integral is given by:

$$
\int_{\Omega} (\mathcal{H}[q])^2 = \int_{\Omega} \mathbb{E}(\eta^2(\mathbf{x}, t) | q(\bullet, t)) = \sum_{k \ge 0} |\lambda_k|^2 = \int_{\Omega} \|\mathbf{H}[q]\|^2. \quad (10)
$$

Therefore, equation (7) is valid and the stochastic forcing $\eta$ defined by (8) and (9) meets our target energy balance (4).

Therefore, given a numerical dissipation $\mathcal{L}$ together with a chosen orthonor-
mal 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:

$$
(\mathcal{H}[q])^2 = \sum_{k \ge 0} |\lambda_k|^2 |e_k|^2 \neq \left| \sum_{k \ge 0} \lambda_k e_k \right|^2 = \|\mathbf{H}[q]\|^2, \quad (11)
$$
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] \triangleq \begin{cases} \tilde{\zeta} \alpha \nabla q & \text{if } \mathcal{L}[q] = \nabla \cdot (\alpha \alpha^T \nabla q), \\ \tilde{\zeta} \nabla \cdot (\alpha \nabla q) & \text{if } \mathcal{L}[q] = -\nabla \cdot (\alpha \nabla (\nabla \cdot (\alpha \nabla q))), \\ \tilde{\zeta} \alpha \nabla \Delta^p q & \text{if } \mathcal{L}[q] = \alpha^2 \Delta^{2p+1} q, \\ \tilde{\zeta} \alpha \Delta^p q & \text{if } \mathcal{L}[q] = -\alpha^2 \Delta^{2p} q, \end{cases} \quad (6) $$

where $\tilde{\zeta} = \sqrt{\frac{2\zeta}{dt}}$, $p$ is a positive integer, $\alpha$ a constant and $\alpha$ a – possibly spatially-varying – matrix.

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

$$ \int_{\Omega} (\mathcal{H}[q])^2 = \int_{\Omega} \|\mathcal{H}[q]\|^2 = \int_{\Omega} q H^* H[q] = \int_{\Omega} q \left(-\frac{2\zeta}{dt} \mathcal{L}[q]\right). \quad (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(\Omega)$ (e.g. a wavelet basis) denoted $\{e_k\}_{k \ge 0}$:

$$ \eta(\mathbf{x}, t) = \sum_{k \ge 0} \lambda_k(t) \xi_k(t) e_k(\mathbf{x}), \quad (8) $$

where the $\xi_k$ are independent Gaussian white noises of unit variance and:

$$ \lambda_k = \int_{\Omega} \|\mathcal{H}[q]\| e_k. \quad (9) $$

Then, by definition of $\mathcal{H}$ and by the Parseval theorem, the random forcing (conditional) variance integral is given by:

$$ \int_{\Omega} (\mathcal{H}[q])^2 = \int_{\Omega} \mathbb{E}(\eta^2(\mathbf{x}, t) | q(\bullet, t)) = \sum_{k \ge 0} |\lambda_k|^2 = \int_{\Omega} \|\mathcal{H}[q]\|^2. \quad (10) $$

Therefore, equation (7) is valid and the stochastic forcing $\eta$ defined by (8) and (9) meets our target energy balance (4).

Therefore, given a numerical dissipation $\mathcal{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:

$$ (\mathcal{H}[q])^2 = \sum_{k \ge 0} |\lambda_k|^2 |e_k|^2 \neq \left| \sum_{k \ge 0} \lambda_k e_k \right|^2 = \|\mathcal{H}[q]\|^2, \quad (11) $$
and moreover,

$$ \|H[q]\|^2 \neq qH^*H[q] = q\left(-\frac{2\zeta}{dt}\mathcal{L}[q]\right). \qquad (12) $$

The choice of a given basis {$e_k$}$_{k \ge 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| = \text{cst}$). In this case, the inequality (11) would be $(\mathcal{H}[q])^2 = \frac{1}{\mu(\Omega)} \int_{\Omega} \|H[q]\|^2 \neq \|H[q]\|^2(\boldsymbol{x}, t)$, where $\mu(\Omega)$ is the measure of the spatial domain $\Omega$. In contrast a basis of regularized Dirac functions $e_k(\boldsymbol{x}) = \delta(\boldsymbol{x} - \boldsymbol{x}_k)$ – defined on the grid points $\boldsymbol{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.

## 5.4 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 (Berner et al., 2015; Craig et al., 2016). Yet as Navier-Stokes equations are deterministic, this path is not straightforward.

### 5.4.1 Time-scale separation

#### Averaging and homogenization

To derive large-scale fluid dynamical model, the time-scale separation assumption is convenient. In the seventies, Hasselmann (1976) 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 (Kurtz, 1973; Papanicolaou and Kohler, 1974; Givon et al., 2004; Mitchell and Gottwald, 2012; Gottwald and Melbourne, 2013; Gottwald and Harlim, 2013; Franzke et al., 2015; Gottwald et al., 2015). 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
and moreover,

$$ \|H[q]\|^2 \neq qH^*H[q] = q\left(-\frac{2\zeta}{dt}\mathcal{L}[q]\right). \qquad (12) $$

The choice of a given basis {$e_k$}$_{k \ge 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| = \text{cst}$). In this case, the inequality (11) would be $(\mathcal{H}[q])^2 = \frac{1}{\mu(\Omega)} \int_{\Omega} \|H[q]\|^2 \neq \|H[q]\|^2(x, t)$, where $\mu(\Omega)$ is the measure of the spatial domain $\Omega$. In contrast a basis of regularized Dirac functions $e_k(x) = \delta(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.

## 5.4 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 (Berner et al., 2015; Craig et al., 2016). Yet as Navier-Stokes equations are deterministic, this path is not straightforward.

### 5.4.1 Time-scale separation

#### Averaging and homogenization

To derive large-scale fluid dynamical model, the time-scale separation assumption is convenient. In the seventies, Hasselmann (1976) 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 (Kurtz, 1973; Papanicolaou and Kohler, 1974; Givon et al., 2004; Mitchell and Gottwald, 2012; Gottwald and Melbourne, 2013; Gottwald and Harlim, 2013; Franzke et al., 2015; Gottwald et al., 2015). 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 (Gottwald and Melbourne, 2013; Gottwald et al., 2015).

A successful application of the homogenization theory in geophysics is the MTV scheme (Majda et al., 1999, 2001; Franzke et al., 2005; Majda et al., 2008). 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. Peavoy et al. (2015) 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 (Givon et al., 2004; Gottwald and Harlim, 2013). This scheme provided good UQ skills especially because of the multiplicative noise appearing in the limit dynamics. Chekroun et al. (2017) 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 (Frank and Gottwald, 2013; Jain et al., 2014). 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 homogenization procedures, Edgeworth expansions goes to higher order in time scale ratio (Wouters and Gottwald, 2018).

### Velocity time scale separation and skew-symmetric noises

The Kraichnan model is an idealization of passive tracer turbulence (Kraichnan, 1968, 1994; Gawędzki and Kupiainen, 1995). A tracer is forced and advected by two independent spatially-homogeneous time-uncorrelated Gaussian
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 (Gottwald and Melbourne, 2013; Gottwald et al., 2015).

A successful application of the homogenization theory in geophysics is the MTV scheme (Majda et al., 1999, 2001; Franzke et al., 2005; Majda et al., 2008). 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. Peavoy et al. (2015) 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 (Givon et al., 2004; Gottwald and Harlim, 2013). This scheme provided good UQ skills especially because of the multiplicative noise appearing in the limit dynamics. Chekroun et al. (2017) 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 (Frank and Gottwald, 2013; Jain et al., 2014). 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 homogenization procedures, Edgeworth expansions goes to higher order in time scale ratio (Wouters and Gottwald, 2018).

### Velocity time scale separation and skew-symmetric noises

The Kraichnan model is an idealization of passive tracer turbulence (Kraichnan, 1968, 1994; Gawędzki and Kupiainen, 1995). 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 (Gawędzki and Kupiainen, 1995; Klyatskin et al., 1996; Majda et al., 1999; Falkovich et al., 2001). 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. Klyatskin et al., 1996).

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 (Kunita, 1997). Brzezniak et al. (1991) first introduces the idea. Then, Mikulevicius and Rozovskii (2004) and Flandoli (2011) modify the formula. Their works have focused on pure mathematical aims: existence and uniqueness of SPDE solutions. Neves and Olivera (2015) also studied the well-posedness 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; Chapron et al., 2018; Cai et al., 2018; Resseguier et al., 2019) and Stochastic Advection by Lie Transport (SALT) (Holm, 2015; Crisan et al., 2017; Gay-Balmaz and Holm, 2018; Cotter et al., 2018c,a,b; Resseguier et al., 2019) are new types of stochastic subgrid tensors based on that stochastic transport. These schemes will be deeply detailed in section 6. Sardeshmukh and Sura (2009); Sardeshmukh et al. (2015) 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 (Cotter et al., 2017). 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.

### 5.4.2 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 (Givon et al., 2004; Gottwald et al., 2015). Indeed, this explicit expression of the large-scale observables of interest involves a memory term. Accordingly, some authors have proposed non-Markovian stochastic
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 (Gawędzki and Kupiainen, 1995; Klyatskin et al., 1996; Majda et al., 1999; Falkovich et al., 2001). 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. Klyatskin et al., 1996).

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 (Kunita, 1997). Brzezniak et al. (1991) first introduces the idea. Then, Mikulevicius and Rozovskii (2004) and Flandoli (2011) modify the formula. Their works have focused on pure mathematical aims: existence and uniqueness of SPDE solutions. Neves and Olivera (2015) also studied the well-posedness 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; Chapron et al., 2018; Cai et al., 2018; Resseguier et al., 2019) and Stochastic Advection by Lie Transport (SALT) (Holm, 2015; Crisan et al., 2017; Gay-Balmaz and Holm, 2018; Cotter et al., 2018c,a,b; Resseguier et al., 2019) are new types of stochastic subgrid tensors based on that stochastic transport. These schemes will be deeply detailed in section 6. Sardeshmukh and Sura (2009); Sardeshmukh et al. (2015) 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 (Cotter et al., 2017). 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.

### 5.4.2 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 (Givon et al., 2004; Gottwald et al., 2015). 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 Chekroun et al., 2011; Kondrashov et al., 2015; Kondrashov and Berloff, 2015; Lu et al., 2017). 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.

### 5.4.3 Other approaches based on statistical physics

Using Ruelle response theory, Wouters and Lucarini (2012) 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), Frederiksen (1999) 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.

Plant and Craig (2008) have proposed a physically-based stochastic subgrid parameterization to model mass flux induced by clouds updrafts and down-drafts. Nevertheless, its application is a priori limited to weather forecasts.

## 6 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. Holm, 2015). 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.

### 6.1 Stochastic transport

#### 6.1.1 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
parameterizations (e.g Chekroun et al., 2011; Kondrashov et al., 2015; Kondrashov and Berloff, 2015; Lu et al., 2017). 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.

### 5.4.3 Other approaches based on statistical physics

Using Ruelle response theory, Wouters and Lucarini (2012) 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), Frederiksen (1999) 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.

Plant and Craig (2008) 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.

## 6 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. Holm, 2015). 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.

### 6.1 Stochastic transport

#### 6.1.1 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 intermit-
tent energy bursts, possibly associated with infinitesimal characteristic time-
scales, 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 $\mathbf{w}$, repre-
sents a large-scale, possibly random, component continuous in time. The un-
resolved contribution, expressed as $\boldsymbol{\sigma}\dot{\mathbf{B}}$, 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(\mathbf{X}_{t+\Delta t}, t + \Delta t) = q(\mathbf{X}_t, t) \quad (13)
$$

reads in an informal stochastic way as

$$
\partial_t q + \underbrace{\mathbf{w}^* \cdot \nabla q}_{\text{Corrected advection}} = \underbrace{\nabla \cdot (\frac{1}{2} a \nabla q)}_{\text{Diffusion}} - \underbrace{\sigma \dot{\mathbf{B}} \cdot \nabla q}_{\text{Random forcing}} \quad (14)
$$

with a drift velocity corrected as

$$
\mathbf{w}^* = \mathbf{w} - \frac{1}{2} (\nabla \cdot \mathbf{a})^T + \sigma (\nabla \cdot \mathbf{\sigma})^T. \quad (15)
$$

Hence, the deterministic (for a fixed realization of velocity field **w**) evolution
equation is replaced by a stochastic equation with respect to **σ**̇**B**. 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 in-
deed essential to take into account essential physical effects due the small-scale
velocity inhomogeneity (Chandramouli et al., 2019; Pinier et al., 2019). The
random forcing term in (14) relates to the advection by the unresolved velocity
**σ**̇**B** = **σ**\*d**B**/dt. This term continuously backscatters random energy to the sys-
tem. 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 mod-
els introduced in classical large-scale representations (Gent and Mcwilliams,
1990; Germano et al., 1991; Lilly, 1966; Smagorinsky, 1963), in analogy with
the molecular diffusion mechanism (Boussinesq, 1877). In particular, for a ho-
mogeneous and isotropic small-scale velocity the turbulent diffusion is also
homogeneous and isotropic:

$$
a = a_0 \mathbb{I}_d \text{ and } \nabla \cdot \left( \frac{1}{2} a \nabla q \right) = \frac{1}{2} a_0 \Delta q. \quad (16)
$$

The inhomogeneous and anisotropic diffusion coefficient matrix, **a**, is then
defined by the one-point one-time covariance of the unresolved displacement
small-scale fluid flow velocities can be characterized by local and intermit-
tent energy bursts, possibly associated with infinitesimal characteristic time-
scales, 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 $\mathbf{w}$, repre-
sents a large-scale, possibly random, component continuous in time. The un-
resolved contribution, expressed as $\boldsymbol{\sigma}\dot{\mathbf{B}}$, 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(\mathbf{X}_{t+\Delta t}, t+\Delta t) = q(\mathbf{X}_t, t) \qquad (13)
$$

reads in an informal stochastic way as

$$
\partial_t q + \underbrace{\mathbf{w}^* \cdot \nabla q}_{\text{Corrected advection}} = \underbrace{\nabla \cdot (\frac{1}{2} a \nabla q)}_{\text{Diffusion}} - \underbrace{\sigma \dot{\mathbf{B}} \cdot \nabla q}_{\text{Random forcing}}, \quad (14)
$$

with a drift velocity corrected as

$$
\mathbf{w}^* = \mathbf{w} - \frac{1}{2} (\nabla \cdot \mathbf{a})^T + \sigma (\nabla \cdot \mathbf{\sigma})^T. \quad (15)
$$

Hence, the deterministic (for a fixed realization of velocity field $\mathbf{w}$) evolution
equation is replaced by a stochastic equation with respect to $\boldsymbol{\sigma}\dot{\mathbf{B}}$. As a result,
the conserved quantity is now advected by an “effective” velocity, $\mathbf{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 in-
deed essential to take into account essential physical effects due the small-scale
velocity inhomogeneity (Chandramouli et al., 2019; Pinier et al., 2019). The
random forcing term in (14) relates to the advection by the unresolved velocity
$\boldsymbol{\sigma}\dot{\mathbf{B}} = \boldsymbol{\sigma}\frac{dB_t}{dt}$. This term continuously backscatters random energy to the sys-
tem. 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 mod-
els introduced in classical large-scale representations (Gent and Mcwilliams,
1990; Germano et al., 1991; Lilly, 1966; Smagorinsky, 1963), in analogy with
the molecular diffusion mechanism (Boussinesq, 1877). In particular, for a ho-
mogeneous and isotropic small-scale velocity the turbulent diffusion is also
homogeneous and isotropic:

$$
\mathbf{a} = a_0 \mathbb{I}_d \text{ and } \nabla \cdot \left( \frac{1}{2} \mathbf{a} \nabla q \right) = \frac{1}{2} a_0 \Delta q. \qquad (16)
$$

The inhomogeneous and anisotropic diffusion coefficient matrix, $\mathbf{a}$, is then
defined by the one-point one-time covariance of the unresolved displacement
per unit of time:

$$ \boldsymbol{a} = \frac{\mathbb{E}\{\boldsymbol{\sigma} \mathrm{d}\boldsymbol{B}_t (\boldsymbol{\sigma} \mathrm{d}\boldsymbol{B}_t)^T\}}{\mathrm{d}t}. \qquad (17) $$

This ensures an exact energy balance between the amount of diffusion and the random forcing.

### 6.1.2 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 $\boldsymbol{X}_t$ writes:

$$ \mathrm{d}\boldsymbol{X}_t = \boldsymbol{w}(\boldsymbol{X}_t, t)\mathrm{d}t + \boldsymbol{\sigma}(\boldsymbol{X}_t, t)\mathrm{d}\boldsymbol{B}_t. \qquad (18) $$

In this equation, the second term explicitly figures the flow location uncertainty. Formally, this random field is defined over the fluid domain, $\Omega \subset \mathbb{R}^d$, from a $d$-dimensional Brownian function $\boldsymbol{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 Prato and Zabczyk, 1992; Prévôt and Röckner, 2007, 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 $\boldsymbol{\sigma}\dot{\boldsymbol{B}} = \boldsymbol{\sigma}\frac{\mathrm{d}\boldsymbol{B}_t}{\mathrm{d}t}$, and is a white noise distribution. The spatial correlations of the flow uncertainty are specified through the diffusion operator $\boldsymbol{\sigma}(\cdot, t)$, defined for any vectorial function, $\boldsymbol{f}$, through the matrix kernel $\tilde{\boldsymbol{\sigma}}(\cdot, \cdot, t)$:

$$ \boldsymbol{\sigma}(\boldsymbol{x}, t)\boldsymbol{f} \triangleq \int_{\Omega} \tilde{\boldsymbol{\sigma}}(\boldsymbol{x}, z, t)\boldsymbol{f}(z, t)\mathrm{d}z. \qquad (19) $$

This quantity is assumed to have a finite norm. More precisely, the operator $\boldsymbol{\sigma}$ 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, $\boldsymbol{\sigma}(\boldsymbol{x}, t)\mathrm{d}\boldsymbol{B}_t$, is a centered vectorial Gaussian function, correlated in space and uncorrelated in time with covariance tensor:

$$ \mathrm{Cov}(\boldsymbol{x}, \boldsymbol{y}, t, t') \triangleq \mathbb{E}\{(\boldsymbol{\sigma}(\boldsymbol{x}, t)\mathrm{d}\boldsymbol{B}_t)(\boldsymbol{\sigma}(\boldsymbol{y}, t')\mathrm{d}\boldsymbol{B}_{t'} )^T\}, \qquad (20) $$

$$ = \int_{\Omega} \tilde{\boldsymbol{\sigma}}(\boldsymbol{x}, \boldsymbol{z}, t) \tilde{\boldsymbol{\sigma}}^{T}(\boldsymbol{y}, \boldsymbol{z}, t) \mathrm{d}\boldsymbol{z} \delta(t - t') \mathrm{d}t. \qquad (21) $$

Hereafter, the diagonal of the covariance tensor, $\boldsymbol{a}$, – also denoted as $\boldsymbol{\sigma}\boldsymbol{\sigma}^T$ – will be referred to as the variance tensor:

$$ \boldsymbol{a}(\boldsymbol{x}, t)\delta(t - t')\mathrm{d}t = \mathrm{Cov}(\boldsymbol{x}, \boldsymbol{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.
per unit of time:

$$ \boldsymbol{a} = \frac{\mathbb{E}\{\boldsymbol{\sigma} \mathrm{d}\boldsymbol{B}_t (\boldsymbol{\sigma} \mathrm{d}\boldsymbol{B}_t)^T\}}{\mathrm{d}t}. \qquad (17) $$

This ensures an exact energy balance between the amount of diffusion and the random forcing.

### 6.1.2 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 $\boldsymbol{X}_t$ writes:

$$ \mathrm{d}\boldsymbol{X}_t = \boldsymbol{w}(\boldsymbol{X}_t, t)\mathrm{d}t + \boldsymbol{\sigma}(\boldsymbol{X}_t, t)\mathrm{d}\boldsymbol{B}_t. \qquad (18) $$

In this equation, the second term explicitly figures the flow location uncertainty. Formally, this random field is defined over the fluid domain, $\Omega \subset \mathbb{R}^d$, from a $d$-dimensional Brownian function $\boldsymbol{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 Prato and Zabczyk, 1992; Prévôt and Röckner, 2007, 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 $\boldsymbol{\sigma}\dot{\boldsymbol{B}} = \boldsymbol{\sigma}\frac{\mathrm{d}\boldsymbol{B}_t}{\mathrm{d}t}$, and is a white noise distribution. The spatial correlations of the flow uncertainty are specified through the diffusion operator $\boldsymbol{\sigma}(., t)$, defined for any vectorial function, $\boldsymbol{f}$, through the matrix kernel $\tilde{\boldsymbol{\sigma}}(., ., t)$:

$$ \boldsymbol{\sigma}(\boldsymbol{x}, t)\boldsymbol{f} \triangleq \int_{\Omega} \tilde{\boldsymbol{\sigma}}(\boldsymbol{x}, z, t)\boldsymbol{f}(z, t)\mathrm{d}z. \qquad (19) $$

This quantity is assumed to have a finite norm. More precisely, the operator $\boldsymbol{\sigma}$ 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, $\boldsymbol{\sigma}(\boldsymbol{x}, t)\mathrm{d}\boldsymbol{B}_t$, is a centered vectorial Gaussian function, correlated in space and uncorrelated in time with covariance tensor:

$$ \mathrm{Cov}(\boldsymbol{x}, \boldsymbol{y}, t, t') \triangleq \mathbb{E}\{(\boldsymbol{\sigma}(\boldsymbol{x}, t)\mathrm{d}\boldsymbol{B}_t)(\boldsymbol{\sigma}(\boldsymbol{y}, t')\mathrm{d}\boldsymbol{B}_{t'} )^T\}, \qquad (20) $$

$$ = \int_{\Omega} \tilde{\boldsymbol{\sigma}}(\boldsymbol{x}, \boldsymbol{z}, t) \tilde{\boldsymbol{\sigma}}^{T}(\boldsymbol{y}, \boldsymbol{z}, t) \mathrm{d}\boldsymbol{z} \delta(t - t') \mathrm{d}t. \qquad (21) $$

Hereafter, the diagonal of the covariance tensor, $\boldsymbol{a}$, – also denoted as $\boldsymbol{\sigma}\boldsymbol{\sigma}^T$ – will be referred to as the variance tensor:

$$ \boldsymbol{a}(\boldsymbol{x}, t)\delta(t - t')\mathrm{d}t = \mathrm{Cov}(\boldsymbol{x}, \boldsymbol{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.
### 6.1.3 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 \stackrel{\Delta}{=} (q(\mathbf{X}_{t+dt}, t + dt) - q(\mathbf{X}_t, t))|_{\mathbf{X}_t = \mathbf{x}} = \underbrace{\mathbf{d}_t q}_{\substack{\text{Advection} \\ \text{Time increment}}} + \underbrace{(\mathbf{w}^* dt + \boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t) \cdot \nabla q}_{\substack{\text{Turbulent} \\ \text{diffusion}}} - \underbrace{\nabla \cdot (\frac{1}{2} \boldsymbol{a} \nabla q)}_{\substack{\text{Coupling} \\ \text{turbulence-forcing}}} dt + \underbrace{\mathrm{tr}((\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T)}_{\substack{\text{Turbulent} \\ \text{diffusion-forcing}}} dt, \quad (22) $$

where the quantity $q$ is forced as follows

$$ D_t q = F dt + \mathbf{H}^T \mathrm{d}\mathbf{B}_t. \quad (23) $$

The time increment $d_t q$ can be interpreted as the analog of the partial time derivative $\partial_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 Kunita, 1997; Chow, 2014).

#### Stratonovich notations:

Equations (22) and (23) rely on the so-called Ito notations. They are convenient to derive ensemble mean – since terms in $d\mathbf{B}_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 "$\circ$d$\mathbf{B}_t$" instead of "$d\mathbf{B}_t$"):

$$ D_t^S q = \underbrace{\mathbf{d}_t q}_{\substack{\text{Centered} \\ \text{time increment}}} + (\mathbf{w}_S dt + \boldsymbol{\sigma} \circ \mathrm{d}\mathbf{B}_t) \cdot \nabla q, \quad (24) $$

where

$$ D_t^S q \stackrel{\Delta}{=} \left( q \left( \mathbf{X}_{t+\frac{dt}{2}}, t+\frac{dt}{2} \right) - q \left( \mathbf{X}_{t-\frac{dt}{2}}, t-\frac{dt}{2} \right) \right) \Big|_{\mathbf{X}_{t=\mathbf{x}}}, \quad (25) $$

$$ = D_t q - \frac{1}{2} \mathrm{tr}((\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T) dt, \quad (26) $$

$$ = \left( F - \frac{1}{2} \frac{d}{dt} \langle \mathbf{H}^T, \mathbf{B}_t \rangle - \frac{1}{2} \mathrm{tr}((\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T) \right) dt + \mathbf{H}^T \circ \mathrm{d}\mathbf{B}_t, \quad (27) $$

is the Stratonovich material derivative and

$$ \mathbf{w}_S = \mathbf{w} - \frac{1}{2} \sum_{i=1}^{d} \partial_i \boldsymbol{\sigma} \boldsymbol{\sigma}_{i\bullet}^T = \mathbf{w}^* - \frac{1}{2} \boldsymbol{\sigma} (\nabla \cdot \boldsymbol{\sigma})^T, \quad (28) $$

is the Stratonovich drift of the flow (18). The Stratonovich material derivative, $D_t^S q$, coincides with the Ito material derivative, $D_t q$, for $\mathbf{H} = 0$, i.e.
### 6.1.3 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 \stackrel{\Delta}{=} (q(\mathbf{X}_{t+dt}, t+dt) - q(\mathbf{X}_t, t))|_{\mathbf{X}_t=\mathbf{x}} = \underbrace{\mathbf{d}_t q}_{\substack{\text{Advection} \\ \text{Time increment}}} + \underbrace{(\mathbf{w}^* dt + \boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t) \cdot \nabla q}_{\substack{\text{Turbulent} \\ \text{diffusion}}} - \underbrace{\nabla \cdot (\frac{1}{2} \boldsymbol{a} \nabla q)}_{\substack{\text{Coupling} \\ \text{turbulence-forcing}}} dt + \underbrace{\mathrm{tr}((\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T)}_{\substack{\text{Turbulent} \\ \text{diffusion-forcing}}} dt, \quad (22) $$

where the quantity $q$ is forced as follows

$$ D_t q = F dt + \mathbf{H}^T \mathrm{d}\mathbf{B}_t. \quad (23) $$

The time increment $d_t q$ can be interpreted as the analog of the partial time derivative $\partial_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 Kunita, 1997; Chow, 2014).

#### Stratonovich notations:

Equations (22) and (23) rely on the so-called Ito notations. They are convenient to derive ensemble mean – since terms in $\mathrm{d}\mathbf{B}_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 "$\circ$d$\mathbf{B}_t$" instead of "$\mathrm{d}\mathbf{B}_t$"):

$$ D_t^S q = \underbrace{\mathbf{d}_t q}_{\substack{\text{Centered} \\ \text{time increment}}} + (\mathbf{w}_S dt + \boldsymbol{\sigma} \circ \mathrm{d}\mathbf{B}_t) \cdot \nabla q, \quad (24) $$

where

$$ D_t^S q \stackrel{\Delta}{=} \left( q \left( \mathbf{X}_{t+\frac{dt}{2}}, t+\frac{dt}{2} \right) - q \left( \mathbf{X}_{t-\frac{dt}{2}}, t-\frac{dt}{2} \right) \right) \Big|_{\mathbf{X}_{t=\mathbf{x}}} , \quad (25) $$

$$ = D_t q - \frac{1}{2} \mathrm{tr}((\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T) dt, \quad (26) $$

$$ = \left( F - \frac{1}{2} \frac{d}{dt} \langle \mathbf{H}^T, \mathbf{B}_t \rangle - \frac{1}{2} \mathrm{tr}((\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T) \right) dt + \mathbf{H}^T \circ \mathrm{d}\mathbf{B}_t, \quad (27) $$

is the Stratonovich material derivative and

$$ \mathbf{w}_S = \mathbf{w} - \frac{1}{2} \sum_{i=1}^{d} \partial_i \boldsymbol{\sigma} \boldsymbol{\sigma}_{i\bullet}^T = \mathbf{w}^* - \frac{1}{2} \boldsymbol{\sigma} (\nabla \cdot \boldsymbol{\sigma})^T, \quad (28) $$

is the Stratonovich drift of the flow (18). The Stratonovich material derivative, $D_t^S q$, coincides with the Ito material derivative, $D_t q$, for $\mathbf{H} = 0$, i.e.
for time correlated forcing in the transport equation (23) – which is gener-
ally 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 (Kunita,
1997; Chow, 2014). The reader not familiar with stochastic partial differential
equations can also refer to section 10.1.2 of Resseguier (2017) for more de-
tails. 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).

6.1.4 Energy

Energy conservation:

To ensure a stochastic isochoric flow, incompressibility constraints on the mod-
ified drift, $\nabla \cdot w^* = 0$ and on the small-scale velocity, $\nabla \cdot \sigma = 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)):

$$
\frac{d}{dt} \int_{\Omega} \frac{1}{2} q^2 = 0. \tag{29}
$$

The noise energy intake is exactly compensated by the diffusion term. In par-
ticular, the mean energy, $\mathbb{E} \int_{\Omega} \frac{1}{2} q^2$, is also conserved. So, the LU-SALT also
meets the balance (4) of section 5.3.2 (with $\zeta = 1$, $\mathcal{L}[q] = \nabla \cdot (\frac{1}{2} a \nabla q)$ and
$(\mathcal{H}[q])^2 = \|H[q]\|^2 = \frac{1}{dt} \nabla q^T a \nabla 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:

$$
\frac{d}{dt} \int_{\Omega} \frac{1}{2} \mathrm{Var}(q) = - \frac{d}{dt} \int_{\Omega} \frac{1}{2} (\mathbb{E}(q))^2, \quad (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 straight-
forward that:

$$
D_t \left( \frac{1}{2} q^2 \right) = q D_t q = 0. \tag{31}
$$

Hence, all the discussions about the tracer $q$ are also valid for the local energy
$\frac{1}{2}q^2$.
for time correlated forcing in the transport equation (23) – which is gener-
ally 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 (Kunita,
1997; Chow, 2014). The reader not familiar with stochastic partial differential
equations can also refer to section 10.1.2 of Resseguier (2017) for more de-
tails. 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).

6.1.4 Energy

Energy conservation:

To ensure a stochastic isochoric flow, incompressibility constraints on the mod-
ified drift, $\nabla \cdot w^* = 0$ and on the small-scale velocity, $\nabla \cdot \sigma = 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)):

$$
\frac{d}{dt} \int_{\Omega} \frac{1}{2} q^2 = 0. \qquad (29)
$$

The noise energy intake is exactly compensated by the diffusion term. In par-
ticular, the mean energy, $\mathbb{E} \int_{\Omega} \frac{1}{2} q^2$, is also conserved. So, the LU-SALT also
meets the balance (4) of section 5.3.2 (with $\zeta = 1$, $\mathcal{L}[q] = \nabla \cdot (\frac{1}{2} a \nabla q)$ and
$(\mathcal{H}[q])^2 = \|H[q]\|^2 = \frac{1}{dt} \nabla q^T a \nabla 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:

$$
\frac{d}{dt} \int_{\Omega} \frac{1}{2} \mathrm{Var}(q) = - \frac{d}{dt} \int_{\Omega} \frac{1}{2} (\mathbb{E}(q))^2, \quad (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 straight-
forward that:

$$
D_t \left( \frac{1}{2} q^2 \right) = q D_t q = 0. \tag{31}
$$

Hence, all the discussions about the tracer $q$ are also valid for the local energy
$\frac{1}{2}q^2$.
Local energy in Fourier space (spectrum):

Here, $\hat{q}(\mathbf{k}, t) \triangleq \int_{\Omega} q(\mathbf{x}, t) e^{-i\mathbf{k} \cdot \mathbf{x}} d\mathbf{x}$ 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 $\Gamma(\mathbf{k}, t) = \mathbb{E}|\hat{q}(\mathbf{k}, t)|^2$. Indeed, if the small-scale velocity is homogeneous and isotropic in space, we can show (see appendix A) that:

$$ \partial_t \Gamma(\mathbf{k}, t) = \underbrace{\frac{2}{(2\pi)^d} \mathbb{E} \Re \{ \widehat{w \cdot \nabla q}(\mathbf{k}, t) \widehat{q}(\mathbf{k}, t) \}}_{\text{Usual spectrum variation}} + \underbrace{\int_{\mathbb{R}^d} F_\sigma(\mathbf{k}, \mathbf{k}', t) d\mathbf{k}',}_{\substack{\text{Noise-induced} \\ \text{spectrum flux}}} \quad (32) $$

where $\Re\{z\}$ is the real part of $z$ and the spectrum flux – induced by the small-scale velocity – from the wave-vector $\mathbf{k} - \mathbf{k}'$ toward the wave-vector $\mathbf{k}$ is:

$$ F_{\sigma}(\mathbf{k}, \mathbf{k}', t) = \frac{\|\mathbf{k}\|^2}{(2\pi)^d} \Gamma_{\sigma}(\mathbf{k}', t) \left( \sin^2((\widehat{\mathbf{k}}, \mathbf{k}')) \Gamma(\mathbf{k} - \mathbf{k}', t) - \frac{1}{d} \Gamma(\mathbf{k}, t) \right), \quad (33) $$

where

$$ \Gamma_{\sigma}(\mathbf{k}', t) = \frac{1}{dt} \mathbb{E} \| \boldsymbol{\sigma} d\widehat{\mathbf{B}}_t(\mathbf{k}', t) \|^2. \quad (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 depend on the tracer phases, but only on the tracer spectrum. Furthermore, the tracer mode $\hat{q}(\mathbf{k} - \mathbf{k}', t)$ gives energy to $\hat{q}(\mathbf{k}, t)$, by the intermediate of $\boldsymbol{\sigma} d\widehat{\mathbf{B}}_t$ if and only if:

$$ \Gamma_{\sigma}(\mathbf{k}', t) \neq 0 \text{ and } \frac{|\hat{q}(\mathbf{k}, t)|}{|\hat{q}(\mathbf{k} - \mathbf{k}', t)|} < \sqrt{d} \left| \sin((\widehat{\mathbf{k}}, \mathbf{k}')) \right|. \quad (35) $$

In this case, the energy transfer from $|\hat{q}(\mathbf{k} - \mathbf{k}', t)|^2$ toward $|\hat{q}(\mathbf{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 $\|\mathbf{k}\|$ and that the dimension $d = 2$, a backscattering of energy from the small-scale Fourier mode $\hat{q}(\mathbf{k} - \mathbf{k}', t)$ to the large-scale Fourier mode $|\hat{q}(\mathbf{k}, t)| > |\hat{q}(\mathbf{k} - \mathbf{k}', t)|$ implies:

$$ \left| \cos \left( \frac{\pi}{2} - \widehat{(\mathbf{k}, \mathbf{k}')} \right) \right| > \frac{1}{\sqrt{2}} = \cos \left( \frac{\pi}{4} \right) \quad (36) $$

$$ \text{i.e. } (\widehat{\mathbf{k}}, \mathbf{k}') \text{ modulo } \pi = - (\widehat{\mathbf{k}}, (\mathbf{k} - \mathbf{k}')) \text{ modulo } \pi \in [\frac{\pi}{4}, \frac{3\pi}{4}]. \quad (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.
Local energy in Fourier space (spectrum):

Here, $\hat{q}(\mathbf{k}, t) \triangleq \int_{\Omega} q(\mathbf{x}, t) e^{-i\mathbf{k} \cdot \mathbf{x}} d\mathbf{x}$ 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 $\Gamma(\mathbf{k}, t) = \mathbb{E}|\hat{q}(\mathbf{k}, t)|^2$. Indeed, if the small-scale velocity is homogeneous and isotropic in space, we can show (see appendix A) that:

$$ \partial_t \Gamma(\mathbf{k}, t) = \underbrace{\frac{2}{(2\pi)^d} \mathbb{E} \Re \{ \widehat{w \cdot \nabla q}(\mathbf{k}, t) \widehat{q}(\mathbf{k}, t) \}}_{\text{Usual spectrum variation}} + \underbrace{\int_{\mathbb{R}^d} F_\sigma(\mathbf{k}, \mathbf{k}', t) d\mathbf{k}',}_{\substack{\text{Noise-induced} \\ \text{spectrum flux}}} \quad (32) $$

where $\Re\{z\}$ is the real part of $z$ and the spectrum flux – induced by the small-scale velocity – from the wave-vector $\mathbf{k} - \mathbf{k}'$ toward the wave-vector $\mathbf{k}$ is:

$$ F_{\sigma}(\mathbf{k}, \mathbf{k}', t) = \frac{\|\mathbf{k}\|^2}{(2\pi)^d} \Gamma_{\sigma}(\mathbf{k}', t) \left( \sin^2((\widehat{\mathbf{k}}, \mathbf{k}')) \Gamma(\mathbf{k} - \mathbf{k}', t) - \frac{1}{d} \Gamma(\mathbf{k}, t) \right), \quad (33) $$

where

$$ \Gamma_{\sigma}(\mathbf{k}', t) = \frac{1}{dt} \mathbb{E} \| \boldsymbol{\sigma} d\widehat{\mathbf{B}}_t(\mathbf{k}', t) \|^2. \quad (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 depend on the tracer phases, but only on the tracer spectrum. Furthermore, the tracer mode $\hat{q}(\mathbf{k} - \mathbf{k}', t)$ gives energy to $\hat{q}(\mathbf{k}, t)$, by the intermediate of $\boldsymbol{\sigma} d\widehat{\mathbf{B}}_t$ if and only if:

$$ \Gamma_{\sigma}(\mathbf{k}', t) \neq 0 \text{ and } \frac{|\hat{q}(\mathbf{k}, t)|}{|\hat{q}(\mathbf{k} - \mathbf{k}', t)|} < \sqrt{d} \left| \sin((\widehat{\mathbf{k}}, \mathbf{k}')) \right|. \quad (35) $$

In this case, the energy transfer from $|\hat{q}(\mathbf{k} - \mathbf{k}', t)|^2$ toward $|\hat{q}(\mathbf{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 $\|\mathbf{k}\|$ and that the dimension $d = 2$, a backscattering of energy from the small-scale Fourier mode $\hat{q}(\mathbf{k} - \mathbf{k}', t)$ to the large-scale Fourier mode $|\hat{q}(\mathbf{k}, t)| > |\hat{q}(\mathbf{k} - \mathbf{k}', t)|$ implies:

$$ \left| \cos \left( \frac{\pi}{2} - \widehat{(\mathbf{k}, \mathbf{k}')} \right) \right| > \frac{1}{\sqrt{2}} = \cos \left( \frac{\pi}{4} \right) \quad (36) $$

$$ \text{i.e. } (\widehat{\mathbf{k}}, \mathbf{k}') \text{ modulo } \pi = - (\widehat{\mathbf{k}}, (\mathbf{k} - \mathbf{k}')) \text{ modulo } \pi \in [\frac{\pi}{4}, \frac{3\pi}{4}]. \quad (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 d$B_t$ have zero-mean, and the mean passive scalar evolution can be immediately derived taking the conditional expectation of the stochastic transport:

$$ \partial_t \mathbb{E}(q) + \underbrace{\boldsymbol{w}^* \cdot \nabla \mathbb{E}(q)}_{\text{Advection}} = \underbrace{\nabla \cdot \left( \frac{1}{2} \boldsymbol{a} \nabla \mathbb{E}(q) \right)}_{\text{Diffusion}}. \quad (38) $$

Since $\boldsymbol{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:

$$ \partial_t \mathrm{Var}(q) + \underbrace{\boldsymbol{w}^* \cdot \nabla \mathrm{Var}(q)}_{\text{Advection}} = \underbrace{\nabla \cdot \left( \frac{1}{2} \boldsymbol{a} \nabla \mathrm{Var}(q) \right)}_{\text{Diffusion}} + \underbrace{(\nabla \mathbb{E}(q))^T \boldsymbol{a} \nabla \mathbb{E}(q)}_{\text{Variance intake}}. \quad (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

$$ \frac{d}{dt} \int_{\Omega} \mathrm{Var}(q) = \int_{\Omega} (\nabla \mathbb{E}(q))^T \boldsymbol{a} \nabla \mathbb{E}(q) \geq 0. \quad (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 ensemble-based 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 (Berner et al., 2011; Gottwald and Harlim, 2013; Snyder et al., 2015).

### 6.1.5 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(\boldsymbol{x}, t)$, integrated over a material volume, $V(t)$, transported by the random flow (18):

$$ d \int_{V(t)} q = \int_{V(t)} \left[ D_t q + \nabla \cdot (\boldsymbol{w}^* dt + \boldsymbol{\sigma} d\boldsymbol{B}_t) q + d \left\langle \int_0^t D_{t'} q, \int_0^t \nabla \cdot \boldsymbol{\sigma} d\boldsymbol{B}_{t'} \right\rangle \right]. \quad (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, $\boldsymbol{\sigma}\dot{\boldsymbol{B}}$, is divergence-free and, the source of variations of the extensive property $\int_{V(t)} q$ is time-differentiable, i.e. with a differential of the form $d\int_{V(t)} q = \mathcal{F}dt$. In such a case, for an arbitrary volume, the transport theorem takes the form $D_t q = fdt$, and the material derivative can be replaced by the equation (22) to provide an intrinsic expression of this stochastic transport theorem.
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 d$B_t$ have zero-mean, and the mean passive scalar evolution can be immediately derived taking the conditional expectation of the stochastic transport:

$$ \partial_t \mathbb{E}(q) + \underbrace{\boldsymbol{w}^* \cdot \nabla \mathbb{E}(q)}_{\text{Advection}} = \underbrace{\nabla \cdot \left( \frac{1}{2} \boldsymbol{a} \nabla \mathbb{E}(q) \right)}_{\text{Diffusion}}. \quad (38) $$

Since $\boldsymbol{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:

$$ \partial_t \mathrm{Var}(q) + \underbrace{\boldsymbol{w}^* \cdot \nabla \mathrm{Var}(q)}_{\text{Advection}} = \underbrace{\nabla \cdot \left( \frac{1}{2} \boldsymbol{a} \nabla \mathrm{Var}(q) \right)}_{\text{Diffusion}} + \underbrace{(\nabla \mathbb{E}(q))^T \boldsymbol{a} \nabla \mathbb{E}(q)}_{\text{Variance intake}}. \quad (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

$$ \frac{d}{dt} \int_{\Omega} \mathrm{Var}(q) = \int_{\Omega} (\nabla \mathbb{E}(q))^T \boldsymbol{a} \nabla \mathbb{E}(q) \geq 0. \quad (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 ensemble-based 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 (Berner et al., 2011; Gottwald and Harlim, 2013; Snyder et al., 2015).

### 6.1.5 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(\boldsymbol{x}, t)$, integrated over a material volume, $V(t)$, transported by the random flow (18):

$$ d \int_{V(t)} q = \int_{V(t)} \left[ D_t q + \nabla \cdot (\boldsymbol{w}^* dt + \boldsymbol{\sigma} d\boldsymbol{B}_t) q + d \left\langle \int_0^t D_{t'} q, \int_0^t \nabla \cdot \boldsymbol{\sigma} d\boldsymbol{B}_{t'} \right\rangle \right]. \quad (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, $\boldsymbol{\sigma}\dot{\boldsymbol{B}}$, is divergence-free and, the source of variations of the extensive property $\int_{V(t)} q$ is time-differentiable, i.e. with a differential of the form $d\int_{V(t)} q = Fdt$. In such a case, for an arbitrary volume, the transport theorem takes the form $D_t q = fdt$, and the material derivative can be replaced by the equation (22) to provide an intrinsic expression of this stochastic transport theorem.
## 6.2 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 divergence-free constraints (Mémin, 2014; Resseguier et al., 2017a):

$$ \nabla \cdot w^* = 0, \qquad (42) $$

$$ \nabla \cdot \sigma = 0. \qquad (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; Resseguier et al., 2015; Harouna and Mémin, 2017; Resseguier et al., 2017a,c; Chandramouli et al., 2018; Cintolesi and Mémin, 2019a,b). Therefore, we briefly explain it here for disambiguation.

### 6.2.1 Euler LU

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

$$ D_t w = D_t^S \underbrace{w}_{\substack{\text{Due to the} \\ \text{transport} \\ \text{of } \rho w}} = \left( g - \frac{1}{\rho} \nabla p \right) dt. \qquad (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, $\int_{\Omega} \frac{1}{2} \|w\|^2$ (chapter 8 of Resseguier, 2017).

### 6.2.2 SALT Euler

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

$$ D_t^S \underbrace{w^*}_{\substack{\text{Due to the} \\ \text{transport} \\ \text{of } \rho w^*}} + \underbrace{\sum_{j=1}^d \nabla(\sigma \circ dB_t)_j w_j^*}_{\text{Additional term}} = \left( g - \frac{1}{\rho} \nabla p \right) dt. \qquad (45) $$

The reader may refer to (Holm, 2015) for a full presentation of the derivation. The SALT Euler does not conserve kinetic energy (chapter 10 of Resseguier,
## 6.2 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 $\rho$ is assumed to be constant. Using the mass conservation, this implies the following divergence-free constraints (Mémin, 2014; Resseguier et al., 2017a):

$$ \nabla \cdot w^* = 0, \qquad (42) $$

$$ \nabla \cdot \sigma = 0. \qquad (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; Resseguier et al., 2015; Harouna and Mémin, 2017; Resseguier et al., 2017a,c; Chandramouli et al., 2018; Cintolesi and Mémin, 2019a,b). Therefore, we briefly explain it here for disambiguation.

### 6.2.1 Euler LU

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

$$ D_t w = D_t^S \underbrace{w}_{\substack{\text{Due to the} \\ \text{transport} \\ \text{of } \rho w}} = \left( g - \frac{1}{\rho} \nabla p \right) dt. \qquad (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, $\int_{\Omega} \frac{1}{2} \|w\|^2$ (chapter 8 of Resseguier, 2017).

### 6.2.2 SALT Euler

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

$$ D_t^S \underbrace{w^*}_{\substack{\text{Due to the} \\ \text{transport} \\ \text{of } \rho w^*}} + \underbrace{\sum_{j=1}^d \nabla(\sigma \circ dB_t)_j w_j^*}_{\text{Additional term}} = \left( g - \frac{1}{\rho} \nabla p \right) dt. \qquad (45) $$

The reader may refer to (Holm, 2015) 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 = \int_{\Omega} w^* \cdot \nabla \times w^*$ (Holm, 2015; Resseguier et al., 2019). 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. Crisan et al. (2017) 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.

### 6.2.3 Pseudo-stochastic Euler LU

The pseudo-stochastic Euler model under location uncertainty writes:

$$ \partial_t w + (w^* \cdot \nabla) w - \nabla \cdot (\frac{1}{2} a \nabla w^T)^T = g - \frac{1}{\rho} \nabla p. \quad (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 $(\sigma dB_t \cdot \nabla) 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 $-(\sigma dB_t \cdot \nabla) 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 (Pinier et al., 2019; Resseguier et al., 2017c). This ability is a strong asset of the LU formalism for reduced order modeling.

## 6.3 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 $\sigma$ or choosing its kernel $\tilde{\sigma}$. 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.
2017) but conserves the (Stratonovich) helicity, $H = \int_{\Omega} w^* \cdot \nabla \times w^*$ (Holm, 2015; Resseguier et al., 2019). 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. Crisan et al. (2017) 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.

### 6.2.3 Pseudo-stochastic Euler LU

The pseudo-stochastic Euler model under location uncertainty writes:

$$ \partial_t w + (w^* \cdot \nabla) w - \nabla \cdot (\frac{1}{2} a \nabla w^T)^T = g - \frac{1}{\rho} \nabla p. \quad (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 $(\sigma dB_t \cdot \nabla) 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 $-(\sigma dB_t \cdot \nabla) 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 (Pinier et al., 2019; Resseguier et al., 2017c). This ability is a strong asset of the LU formalism for reduced order modeling.

## 6.3 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 $\sigma$ or choosing its kernel $\tilde{\sigma}$. 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.
### 6.3.1 Kraichnan model

Similarities between LU-SALT dynamics and the Kraichnan model (Kraichnan, 1968, 1994; Gawędzki and Kupiainen, 1995; Majda et al., 1999) were already highlighted by Mémin (2014). Kraichnan (1968) focuses on the passive scalar transport by a homogeneous ($Cov(x, y, t, t') = Cov(x - y, t, t')$), isotropic ($Cov(x, y, t, t') = Cov(\|x - y\|, t, t')$) and divergence-free velocity field uncorrelated in time. With LU notations, it implies:

$$w = 0, \qquad (47)$$

$$\nabla \cdot \sigma = 0, \qquad (48)$$

$$\sigma d\mathbf{B}_t = \tilde{\sigma} * d\mathbf{B}_t, \qquad (49)$$

where * denotes a convolution and $\tilde{\sigma}$ is a particular (isotropic) filter. The homogeneity condition (49) makes the parameterization easier in Fourier space:

$$\mathbb{E} \left\{ \left( \widehat{\sigma d\mathbf{B}_t} \right) (\mathbf{k}_1, t) \left( \widehat{\sigma d\mathbf{B}_{t'}} \right)^T (\mathbf{k}_2, t') \right\} = \Gamma_{\sigma}(\mathbf{k}_1) \delta(\mathbf{k}_1 - \mathbf{k}_2) \delta(t - t') dt, (50)$$

where $\Gamma_{\sigma}$ can be taken to mimic the self-similarities observed in real turbulence flow (Gawędzki and Kupiainen, 1995):

$$\Gamma_{\sigma}(\mathbf{k}) = A^2 \hat{\mathcal{P}}(\mathbf{k}) \left( 1 + \frac{\|\mathbf{k}\|^2}{\kappa_m^2} \right)^{-2\alpha}, \qquad (51)$$

with A is a constant, $1/\kappa_m$ specifies the velocity correlation length and

$$\hat{\mathcal{P}}(\mathbf{k}) = \mathbb{I}_d - \frac{\mathbf{k}\mathbf{k}^T}{\|\mathbf{k}\|^2}, \qquad (52)$$

is the projection – expressed in Fourier space – onto the space of solenoidal functions. The Fourier transform of the associated filter $\tilde{\sigma}$ (appearing in (49)) is:

$$\hat{\tilde{\sigma}}(\mathbf{k}) = A \hat{\mathcal{P}}(\mathbf{k}) \left( 1 + \frac{\|\mathbf{k}\|^2}{\kappa_m^2} \right)^{-\alpha}. \qquad (53)$$

Finally, the Fourier transform of the small-scale velocity writes:

$$\hat{\sigma \tilde{B}(\mathbf{k})} = \frac{A}{\sqrt{\Delta t}} \hat{\mathcal{P}}(\mathbf{k}) \left( 1 + \frac{\|\mathbf{k}\|^2}{\kappa_m^2} \right)^{-\alpha} \hat{\frac{d\mathbf{B}_t}{\sqrt{\Delta t}}} (\mathbf{k}), \qquad (54)$$

where $\hat{d\mathbf{B}_t}$ is the spatial Fourier transform of $d\mathbf{B}_t$, with $d\mathbf{B}_t/\sqrt{\Delta 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 $d\mathbf{B}_t/\sqrt{\Delta t}$, to get $\hat{d\mathbf{B}_t}/\sqrt{\Delta t}$, and finally $\hat{\sigma \tilde{B}(\mathbf{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 $\hat{\mathcal{M}}(\mathbf{k}) \triangleq (1 + \frac{\|\mathbf{k}\|^2}{\kappa_m^2})^{-\alpha}$ is
### 6.3.1 Kraichnan model

Similarities between LU-SALT dynamics and the Kraichnan model (Kraichnan, 1968, 1994; Gawędzki and Kupiainen, 1995; Majda et al., 1999) were already highlighted by Mémin (2014). Kraichnan (1968) focuses on the passive scalar transport by a homogeneous ($Cov(x, y, t, t') = Cov(x - y, t, t')$), isotropic ($Cov(x, y, t, t') = Cov(\|x - y\|, t, t')$) and divergence-free velocity field uncorrelated in time. With LU notations, it implies:

$$w = 0, \qquad (47)$$

$$\nabla \cdot \sigma = 0, \qquad (48)$$

$$\sigma d\mathbf{B}_t = \tilde{\sigma} * d\mathbf{B}_t, \qquad (49)$$

where * denotes a convolution and $\tilde{\sigma}$ is a particular (isotropic) filter. The homogeneity condition (49) makes the parameterization easier in Fourier space:

$$\mathbb{E} \left\{ \left( \widehat{\sigma d\mathbf{B}_t} \right) (\mathbf{k}_1, t) \left( \widehat{\sigma d\mathbf{B}_{t'}} \right)^T (\mathbf{k}_2, t') \right\} = \Gamma_{\sigma}(\mathbf{k}_1) \delta(\mathbf{k}_1 - \mathbf{k}_2) \delta(t - t') dt, (50)$$

where $\Gamma_{\sigma}$ can be taken to mimic the self-similarities observed in real turbulence flow (Gawędzki and Kupiainen, 1995):

$$\Gamma_{\sigma}(\mathbf{k}) = A^2 \hat{\mathcal{P}}(\mathbf{k}) \left( 1 + \frac{\|\mathbf{k}\|^2}{\kappa_m^2} \right)^{-2\alpha}, \qquad (51)$$

with A is a constant, $1/\kappa_m$ specifies the velocity correlation length and

$$\hat{\mathcal{P}}(\mathbf{k}) = \mathbb{I}_d - \frac{\mathbf{k}\mathbf{k}^T}{\|\mathbf{k}\|^2}, \qquad (52)$$

is the projection – expressed in Fourier space – onto the space of solenoidal functions. The Fourier transform of the associated filter $\tilde{\sigma}$ (appearing in (49)) is:

$$\hat{\tilde{\sigma}}(\mathbf{k}) = A \hat{\mathcal{P}}(\mathbf{k}) \left( 1 + \frac{\|\mathbf{k}\|^2}{\kappa_m^2} \right)^{-\alpha}. \qquad (53)$$

Finally, the Fourier transform of the small-scale velocity writes:

$$\hat{\sigma \tilde{B}(\mathbf{k})} = \frac{A}{\sqrt{\Delta t}} \hat{\mathcal{P}}(\mathbf{k}) \left( 1 + \frac{\|\mathbf{k}\|^2}{\kappa_m^2} \right)^{-\alpha} \hat{\frac{d\mathbf{B}_t}{\sqrt{\Delta t}}} (\mathbf{k}), \qquad (54)$$

where $\hat{d\mathbf{B}_t}$ is the spatial Fourier transform of $d\mathbf{B}_t$, with $d\mathbf{B}_t/\sqrt{\Delta 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 $d\mathbf{B}_t/\sqrt{\Delta t}$, to get $\hat{d\mathbf{B}_t}/\sqrt{\Delta t}$, and finally $\hat{\sigma \tilde{B}(\mathbf{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 $\hat{\mathcal{M}}(\mathbf{k}) \triangleq (1 + \frac{\|\mathbf{k}\|^2}{\kappa_m^2})^{-\alpha}$ is
the Fourier transform of a Matérn covariance $\mathcal{M}$ of degree $\nu = \alpha - 1$ and of range $\rho = \sqrt{2(\alpha - 1)}/\kappa_m$. Therefore, the unresolved velocity is:

$$ \sigma d\mathbf{B}_t = AP\{\mathcal{M} * d\mathbf{B}_t\}, \qquad (55) $$

with

$$ \mathcal{P} = \mathbb{I}_d - \Delta^{-1}\nabla\nabla^T. \qquad (56) $$

### 6.3.2 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 velocity leads to a heterogeneous non-Gaussian random forcing, due to its multiplicative structure.

Rességuier 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 (Blumen, 1978; Held et al., 1995; Lapeyre, 2017). For this purpose, a modification of (53) is considered. Instead of involving the operator $\mathcal{P} = \mathbb{I}_d - \Delta^{-1}\nabla\nabla^T$, Rességuier et al. (2017b) work on the streamfunction $\psi_\sigma d\mathbf{B}_t$. As such no additional divergence-free constraint is needed. The streamfunction kernel, $\tilde{\psi}_\sigma$ is defined such as:

$$ \sigma(\mathbf{x})d\mathbf{B}_t = \nabla^\perp \psi_\sigma(\mathbf{x})d\mathbf{B}_t = (\nabla^\perp \tilde{\psi}_\sigma \star d\mathbf{B}_t)(\mathbf{x}). \qquad (57) $$

with $\nabla^\perp = (-\partial_y, \partial_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 $\mathbb{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:

$$ \mathbf{k} \mapsto \frac{1}{\mu(\Omega)} \mathbb{E} \oint_{[0,2\pi]} d\theta_k \|\mathbf{k}\| \left\| \widehat{\sigma \mathbf{B}(\mathbf{k})} \right\|^2 = \frac{2\pi}{\Delta t} \|\mathbf{k}\|^3 \left| \widehat{\psi_\sigma}(\|\mathbf{k}\|) \right|^2, \qquad (58) $$

where $\mu(\Omega)$ is the surface of the spatial domain $\Omega$, $\theta_k$ is the angle of the wave-vector $\mathbf{k}$ and $\Delta 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 $\kappa_m$ and $\kappa_M$. The parameter $\kappa_m$ is inversely related to the spatial correlation length of the unresolved component. $\kappa_M$ is set to the theoretical resolution, $\pi/\Delta x$, and $\kappa_m$ to the estimated effective
the Fourier transform of a Matérn covariance $\mathcal{M}$ of degree $\nu = \alpha - 1$ and of range $\rho = \sqrt{2(\alpha - 1)}/\kappa_m$. Therefore, the unresolved velocity is:

$$ \sigma d\mathbf{B}_t = AP\{\mathcal{M} * d\mathbf{B}_t\}, \qquad (55) $$

with

$$ \mathcal{P} = \mathbb{I}_d - \Delta^{-1}\nabla\nabla^T. \qquad (56) $$

### 6.3.2 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 velocity leads to a heterogeneous non-Gaussian random forcing, due to its multiplicative structure.

Rességuier 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 (Blumen, 1978; Held et al., 1995; Lapeyre, 2017). For this purpose, a modification of (53) is considered. Instead of involving the operator $\mathcal{P} = \mathbb{I}_d - \Delta^{-1}\nabla\nabla^T$, Rességuier et al. (2017b) work on the streamfunction $\psi_\sigma d\mathbf{B}_t$. As such no additional divergence-free constraint is needed. The streamfunction kernel, $\tilde{\psi}_\sigma$ is defined such as:

$$ \sigma(\mathbf{x})d\mathbf{B}_t = \nabla^\perp \psi_\sigma(\mathbf{x})d\mathbf{B}_t = (\nabla^\perp \tilde{\psi}_\sigma \star d\mathbf{B}_t)(\mathbf{x}). \qquad (57) $$

with $\nabla^\perp = (-\partial_y, \partial_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 $\mathbb{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:

$$ \mathbf{k} \mapsto \frac{1}{\mu(\Omega)} \mathbb{E} \oint_{[0,2\pi]} d\theta_k \|\mathbf{k}\| \left\| \widehat{\sigma \mathbf{B}(\mathbf{k})} \right\|^2 = \frac{2\pi}{\Delta t} \|\mathbf{k}\|^3 \left| \widehat{\tilde{\psi}_\sigma}(\|\mathbf{k}\|) \right|^2, \qquad (58) $$

where $\mu(\Omega)$ is the surface of the spatial domain $\Omega$, $\theta_k$ is the angle of the wave-vector $\mathbf{k}$ and $\Delta 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 $\kappa_m$ and $\kappa_M$. The parameter $\kappa_m$ is inversely related to the spatial correlation length of the unresolved component. $\kappa_M$ is set to the theoretical resolution, $\pi/\Delta x$, and $\kappa_m$ to the estimated effective
resolution ($\kappa_m = \kappa_M/2$ in Resseguier et al. (2017b)). The small scales' energy is specified by the diffusion coefficient $a$ and the simulation time step:

$$ \mathbb{E} (\boldsymbol{\sigma} \dot{\boldsymbol{B}}) (\boldsymbol{\sigma} \dot{\boldsymbol{B}})^T = \frac{1}{\Delta t} a = \frac{1}{\Delta t} \begin{pmatrix} a_0 & 0 \\ 0 & a_0 \end{pmatrix}. \qquad (59) $$

The diagonal structure of the variance tensor is due both to incompressibility 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:

$$ \widehat{\boldsymbol{\sigma} \dot{\boldsymbol{B}}}(\boldsymbol{k}) \triangleq \frac{A}{\sqrt{\Delta t}} i \boldsymbol{k}^{\perp} f_{\text{BP}} (\|\boldsymbol{k}\|) \|\boldsymbol{k}\|^{-\alpha} \frac{\widehat{d B_t}}{\sqrt{\Delta t}}(\boldsymbol{k}) \quad \text{with } s = 3 - 2\alpha = -\frac{5}{3}, \quad (60) $$

where $A$ is a constant to ensure $\mathbb{E}\|\boldsymbol{\sigma} \dot{\boldsymbol{B}}\|^2 = \text{tr}(\boldsymbol{a})/(\Delta t) = 2a_0/(\Delta t)$ (see (59) above), $\widehat{d B_t}$ is the spatial Fourier transform of $d B_t$, with $d B_t/\sqrt{\Delta t}$, a discrete scalar white noise process of unit variance in space and time. To sample the small-scale velocity, we first sample $d B_t/\sqrt{\Delta t}$, to get $\widehat{d B_t}/\sqrt{\Delta t}$, and finally $\widehat{\boldsymbol{\sigma} \dot{\boldsymbol{B}}}(\boldsymbol{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. (2017b) have numerically demonstrated the good UQ skills of this method for SQG flows.

### 6.3.3 The homogeneous non-stationary model

A main drawback of the previous choice of $\boldsymbol{\sigma}$ is the needed tuning of the parameters. Resseguier et al. (2019) propose a more general homogeneous parameterization where no tuning is needed. A new $\boldsymbol{\sigma}$ 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 (Falkovich et al., 2001; Penland, 2003b; Klyatskin, 2005; Vallis, 2006; Keating et al., 2011)) of the unresolved velocity is twice the variance tensor trace $\text{tr}(\boldsymbol{a})$ whereas the unresolved kinetic energy is $\text{tr}(\boldsymbol{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
resolution ($\kappa_m = \kappa_M/2$ in Resseguier et al. (2017b)). The small scales' energy is specified by the diffusion coefficient $a$ and the simulation time step:

$$ \mathbb{E} (\boldsymbol{\sigma} \dot{\boldsymbol{B}}) (\boldsymbol{\sigma} \dot{\boldsymbol{B}})^T = \frac{1}{\Delta t} a = \frac{1}{\Delta t} \begin{pmatrix} a_0 & 0 \\ 0 & a_0 \end{pmatrix}. \qquad (59) $$

The diagonal structure of the variance tensor is due both to incompressibility 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:

$$ \widehat{\boldsymbol{\sigma}\dot{\boldsymbol{B}}}(\boldsymbol{k}) \triangleq \frac{A}{\sqrt{\Delta t}} i\boldsymbol{k}^{\perp} f_{\text{BP}} (\|\boldsymbol{k}\|) \|\boldsymbol{k}\|^{-\alpha} \frac{\widehat{dB_t}}{\sqrt{\Delta t}}(\boldsymbol{k}) \quad \text{with } s = 3 - 2\alpha = -\frac{5}{3}, \quad (60) $$

where $A$ is a constant to ensure $\mathbb{E}\|\boldsymbol{\sigma}\dot{\boldsymbol{B}}\|^2 = \text{tr}(\boldsymbol{a})/(\Delta t) = 2a_0/(\Delta t)$ (see (59) above), $\widehat{dB_t}$ is the spatial Fourier transform of $dB_t$, with $dB_t/\sqrt{\Delta 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/\sqrt{\Delta t}$, to get $\widehat{dB_t}/\sqrt{\Delta t}$, and finally $\widehat{\boldsymbol{\sigma}\dot{\boldsymbol{B}}}(\boldsymbol{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. (2017b) have numerically demonstrated the good UQ skills of this method for SQG flows.

### 6.3.3 The homogeneous non-stationary model

A main drawback of the previous choice of $\boldsymbol{\sigma}$ is the needed tuning of the parameters. Resseguier et al. (2019) propose a more general homogeneous parameterization where no tuning is needed. A new $\boldsymbol{\sigma}$ 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 (Falkovich et al., 2001; Penland, 2003b; Klyatskin, 2005; Vallis, 2006; Keating et al., 2011)) of the unresolved velocity is twice the variance tensor trace $\text{tr}(\boldsymbol{a})$ whereas the unresolved kinetic energy is $\text{tr}(\boldsymbol{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(\kappa)$ at the wave-number $\kappa$. We will rely on this ADSD rather than on the KE spectrum, $E(\kappa)$. Since the absolute diffusivity is the variance multiplied by the correlation time, it is naturally to defined the ADSD as follows:

$$A(\kappa) \triangleq E(\kappa)\tau(\kappa) \quad \text{where} \quad \tau(k) = \frac{1/\kappa}{v_{\kappa}} = \frac{1/\kappa}{\sqrt{\kappa E(\kappa)}}, \qquad (61)$$

is the eddy turnover time at the scale $1/\kappa$ and $v_\kappa$ is at characteristic velocity at this scale. Accordingly, we have:

$$A(\kappa) = \kappa^{-3/2} E^{1/2}(\kappa). \qquad (62)$$

If in addition we assume a KE self-similar distribution,

$$E(\kappa) = C^2 \kappa^{-s}, \qquad (63)$$

we obtain:

$$A(\kappa) = C\kappa^{-r}, \qquad (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 large-scale 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 an inertial range of fully-resolved scales (i.e. before the spectrum roll-off).

From there, we can define the operator $\sigma$ 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, $\sigma$ is chosen such as the unresolved ADSD compensates the resolved ADSD roll off – introduced by the deterministic subgrid parameterization – at small scales. Specifically, the unresolved ADSD is set to:

$$A_{\sigma \dot{B}}(\kappa) = \max(0, C_w \kappa^{-r_w} - A_w(\kappa)) f_{BP}^2(\kappa). \qquad (65)$$

As previously, $f_{BP}^2$ is a band-pass filter between $\kappa_m$ and $\kappa_M$. In practice, we set $\kappa_M$ to the theoretical resolution:

$$\kappa_M = \frac{\pi}{\Delta x}, \qquad (66)$$

and $\kappa_m$ to the effective resolution which is estimated as follows:

$$\kappa_m = \left( \frac{\ln(0.95)}{\ln(0.1)} \right)^{1/p} \kappa_M, \qquad (67)$$
(ADSD) denoted $A(\kappa)$ at the wave-number $\kappa$. We will rely on this ADSD rather than on the KE spectrum, $E(\kappa)$. Since the absolute diffusivity is the variance multiplied by the correlation time, it is naturally to defined the ADSD as follows:

$$A(\kappa) \triangleq E(\kappa)\tau(\kappa) \quad \text{where} \quad \tau(k) = \frac{1/\kappa}{v_\kappa} = \frac{1/\kappa}{\sqrt{\kappa E(\kappa)}}, \qquad (61)$$

is the eddy turnover time at the scale $1/\kappa$ and $v_\kappa$ is at characteristic velocity at this scale. Accordingly, we have:

$$A(\kappa) = \kappa^{-3/2} E^{1/2}(\kappa). \qquad (62)$$

If in addition we assume a KE self-similar distribution,

$$E(\kappa) = C^2 \kappa^{-s}, \qquad (63)$$

we obtain:

$$A(\kappa) = C\kappa^{-r}, \qquad (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 large-scale 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 an inertial range of fully-resolved scales (i.e. before the spectrum roll-off).

From there, we can define the operator $\sigma$ 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, $\sigma$ is chosen such as the unresolved ADSD compensates the resolved ADSD roll off – introduced by the deterministic subgrid parameterization – at small scales. Specifically, the unresolved ADSD is set to:

$$A_{\sigma \dot{B}}(\kappa) = \max(0, C_w \kappa^{-r_w} - A_w(\kappa)) f_{BP}^2(\kappa). \qquad (65)$$

As previously, $f_{BP}^2$ is a band-pass filter between $\kappa_m$ and $\kappa_M$. In practice, we set $\kappa_M$ to the theoretical resolution:

$$\kappa_M = \frac{\pi}{\Delta x}, \qquad (66)$$

and $\kappa_m$ to the effective resolution which is estimated as follows:

$$\kappa_m = \left( \frac{\ln(0.95)}{\ln(0.1)} \right)^{1/p} \kappa_M, \qquad (67)$$
where $p$ is the order of the Laplacian used as deterministic subgrid tensor (i.e. $D_t b = -\nu(-\Delta)^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 $\kappa_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-numbers – scales smaller than $\kappa_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_{\sigma \hat{B}}$ through (65) at large wave-numbers – larger than $\kappa_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 $\hat{\sigma}$ which defines the unresolved velocity, we note that

$$A_{\sigma \hat{B}}(\kappa) = E_{\sigma \hat{B}}(\kappa) dt = \frac{1}{\mu(\Omega)} \oint_{[0,2\pi]} d\theta_k \kappa \left\| \hat{\sigma}(t, \kappa) \right\|^2 = 2\pi\kappa^3 \left| \hat{\psi}_{\sigma}(t, \kappa) \right|^2 \quad (68)$$

From formulas (65)-(68), we can finally express the unresolved velocity as follows:

$$\hat{\sigma \hat{B}}(t, \mathbf{k}) = \frac{1}{\sqrt{\Delta t}} i \mathbf{k}^{\perp} \hat{\psi}_{\sigma}(t, \|\mathbf{k}\|) \frac{\mathrm{d}B_t}{\sqrt{\Delta t}}(\mathbf{k}), \quad (69)$$

$$= \frac{1}{\sqrt{\Delta t}} i \mathbf{k}^{\perp} \sqrt{\frac{\max(0, C_w \|\mathbf{k}\|^{-r_w} - A_w(\|\mathbf{k}\|))}{2\pi \|\mathbf{k}\|^3}} f_{BP}(\|\mathbf{k}\|) \frac{\mathrm{d}B_t}{\sqrt{\Delta t}}(\mathbf{k}). \quad (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 $\tau(\kappa)$ is not restricted to the time step $dt$.

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

### 6.3.4 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
where $p$ is the order of the Laplacian used as deterministic subgrid tensor (i.e. $D_t b = -\nu(-\Delta)^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 $\kappa_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-numbers – scales smaller than $\kappa_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_{\sigma \hat{B}}$ through (65) at large wave-numbers – larger than $\kappa_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 $\hat{\sigma}$ which defines the unresolved velocity, we note that

$$A_{\sigma \hat{B}}(\kappa) = E_{\sigma \hat{B}}(\kappa) dt = \frac{1}{\mu(\Omega)} \oint_{[0,2\pi]} d\theta_k \kappa \left\| \hat{\sigma}(t, \kappa) \right\|^2 = 2\pi\kappa^3 \left| \hat{\psi}_{\sigma}(t, \kappa) \right|^2 \quad (68)$$

From formulas (65)-(68), we can finally express the unresolved velocity as follows:

$$\hat{\sigma \hat{B}}(t, \mathbf{k}) = \frac{1}{\sqrt{\Delta t}} i \mathbf{k}^{\perp} \hat{\psi}_{\sigma}(t, \|\mathbf{k}\|) \frac{\mathrm{d}B_t}{\sqrt{\Delta t}}(\mathbf{k}), \quad (69)$$

$$= \frac{1}{\sqrt{\Delta t}} i \mathbf{k}^{\perp} \sqrt{\frac{\max(0, C_w \|\mathbf{k}\|^{-r_w} - A_w(\|\mathbf{k}\|))}{2\pi \|\mathbf{k}\|^3}} f_{BP}(\|\mathbf{k}\|) \frac{\mathrm{d}B_t}{\sqrt{\Delta t}}(\mathbf{k}). \quad (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 $\tau(\kappa)$ is not restricted to the time step $dt$.

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

### 6.3.4 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 (Ferrari and Nikurashin, 2010), 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 $\kappa$ cannot depend on 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 spatially-homogeneous spectrum slope. The constant may also varies with the time and the wave number. According to the Kolmogorov theory (e.g. Frisch, 1995) and (63):

$$C_w(x, t, \kappa) = \text{cst.} \epsilon_F^p(x, \kappa), \qquad (71)$$

where $\epsilon_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 larger than $1/\kappa$ toward smaller scales and can be computed as follows:

$$\epsilon_F(x, t, \kappa) \triangleq \overline{q_{\kappa}^{<}((\overline{w \cdot \nabla})q)_{\kappa}^{<}}, \qquad (72)$$

where $q$ is the transported (up to possible source terms) quantity, $g_{\kappa}^{<}$ is the low-pass filtered version of $g$ (setting to zero the Fourier modes of $g$ which have frequencies larger than $\kappa$) and $\cdot$ stands for the spatial average (Frisch, 1995). 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 (Frisch, 1995).

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 $\epsilon_F^p/\epsilon_F^{\overline{p}}$ – averaged over the
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 (Ferrari and Nikurashin, 2010), 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 $\kappa$ cannot depend on 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 spatially-homogeneous spectrum slope. The constant may also varies with the time and the wave number. According to the Kolmogorov theory (e.g. Frisch, 1995) and (63):

$$C_w(x, t, \kappa) = \text{cst.} \epsilon_F^p(x, \kappa), \quad (71)$$

where $\epsilon_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 larger than $1/\kappa$ toward smaller scales and can be computed as follows:

$$\epsilon_F(x, t, \kappa) \triangleq \overline{q_{\kappa}^{<}((\overline{w \cdot \nabla})q)_{\kappa}^{<}}, \quad (72)$$

where $q$ is the transported (up to possible source terms) quantity, $g_{\kappa}^{<}$ is the low-pass filtered version of $g$ (setting to zero the Fourier modes of $g$ which have frequencies larger than $\kappa$) and $\cdot$ stands for the spatial average (Frisch, 1995). 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 (Frisch, 1995).

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 $\epsilon_F^p/\epsilon_F^{\overline{p}}$ – averaged over the
resolved inertial range wave-numbers:

$$ \widetilde{\sigma \dot{B}}(t, \mathbf{k}) = \frac{1}{\sqrt{\Delta t}} i \mathbf{k}^{\perp} \sqrt{\frac{\max(0, C_w \|\mathbf{k}\|^{-r_w} - A_w(\|\mathbf{k}\|))}{2\pi \|\mathbf{k}\|^3}} f_{BP}(\|\mathbf{k}\|) \frac{\widehat{dB_t}(\mathbf{k})}{\sqrt{\Delta t}} \quad (73) $$

$$ \sigma \dot{B}(t, \mathbf{x}) = \mathcal{P} \left\{ \sqrt{\underbrace{\frac{\epsilon_F^p(t, \mathbf{x})}{\epsilon_F^p(t)}}_{\text{Heterogeneous modulation}} \underbrace{\widetilde{\sigma \dot{B}}(t, \mathbf{x})}_{\text{Homogeneous velocity}}} \right\}, \quad (74) $$

where $\mathcal{P} = \mathbb{I}_d - \Delta^{-1}\nabla\nabla^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 $\epsilon_F$ and thus larger $C_w$) suggests that the unresolved velocity (large) should maintain 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 $\mathcal{P}$. Because of that we do not consider the advection correction $\mathbf{w}^* - \mathbf{w} = -\frac{1}{2}(\nabla \cdot \mathbf{a})^T$ of the LU formalism. Indeed, here the variance tensor has the simple form $\mathbf{a} = \frac{1}{d}\text{tr}(\mathbf{a})\mathbb{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.

Resseguier et al. (2019) have numerically shown that this parameterization indeed keeps the sharp straight fronts while perturbing the eddies and the meanders in SQG flows.

### 6.3.5 The heterogeneous spatially-uncorrelated model

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

### 6.3.6 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):

$$ \sigma(\mathbf{x}, t) \dot{\mathbf{B}} = \sum_{k=1}^{n_{KL}} \xi_k(\mathbf{x}, t) \dot{W}_k, \quad (75) $$

with $W_k$ are independent one-dimensional Brownian motions and $n_{KL}$ possibly infinite. The $(\xi_k)_k$ are so-called (weighted) Empirical Orthogonal Functions
resolved inertial range wave-numbers:

$$ \widetilde{\sigma \dot{B}}(t, \mathbf{k}) = \frac{1}{\sqrt{\Delta t}} i \mathbf{k}^{\perp} \sqrt{\frac{\max(0, C_w \|\mathbf{k}\|^{-r_w} - A_w(\|\mathbf{k}\|))}{2\pi \|\mathbf{k}\|^3}} f_{BP}(\|\mathbf{k}\|) \frac{\widehat{dB_t}(\mathbf{k})}{\sqrt{\Delta t}} \quad (73) $$

$$ \sigma \dot{B}(t, \mathbf{x}) = \mathcal{P} \left\{ \sqrt{\underbrace{\frac{\epsilon_F^p(t, \mathbf{x})}{\epsilon_F^p(t)}}_{\text{Heterogeneous modulation}} \underbrace{\widetilde{\sigma \dot{B}}(t, \mathbf{x})}_{\text{Homogeneous velocity}}} \right\}, \quad (74) $$

where $\mathcal{P} = \mathbb{I}_d - \Delta^{-1}\nabla\nabla^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 $\epsilon_F$ and thus larger $C_w$) suggests that the unresolved velocity (large) should maintain 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 $\mathcal{P}$. Because of that we do not consider the advection correction $w^* - w = -\frac{1}{2}(\nabla \cdot \mathbf{a})^T$ of the LU formalism. Indeed, here the variance tensor has the simple form $\mathbf{a} = \frac{1}{d}\text{tr}(\mathbf{a})\mathbb{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.

Resseguier et al. (2019) have numerically shown that this parameterization indeed keeps the sharp straight fronts while perturbing the eddies and the meanders in SQG flows.

### 6.3.5 The heterogeneous spatially-uncorrelated model

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

### 6.3.6 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):

$$ \sigma(\mathbf{x}, t) \dot{\mathbf{B}} = \sum_{k=1}^{n_{KL}} \xi_k(\mathbf{x}, t) \dot{W}_k, \quad (75) $$

with $W_k$ are independent one-dimensional Brownian motions and $n_{KL}$ possibly infinite. The $(\xi_k)_k$ are so-called (weighted) Empirical Orthogonal Functions
(EOFs). If we denotes by $\lambda_k$ the $L^2$ norm of the EOFs, their normalized versions $(\tilde{\xi}_k = \xi_k/\lambda_k)_k$ are the eigenvectors of the auto-adjoint operator defined by the kernel $\sigma(\mathbf{x}, t)\sigma^T(\mathbf{y}, t)$. The $(\lambda_k)_k$ are the corresponding eigenvalues. That kernel corresponds – up to the factor $\frac{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 $\sigma\sigma^T$-Wiener processes (Da Prato and Zabczyk, 1992; Prévôt and Röckner, 2007).

First, we will consider the off-line learning of the EOFs $(\xi_k)_k$ from a velocity field set. Specifically, we will assume that we have a history of “small-scale” velocity fields $(\mathbf{v}'(\mathbf{x}_{ij}, m))_{1\le m\le n_O}$ at $n_O$ different times $t_m$ (or $n_O$ different “realizations”) sampled on a spatial grid of $M$ points $(\mathbf{x}_{ij})_{ij} = ((i\Delta x, j\Delta 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, $(\xi_i)_i$, we assume that we can observe increments

$$ \int_{t_m}^{t_m+\Delta T} \sigma \mathrm{d}\mathbf{B}_{t'} = \sum_{k=1}^{n_{KL}} \xi_k(\mathbf{x}_{ij})(W_k(t_m + \Delta T) - W_k(t_m)). \quad (76) $$

at several times $t_m$. Here, $\Delta T$ is the time step of a low-resolution simulation. We will interpret the following snapshots:

$$ \widetilde{\Delta \mathbf{X}}_{ij}^m \triangleq \mathbf{v}'(\mathbf{x}_{ij}, m) \Delta T, \quad (77) $$

as a realizations of the flow increments (76). Since Brownian time increments (76) are independent, by the above interpretation, we implicitly assume that $\mathbf{v}'(\mathbf{x}_{ij}, m)$ are independent for different values of $m$. In this case, it is equivalent to consider them 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 $(\nabla \cdot \xi_i = 0)$. Therefore, after computing the residual flow increments $(\Delta \widetilde{\mathbf{X}}_{ij}^m)_m$, they are centered:

$$ \Delta \mathbf{X}_{ij}^m \triangleq \Delta \widetilde{\mathbf{X}}_{ij}^m - \hat{\mathbb{E}}\{\Delta \widetilde{\mathbf{X}}_{ij}^m\}, \quad (78) $$
(EOFs). If we denotes by $\lambda_k$ the $L^2$ norm of the EOFs, their normalized versions $(\tilde{\xi}_k = \xi_k/\lambda_k)_k$ are the eigenvectors of the auto-adjoint operator defined by the kernel $\sigma(\mathbf{x}, t)\sigma^T(\mathbf{y}, t)$. The $(\lambda_k)_k$ are the corresponding eigenvalues. That kernel corresponds – up to the factor $\frac{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 $\sigma\sigma^T$-Wiener processes (Da Prato and Zabczyk, 1992; Prévôt and Röckner, 2007).

First, we will consider the off-line learning of the EOFs $(\xi_k)_k$ from a velocity field set. Specifically, we will assume that we have a history of “small-scale” velocity fields $(\mathbf{v}'(\mathbf{x}_{ij}, m))_{1\le m\le n_O}$ at $n_O$ different times $t_m$ (or $n_O$ different “realizations”) sampled on a spatial grid of $M$ points $(\mathbf{x}_{ij})_{ij} = ((i\Delta x, j\Delta 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, $(\xi_i)_i$, we assume that we can observe increments

$$ \int_{t_m}^{t_m+\Delta T} \sigma \mathrm{d}\mathbf{B}_{t'} = \sum_{k=1}^{n_{KL}} \xi_k(\mathbf{x}_{ij})(W_k(t_m + \Delta T) - W_k(t_m)). \quad (76) $$

at several times $t_m$. Here, $\Delta T$ is the time step of a low-resolution simulation. We will interpret the following snapshots:

$$ \widetilde{\Delta \mathbf{X}}_{ij}^m \triangleq \mathbf{v}'(\mathbf{x}_{ij}, m) \Delta T, \quad (77) $$

as a realizations of the flow increments (76). Since Brownian time increments (76) are independent, by the above interpretation, we implicitly assume that $\mathbf{v}'(\mathbf{x}_{ij}, m)$ are independent for different values of $m$. In this case, it is equivalent to consider them 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 $(\nabla \cdot \xi_i = 0)$. Therefore, after computing the residual flow increments $(\widetilde{\Delta \mathbf{X}}_{ij}^m)_m$, they are centered:

$$ \Delta \mathbf{X}_{ij}^m \triangleq \widetilde{\Delta \mathbf{X}}_{ij}^m - \hat{\mathbb{E}}\{\widetilde{\Delta \mathbf{X}}_{ij}^m\}, \quad (78) $$
with the estimator

$$ \hat{\mathbb{E}}\{\widetilde{\Delta\mathbf{X}}_{ij}^m\} \triangleq \frac{1}{n_O} \sum_{m=1}^{n_O} \widetilde{\Delta\mathbf{X}}_{ij}^m, \qquad (79) $$

and projected onto the space of divergence-free functions:

$$ \Delta\mathbf{X}_{ij}^m \triangleq \mathcal{P}\{\Delta\mathbf{X}_{ij}^{m'}\} \quad \text{with} \quad \mathcal{P} = \mathbb{I}_d - \Delta^{-1}\nabla\nabla^T. \qquad (80) $$

Note that the projection can also be applied directly on the final EOFs $\xi_k$ – at the end of the estimation process – rather than on each increments $\Delta\mathbf{X}_{ij}^m$. Besides, if $\nabla \cdot v' = 0$ then the centred flow increments are already divergence-free 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$):

$$ \gamma_{\sigma}(\mathbf{x}_{ij}, \mathbf{x}_{pq}) \approx \frac{1}{(n_O - 1)\Delta T} \sum_{m=1}^{n_O} \Delta\mathbf{X}_{ij}^m (\Delta\mathbf{X}_{pq}^m)^T. \qquad (81) $$

Besides, we have:

$$ \gamma_{\sigma}(\mathbf{x}_{ij}, \mathbf{x}_{pq}) = \sum_{k=1}^{n_{KL}} \xi_k(\mathbf{x}_{ij}) \xi_k^T(\mathbf{x}_{pq}). \qquad (82) $$

In order to properly define the EOFs, we must add the following orthogonal constraint:

$$ \sum_{pq} \xi_i(\mathbf{x}_{pq}) \cdot \xi_j(\mathbf{x}_{pq}) = 0 \quad \text{if } i \neq j. \qquad (83) $$

According to equations (82) and (83), the EOFs can be obtained from a diagonalization of the covariance $(\gamma_{\sigma}(\mathbf{x}_{ij}, \mathbf{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 ($Md \times Md$). It is easier to consider the dual problem to estimate the EOFs. As such, the temporal covariance is first computed:

$$ \tilde{\gamma}_{\sigma}(t_p, t_q) \triangleq \frac{1}{M} \sum_{ij} \frac{1}{\Delta T} \Delta\mathbf{X}_{ij}^p \cdot \Delta\mathbf{X}_{ij}^q. \qquad (84) $$

After diagonalizing this $n_O \times n_O$ matrix, the EOFs are obtained by projecting the eigenvectors onto the normalized increments $\frac{1}{\sqrt{\Delta T}} \Delta\mathbf{X}_{ij}^k$.
with the estimator

$$ \hat{\mathbb{E}}\{\widetilde{\Delta \mathbf{X}}_{ij}^m\} \triangleq \frac{1}{n_O} \sum_{m=1}^{n_O} \widetilde{\Delta \mathbf{X}}_{ij}^m, \qquad (79) $$

and projected onto the space of divergence-free functions:

$$ \Delta \mathbf{X}_{ij}^m \triangleq \mathcal{P}\{\Delta \mathbf{X}_{ij}^{m'}\} \quad \text{with} \quad \mathcal{P} = \mathbb{I}_d - \Delta^{-1} \nabla \nabla^T. \qquad (80) $$

Note that the projection can also be applied directly on the final EOFs $\xi_k$ – at the end of the estimation process – rather than on each increments $\Delta \mathbf{X}_{ij}^m$. Besides, if $\nabla \cdot v' = 0$ then the centred flow increments are already divergence-free 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$):

$$ \gamma_{\sigma}(\mathbf{x}_{ij}, \mathbf{x}_{pq}) \approx \frac{1}{(n_O - 1)\Delta T} \sum_{m=1}^{n_O} \Delta \mathbf{X}_{ij}^m (\Delta \mathbf{X}_{pq}^m)^T. \qquad (81) $$

Besides, we have:

$$ \gamma_{\sigma}(\mathbf{x}_{ij}, \mathbf{x}_{pq}) = \sum_{k=1}^{n_{KL}} \xi_k(\mathbf{x}_{ij}) \xi_k^T(\mathbf{x}_{pq}). \qquad (82) $$

In order to properly define the EOFs, we must add the following orthogonal constraint:

$$ \sum_{pq} \xi_i(\mathbf{x}_{pq}) \cdot \xi_j(\mathbf{x}_{pq}) = 0 \quad \text{if } i \neq j. \qquad (83) $$

According to equations (82) and (83), the EOFs can be obtained from a diagonalization of the covariance $(\gamma_{\sigma}(\mathbf{x}_{ij}, \mathbf{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 ($Md \times Md$). It is easier to consider the dual problem to estimate the EOFs. As such, the temporal covariance is first computed:

$$ \tilde{\gamma}_{\sigma}(t_p, t_q) \triangleq \frac{1}{M} \sum_{ij} \frac{1}{\Delta T} \Delta \mathbf{X}_{ij}^p \cdot \Delta \mathbf{X}_{ij}^q. \qquad (84) $$

After diagonalizing this $n_O \times n_O$ matrix, the EOFs are obtained by projecting the eigenvectors onto the normalized increments $\frac{1}{\sqrt{\Delta T}} \Delta \mathbf{X}_{ij}^k$.
Finally, after estimating the $(\xi_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 Chandramouli et al. (2019).

### 6.3.7 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, $(\xi_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 $\Delta\overline{\mathbf{X}}_{ij}^m$.

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 $\overline{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”, $\mathbf{X}_{ij}(t_0, t)$, defined by the high-resolution velocity $v$:

$$ \frac{d\mathbf{X}_{ij}(t)}{dt} = v(t, \mathbf{X}_{ij}(t_0, t)) \quad \text{and} \quad \mathbf{X}_{ij}(t_0, t_0) = \mathbf{x}_{ij}. \qquad (85) $$

- a “low-resolution flow”, $\overline{\mathbf{X}}_{ij}(t_0, t)$, defined by the low-resolution velocity $\overline{v}$:

$$ \frac{d\overline{\mathbf{X}}_{ij}(t_0, t)}{dt} = \overline{v}(t, \overline{\mathbf{X}}_{ij}(t_0, t)) \quad \text{and} \quad \overline{\mathbf{X}}_{ij}(t_0, t_0) = \mathbf{x}_{ij}. \qquad (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):

$$ \Delta\overline{\mathbf{X}}_{ij}^m \triangleq \mathbf{X}_{ij}(t_m, t_m + \Delta T) - \overline{\mathbf{X}}_{ij}(t_m, t_m + \Delta T). \qquad (87) $$

Then, the following steps are the same than in 6.3.6 (preprocessing and covariance diagonalization).
Finally, after estimating the $(\xi_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 Chandramouli et al. (2019).

### 6.3.7 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, $(\xi_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 $\Delta\overline{\mathbf{X}}_{ij}^m$.

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 $\overline{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”, $\mathbf{X}_{ij}(t_0, t)$, defined by the high-resolution velocity $v$:

$$ \frac{d\mathbf{X}_{ij}(t)}{dt} = v(t, \mathbf{X}_{ij}(t_0, t)) \quad \text{and} \quad \mathbf{X}_{ij}(t_0, t_0) = \mathbf{x}_{ij}. \qquad (85) $$

- a “low-resolution flow”, $\overline{\mathbf{X}}_{ij}(t_0, t)$, defined by the low-resolution velocity $\overline{v}$:

$$ \frac{d\overline{\mathbf{X}}_{ij}(t_0, t)}{dt} = \overline{v}(t, \overline{\mathbf{X}}_{ij}(t_0, t)) \quad \text{and} \quad \overline{\mathbf{X}}_{ij}(t_0, t_0) = \mathbf{x}_{ij}. \qquad (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):

$$ \Delta\overline{\mathbf{X}}_{ij}^m \triangleq \mathbf{X}_{ij}(t_m, t_m + \Delta T) - \overline{\mathbf{X}}_{ij}(t_m, t_m + \Delta T). \qquad (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 $\Delta t$ being small, $\Delta\tilde{X}_{ij}^m$ is close to $(v - \bar{v}) (t_m, x_{ij}) \Delta T$. Therefore, the main difference between the two methods is probably in the definition of the “residual velocities”: $v'$ on one hand and $(v - \bar{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 (Resseguier et al., 2019). Resseguier et al. (2019) 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.

### 6.3.8 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. Holmes et al., 1996). By definition, the resolved velocity component lies on this subspace, whereas the unresolved component does not (Mémin, 2014; Resseguier et al., 2015, 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 \mapsto \sigma dB_t \cdot \nabla f, \qquad (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 \times 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_r^2 \times n_r^2$ coefficients.

In comparison, the approximate unresolved velocity Karhunen-Loeve decomposition would lead to stochastic operator

$$f \mapsto \sum_{k=1}^{n_{KL}} (\xi_k \cdot \nabla f) \dot{W}_k, \qquad (89)$$

characterized by $n_r^2 \times n_{KL}$ coefficients. For $n_{KL} < n_r^2$, 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_r^2$). 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).
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 $\Delta t$ being small, $\Delta\tilde{X}_{ij}^m$ is close to $(v - \bar{v})(t_m, x_{ij}) \Delta T$. Therefore, the main difference between the two methods is probably in the definition of the “residual velocities”: $v'$ on one hand and $(v - \bar{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 (Resseguier et al., 2019). Resseguier et al. (2019) 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.

### 6.3.8 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. Holmes et al., 1996). By definition, the resolved velocity component lies on this subspace, whereas the unresolved component does not (Mémin, 2014; Resseguier et al., 2015, 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 \mapsto \sigma dB_t \cdot \nabla f, \qquad (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 \times 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_r^2 \times n_r^2$ coefficients.

In comparison, the approximate unresolved velocity Karhunen-Loeve decomposition would lead to stochastic operator

$$f \mapsto \sum_{k=1}^{n_{KL}} (\xi_k \cdot \nabla f) \dot{W}_k, \qquad (89)$$

characterized by $n_r^2 \times n_{KL}$ coefficients. For $n_{KL} < n_r^2$, 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_r^2$). 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).