Nevertheless, for large state-space dimension $n_r$, this $n_r^2 \times n_r^2$ covariance matrix of the random operator is too large to handle. Hence, additional assumptions on the unresolved small-scale velocity (e.g. homogeneity, parametric model, transport) or a Karhunen-Loeve decomposition become necessary.

### 6.3.9 Transport of the unresolved velocity

Methods 6.3.6 and 6.3.7 learn the unresolved small-scale structure from high-resolution simulations. However, such simulations are not always available. Moreover, methods 6.3.6 and 6.3.7 assume heterogeneity without time dependence. While some turbulence heterogeneities could be fixed in time (e.g. heterogeneity due to boundaries conditions), we can expect that a large part of it moves with the flow.

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

$$D_t^S \xi_k = (\xi_k \cdot \nabla)(w^* dt + \sigma \circ dB_t). \quad (90)$$

As such, the unresolved velocity is both heterogeneous and non-stationary in time, without requiring any learning. However, as far as we know, this method has not been tested in numerical simulations yet.

### 6.3.10 On-line learning of EOF from velocity

We propose here a new approach – based on the EOF decomposition (75) – but where the EOFs are time-dependent. They are estimated on-line from a coarse simulation. For this purpose, we propose to generate pseudo-observations of the small scales directly from the resolved velocity. Then, we will compute an EOF representation as in 6.3.6. Finally, we will adapt the noise variance to smaller scales using a turbulence power-law scaling (originally proposed by Harouna and Mémin, 2017).

#### Pseudo-observations

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

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

$$I(x_{ij}, t) \triangleq \left\{ w(x_{pq}, t) \mid |p - i| \le \frac{n-1}{2}, |q - j| \le \frac{n-1}{2} \right\}. \quad (91)$$
Nevertheless, for large state-space dimension $n_r$, this $n_r^2 \times n_r^2$ covariance matrix of the random operator is too large to handle. Hence, additional assumptions on the unresolved small-scale velocity (e.g. homogeneity, parametric model, transport) or a Karhunen-Loeve decomposition become necessary.

### 6.3.9 Transport of the unresolved velocity

Methods 6.3.6 and 6.3.7 learn the unresolved small-scale structure from high-resolution simulations. However, such simulations are not always available. Moreover, methods 6.3.6 and 6.3.7 assume heterogeneity without time dependence. While some turbulence heterogeneities could be fixed in time (e.g. heterogeneity due to boundaries conditions), we can expect that a large part of it moves with the flow.

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

$$D_t^S \xi_k = (\xi_k \cdot \nabla)(w^* dt + \sigma \circ dB_t). \quad (90)$$

As such, the unresolved velocity is both heterogeneous and non-stationary in time, without requiring any learning. However, as far as we know, this method has not been tested in numerical simulations yet.

### 6.3.10 On-line learning of EOF from velocity

We propose here a new approach – based on the EOF decomposition (75) – but where the EOFs are time-dependent. They are estimated on-line from a coarse simulation. For this purpose, we propose to generate pseudo-observations of the small scales directly from the resolved velocity. Then, we will compute an EOF representation as in 6.3.6. Finally, we will adapt the noise variance to smaller scales using a turbulence power-law scaling (originally proposed by Harouna and Mémin, 2017).

#### Pseudo-observations

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

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

$$I(x_{ij}, t) \triangleq \left\{ w(x_{pq}, t) \mid |p - i| \le \frac{n-1}{2}, |q - j| \le \frac{n-1}{2} \right\}. \quad (91)$$
Then, for each $m \in \{1, ..., n_O\}$, for each point $x_{ij}$ independently, we set the value of the pseudo-observation $v'(x_{ij}, t, m)$ by randomly choosing a value in the set $I(x_{ij}, t)$.

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

## Rescaling

From the SVD, we obtain a set of EOFs $\xi_k^{(L)}$ and hence a model for the unresolved velocity $\sigma^{(L)}(x, t)\dot{B}$:

$$ \sigma^{(L)} \dot{B} = \sum_{k=1}^{n_{KL}} \xi_k^{(L)} \dot{W}_k. \qquad (92) $$

The pseudo-realizations $v' = v'_L$ have been generated at a spatial scale $L = nl$. These fluctuations correspond to a virtual observation scale $L$ ($n \times n$ window) and must be scaled down to the “simulation scale” $l$. Therefore the unresolved velocity variance tensor, $a$, is rescaled by a coefficient proposed by Harouna and Mémin (2017), and which in 2D reads:

$$ a_l = \left( \frac{l}{L} \right)^{2/3} a_L, \qquad (93) $$

where $a_L$ and $a_l$ are the variance tensors at the scales $L$ and $l$ respectively. This scaling relies on the Kolmogorov-Richardson cascade assumption with the velocity fluctuations at scale $l$ as $u_l \propto \epsilon^{1/3}l^{1/3}$. The unresolved velocity, $\sigma\dot{B}$, can be finally simulated – at the “simulation scale” $l$ – as:

$$ \sigma \dot{B} = n^{-1/3} \sigma^{(L)} \dot{B}. \qquad (94) $$

Let us note that in 3D the scaling must be adapted as the power exponent is 4/3 Harouna and Mémin (2017).

The on-line SVD method 6.3.10 will be numerically tested in section 8. There, the LU parameterizations 6.3.2, 6.3.3 and 6.3.10 will be compared to the new scheme 5.3.2 and to ensembles generated by initial condition perturbations.

## 7 Metrics for UQ

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

This section is a short summary about existing verification ensemble tools. More detailed review could be found in (Hamill and Colucci, 1997; Gneiting and Raftery, 2007; Scheuerer and Hamill, 2015; Thorarinsdottir, 2017).
Then, for each $m \in \{1, ..., n_O\}$, for each point $x_{ij}$ independently, we set the value of the pseudo-observation $v'(x_{ij}, t, m)$ by randomly choosing a value in the set $I(x_{ij}, t)$.

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

### Rescaling

From the SVD, we obtain a set of EOFs $\xi_k^{(L)}$ and hence a model for the unresolved velocity $\sigma^{(L)}(x, t)\dot{B}$:

$$ \sigma^{(L)} \dot{B} = \sum_{k=1}^{n_{KL}} \xi_k^{(L)} \dot{W}_k. \qquad (92) $$

The pseudo-realizations $v' = v'_L$ have been generated at a spatial scale $L = nl$. These fluctuations correspond to a virtual observation scale $L$ ($n \times n$ window) and must be scaled down to the “simulation scale” $l$. Therefore the unresolved velocity variance tensor, $a$, is rescaled by a coefficient proposed by Harouna and Mémin (2017), and which in 2D reads:

$$ a_l = \left( \frac{l}{L} \right)^{2/3} a_L, \qquad (93) $$

where $a_L$ and $a_l$ are the variance tensors at the scales $L$ and $l$ respectively. This scaling relies on the Kolmogorov-Richardson cascade assumption with the velocity fluctuations at scale $l$ as $u_l \propto \epsilon^{1/3}l^{1/3}$. The unresolved velocity, $\sigma\dot{B}$, can be finally simulated – at the “simulation scale” $l$ – as:

$$ \sigma \dot{B} = n^{-1/3} \sigma^{(L)} \dot{B}. \qquad (94) $$

Let us note that in 3D the scaling must be adapted as the power exponent is 4/3 Harouna and Mémin (2017).

The on-line SVD method 6.3.10 will be numerically tested in section 8. There, the LU parameterizations 6.3.2, 6.3.3 and 6.3.10 will be compared to the new scheme 5.3.2 and to ensembles generated by initial condition perturbations.

## 7 Metrics for UQ

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

This section is a short summary about existing verification ensemble tools. More detailed review could be found in (Hamill and Colucci, 1997; Gneiting and Raftery, 2007; Scheuerer and Hamill, 2015; Thorarinsdottir, 2017).
## 7.1 Talagrand diagram

A Talagrand diagram (or rank histogram) is a technique used to check the reliability of an ensemble forecast or a set of quantile. The idea is considered the availability of N observation $q_1^o, ..., q_N^o$ in a set of rank r. The sets of rank are built from rank statistic which is a sorted ensemble $\{(q_1^{(i_1)}, ..., q_1^{(i_{n_e})}), ..., (q_N^{(i_1)}, ..., q_N^{(i_{n_e})})\}$ of $n_e$ members (Anderson, 1996; Talagrand, 1999; Hamill and Colucci, 1997):

$$ \forall (j, k) \in \{1, ..., n_e\} \times \{1, ..., N\}, r_j = \frac{1}{N} \sum_{k=1}^{N} \hat{P}(q_k^{(i_{j-1})} \le q_k^o < q_k^{(i_j)}), \quad (95) $$

where $\hat{P}$ is the estimated probability of ranking an observation between two sorted ensemble members. In the context of (95), N represents the spatial points or times steps of the observation.

A calibrated ensemble should result in a flat histogram. However, a flat histogram does not guarantee a calibrated ensemble (Hamill, 2001). A flat histogram mostly indicates that the ensemble and observation are sampled from a common distribution. A U-shaped rank histogram corresponds to an underdispersion or to a conditional bias (Hamill, 2001) of the ensemble members. A dome-shaped shows an overdispersion of the ensemble. A nonsymmetric histogram is the footprint of bias.

## 7.2 Mean absolute and mean squared error

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

$$ \forall (j, k) \in \{1, ..., n_e\} \times \{1, ..., N\}, \text{MAE} = \frac{1}{Nn_e} \sum_{k=1}^{N} \sum_{j=1}^{n_e} |q_k^o - q_k^{(j)}|, \quad (96) $$

and the MSE by this one:

$$ \forall (j, k) \in \{1, ..., n_e\} \times \{1, ..., N\}, \text{MSE} = \frac{1}{Nn_e} \sum_{k=1}^{N} \sum_{j=1}^{n_e} (q_k^o - q_k^{(j)})^2. \quad (97) $$

The goal behind these two equations is to calculate the total error between a sample of observation and ensemble of size N. So, the empirical mean of all error between an observation $q_k^o$ and each member $q_k^{(j)}$ of the ensemble of size $n_e$ is calculated. Both metrics are proper scores. A proper score is negatively oriented, such that a lower score indicated a better ensemble forecast. A proper score converging to zero means that the observation cannot be distinguished from the ensemble members. There are several proper scoring rules.
## 7.1 Talagrand diagram

A Talagrand diagram (or rank histogram) is a technique used to check the reliability of an ensemble forecast or a set of quantile. The idea is considered the availability of N observation $q_1^o, ..., q_N^o$ in a set of rank r. The sets of rank are built from rank statistic which is a sorted ensemble $\{(q_1^{(i_1)}, ..., q_1^{(i_{n_e})}), ..., (q_N^{(i_1)}, ..., q_N^{(i_{n_e})})\}$ of $n_e$ members (Anderson, 1996; Talagrand, 1999; Hamill and Colucci, 1997):

$$ \forall (j, k) \in \{1, ..., n_e\} \times \{1, ..., N\}, r_j = \frac{1}{N} \sum_{k=1}^{N} \hat{P}(q_k^{(i_{j-1})} \le q_k^o < q_k^{(i_j)}), \quad (95) $$

where $\hat{P}$ is the estimated probability of ranking an observation between two sorted ensemble members. In the context of (95), N represents the spatial points or times steps of the observation.

A calibrated ensemble should result in a flat histogram. However, a flat histogram does not guarantee a calibrated ensemble (Hamill, 2001). A flat histogram mostly indicates that the ensemble and observation are sampled from a common distribution. A U-shaped rank histogram corresponds to an underdispersion or to a conditional bias (Hamill, 2001) of the ensemble members. A dome-shaped shows an overdispersion of the ensemble. A nonsymmetric histogram is the footprint of bias.

## 7.2 Mean absolute and mean squared error

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

$$ \forall (j, k) \in \{1, ..., n_e\} \times \{1, ..., N\}, \text{MAE} = \frac{1}{Nn_e} \sum_{k=1}^{N} \sum_{j=1}^{n_e} |q_k^o - q_k^{(j)}|, \quad (96) $$

and the MSE by this one:

$$ \forall (j, k) \in \{1, ..., n_e\} \times \{1, ..., N\}, \text{MSE} = \frac{1}{Nn_e} \sum_{k=1}^{N} \sum_{j=1}^{n_e} (q_k^o - q_k^{(j)})^2. \quad (97) $$

The goal behind these two equations is to calculate the total error between a sample of observation and ensemble of size N. So, the empirical mean of all error between an observation $q_k^o$ and each member $q_k^{(j)}$ of the ensemble of size $n_e$ is calculated. Both metrics are proper scores. A proper score is negatively oriented, such that a lower score indicated a better ensemble forecast. A proper score converging to zero means that the observation cannot be distinguished from the ensemble members. There are several proper scoring rules.
## 7.3 Continuous ranked probability score

Another proper score used for ensemble model evaluation is the continuous ranked probability score (CRPS). This kind of score is characterized by the so-called predictive distribution $F$ – the distribution represented by the ensemble forecast {$q^{(1)}, ..., q^{(n_e)}$} realizations and the observation $q^o$. The CRPS is defined by the following equation (Matheson and Winkler, 1976; Hersbach, 2000):

$$ \mathrm{CRPS}(F, q^o) = \int_{\mathbb{R}} (F(\mathcal{B}) - H(\mathcal{B} - q^o))^2 d\mathcal{B}, \qquad (98) $$

where $H$ is the heaviside function (equals to one if $q^o \le \mathcal{B}$ and zero otherwise), $F$ is the cumulative distribution function (CDF) of the random process $q$ (represented by the ensemble) at **one spatial point $x$ and time $t$**. The CRPS is the distance between the heaviside function and the distribution of the random variable $q$. This heaviside function represents the inequality between an ensemble member and the observation. If $F$ has a finite first-order moment, there is another representation of the CRPS (Gneiting and Raftery, 2007):

$$ \mathrm{CRPS}(F, q^o) = \mathbb{E}_F |q - q^o| - \frac{1}{2} \mathbb{E}_F |q - q'|, \qquad (99) $$

where $q$ and $q'$ are independent realizations of $F$. If a priori assumption are made on the behaviour of the random variable $q$. It is possible to deduce an analytic expression of the CRPS thanks to the equation (98) and (99).

In the case where $F$ is issued from a normal distribution $N(\mu, \sigma^2)$, the CRPS formula is (Gneiting et al., 2005):

$$ \mathrm{CRPS}(F, q^o) = \sigma \left( \omega(2\Phi(\omega) - 1) + 2\phi(\omega) - \frac{1}{\sqrt{\pi}} \right), \qquad (100) $$

where $\omega = \frac{(q^o - \mu)}{\sigma}$, $\phi$ is probability density function (PDF) of the standard normal distribution and $\Phi$ is the CDF of the standard normal distribution.

Also following (Grimit et al., 2006), for an ensemble forecast {$q^{(1)}, ..., q^{(n_e)}$}, a natural fair estimator of the CRPS is given by:

$$ \widehat{\mathrm{CRPS}}(F, q^o) = \frac{1}{n_e} \sum_{i=1}^{n_e} |q^{(i)} - q^o| - \frac{1}{2n_e^2} \sum_{i=1}^{n_e} \sum_{j=1}^{n_e} |q^{(i)} - q^{(j)}|. \qquad (101) $$

The expression (101) allows to compute the CRPS without making a priori assumptions on our random variable. However, this expression has a $O(n_e^2)$ computational complexity.

An algebraically equivalent representation of the CRPS based on the generalized quantile function is proposed by (Laio and Tamea, 2007)

$$ \widehat{\mathrm{CRPS}}(\hat{F}, q^o) = \frac{2}{n_e^2} \sum_{j=1}^{n_e} \left( (q^{(i_j)} - q^o)(n_e H(q^{(i_j)} - q^o) - j + \frac{1}{2}) \right), \qquad (102) $$
## 7.3 Continuous ranked probability score

Another proper score used for ensemble model evaluation is the continuous ranked probability score (CRPS). This kind of score is characterized by the so-called predictive distribution $F$ – the distribution represented by the ensemble forecast {$q^{(1)}, ..., q^{(n_e)}$} realizations and the observation $q^o$. The CRPS is defined by the following equation (Matheson and Winkler, 1976; Hersbach, 2000):

$$ \mathrm{CRPS}(F, q^o) = \int_{\mathbb{R}} (F(\mathcal{B}) - H(\mathcal{B} - q^o))^2 d\mathcal{B}, \qquad (98) $$

where $H$ is the heaviside function (equals to one if $q^o \le \mathcal{B}$ and zero otherwise), $F$ is the cumulative distribution function (CDF) of the random process $q$ (represented by the ensemble) at one spatial point $x$ and time $t$. The CRPS is the distance between the heaviside function and the distribution of the random variable $q$. This heaviside function represents the inequality between an ensemble member and the observation. If $F$ has a finite first-order moment, there is another representation of the CRPS (Gneiting and Raftery, 2007):

$$ \mathrm{CRPS}(F, q^o) = \mathbb{E}_F |q - q^o| - \frac{1}{2} \mathbb{E}_F |q - q'|, \qquad (99) $$

where $q$ and $q'$ are independent realizations of $F$. If a priori assumption are made on the behaviour of the random variable $q$. It is possible to deduce an analytic expression of the CRPS thanks to the equation (98) and (99).

In the case where $F$ is issued from a normal distribution $N(\mu, \sigma^2)$, the CRPS formula is (Gneiting et al., 2005):

$$ \mathrm{CRPS}(F, q^o) = \sigma \left( \omega(2\Phi(\omega) - 1) + 2\phi(\omega) - \frac{1}{\sqrt{\pi}} \right), \qquad (100) $$

where $\omega = \frac{(q^o - \mu)}{\sigma}$, $\phi$ is probability density function (PDF) of the standard normal distribution and $\Phi$ is the CDF of the standard normal distribution.

Also following (Grimit et al., 2006), for an ensemble forecast {$q^{(1)}, ..., q^{(n_e)}$}, a natural fair estimator of the CRPS is given by:

$$ \widehat{\mathrm{CRPS}}(F, q^o) = \frac{1}{n_e} \sum_{i=1}^{n_e} |q^{(i)} - q^o| - \frac{1}{2n_e^2} \sum_{i=1}^{n_e} \sum_{j=1}^{n_e} |q^{(i)} - q^{(j)}|. \qquad (101) $$

The expression (101) allows to compute the CRPS without making a priori assumptions on our random variable. However, this expression has a $O(n_e^2)$ computational complexity.

An algebraically equivalent representation of the CRPS based on the generalized quantile function is proposed by (Laio and Tamea, 2007)

$$ \widehat{\mathrm{CRPS}}(\hat{F}, q^o) = \frac{2}{n_e^2} \sum_{j=1}^{n_e} \left( (q^{(i_j)} - q^o)(n_e H(q^{(i_j)} - q^o) - j + \frac{1}{2}) \right), \qquad (102) $$
where {$q^{(i_1)}, ..., q^{(i_{n_e})}$} is the sorted ensemble forecast. The computational complexity of the expression (102) is $O(n_e \log(n_e))$. The interested reader can refer to (Jordan, 2016) for mathematical details on the discrete form of the CRPS.

## 7.4 Energy score

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

There, our observation $q^o$ and our random variable $q$ is multivariate. The corresponding ensemble {$((q_1^{(1)}, ..., q_N^{(1)}), ..., (q_1^{(n_e)}, ..., q_N^{(n_e)}))$} contains $n_e$ realizations. Each of them represents the field of interest in $N$ time steps or spatial points. $F$ is now the multivariate CDF associated to $q$. We need a generalized form of the CRPS to work in this case. Following (Gneiting and Raftery, 2007; Scheuerer and Hamill, 2015), a possible generalization of the CRPS is the energy score (ES):

$$ \mathrm{ES}(F, q^o) = \mathbb{E}_F \|q - q^o\| - \frac{1}{2} \mathbb{E}_F \|q - q'\|, \quad (103) $$

where $\|\cdot\|$ is the euclidean norm. So, an estimator of the ES would be:

$$ \widehat{\mathrm{ES}}(\hat{F}, q^o) = \frac{1}{n_e} \sum_{i=1}^{n_e} \|q^{(i)} - q^o\| - \frac{1}{2n_e^2} \sum_{i=1}^{n_e} \sum_{j=1}^{n_e} \|q^{(i)} - q^{(j)}\| \quad (104) $$

## 7.5 Variogram score - p

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

One tool used to study these correlations is the variogram of order $p$ (also called structure function). This tool consists in studying the pairwise differences of the component of the multivariate ensemble forecast and observation by considering the statistics of our multivariate observation. Following (Scheuerer and Hamill, 2015), the p-variogram score (VS-p) is an alternative
where {$q^{(i_1)}, ..., q^{(i_{n_e})}$} is the sorted ensemble forecast. The computational complexity of the expression (102) is $O(n_{e}log(n_{e}))$. The interested reader can refer to (Jordan, 2016) for mathematical details on the discrete form of the CRPS.

## 7.4 Energy score

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

There, our observation $q^o$ and our random variable $q$ is multivariate. The corresponding ensemble {$((q_1^{(1)}, ..., q_N^{(1)}), ..., (q_1^{(n_e)}, ..., q_N^{(n_e)}))$} contains $n_e$ realizations. Each of them represents the field of interest in $N$ time steps or spatial points. $F$ is now the multivariate CDF associated to $q$. We need a generalized form of the CRPS to work in this case. Following (Gneiting and Raftery, 2007; Scheuerer and Hamill, 2015), a possible generalization of the CRPS is the energy score (ES):

$$ES(F, q^o) = E_F \|q - q^o\| - \frac{1}{2} E_F \|q - q'\|, \quad (103)$$

where $\|\cdot\|$ is the euclidean norm. So, an estimator of the ES would be:

$$\widehat{ES}(\hat{F}, q^o) = \frac{1}{n_e} \sum_{i=1}^{n_e} \|q^{(i)} - q^o\| - \frac{1}{2n_e^2} \sum_{i=1}^{n_e} \sum_{j=1}^{n_e} \|q^{(i)} - q^{(j)}\| \quad (104)$$

## 7.5 Variogram score - p

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

One tool used to study these correlations is the variogram of order $p$ (also called structure function). This tool consists in studying the pairwise differences of the component of the multivariate ensemble forecast and observation by considering the statistics of our multivariate observation. Following (Scheuerer and Hamill, 2015), the p-variogram score (VS-p) is an alternative
solution which overcomes this drawback. We express it using the same multi-variate random variable $q$ and observation $q^o$:

$$VS^p(F, q^o) = \sum_{i=1}^{N} \sum_{j=1}^{N} \frac{1}{w_{ij}} (|q_i^o - q_j^o|^p - \mathbb{E}_F |q_i - q_j|^p)^2, \quad (105)$$

where the $w_{ij}$ are non negative weights that allows one to emphasize or down-weight pairs of component combinations based on subjective expert decisions. $p$ is the order of the variogram score. $N$ represents the number of spatial points or time steps. $F$ is the multivariate cumulative distribution function of $q$. According to (Scheuerer and Hamill, 2015), a typical choice of weight to discriminate time or space dependencies is the inverse distance weight: $w_{ij} := \frac{1}{\sqrt{|i-j|}}$ and $p := 0.5$. The VS-p score measures the dissimilarity between approximations of the variograms of order $p$ of the observations and forecast over all pairs of the components of our variable of interest.

## 8 Numerical results : Ensemble forecasts verification for SQG dynamics

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

### 8.1 Test flows

#### 8.1.1 Deterministic and stochastic surface quasi-geostrophic (SQG) dynamics

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

$$b \triangleq -g \frac{\rho'}{\rho_0} \quad \text{with} \quad \rho(x, y, z, t) = \rho_0 \left(1 - \frac{N^2}{g} z\right) + \rho'(x, y, z, t), \quad (106)$$

where $\rho$ is the density and $N$ the stratification. In the ocean, the density anomaly is small compared to the total density (Boussinesq approximation) and the flow is approximately isochoric. The conservation of salinity and temperature, with a linearized equation of state, provides the transport of buoyancy (Vallis, 2006). Then, considering rapid rotation, strong stratification and uniform potential vorticity leads to the so-called SQG model (Blumen, 1978; Held et al., 1995; Lapeyre, 2017).
solution which overcomes this drawback. We express it using the same multi-variate random variable $q$ and observation $q^o$:

$$VS^p(F, q^o) = \sum_{i=1}^{N} \sum_{j=1}^{N} \frac{1}{w_{ij}} (|q_i^o - q_j^o|^p - \mathbb{E}_F |q_i - q_j|^p)^2, \quad (105)$$

where the $w_{ij}$ are non negative weights that allows one to emphasize or down-weight pairs of component combinations based on subjective expert decisions. $p$ is the order of the variogram score. $N$ represents the number of spatial points or time steps. $F$ is the multivariate cumulative distribution function of $q$. According to (Scheuerer and Hamill, 2015), a typical choice of weight to discriminate time or space dependencies is the inverse distance weight: $w_{ij} := \frac{1}{\sqrt{|i-j|}}$ and $p := 0.5$. The VS-p score measures the dissimilarity between approximations of the variograms of order $p$ of the observations and forecast over all pairs of the components of our variable of interest.

## 8 Numerical results : Ensemble forecasts verification for SQG dynamics

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

### 8.1 Test flows

#### 8.1.1 Deterministic and stochastic surface quasi-geostrophic (SQG) dynamics

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

$$b \triangleq -g \frac{\rho'}{\rho_0} \quad \text{with} \quad \rho(x, y, z, t) = \rho_0 \left(1 - \frac{N^2}{g} z\right) + \rho'(x, y, z, t), \quad (106)$$

where $\rho$ is the density and $N$ the stratification. In the ocean, the density anomaly is small compared to the total density (Boussinesq approximation) and the flow is approximately isochoric. The conservation of salinity and temperature, with a linearized equation of state, provides the transport of buoyancy (Vallis, 2006). Then, considering rapid rotation, strong stratification and uniform potential vorticity leads to the so-called SQG model (Blumen, 1978; Held et al., 1995; Lapeyre, 2017).
A random version of this model (denoted SQG<sub>MU</sub>) can be derived from the location uncertainty principle (Resseguier et al., 2017b). It will keep the same structure except that the buoyancy is now transported in the stochastic sense (14). The horizontal velocity **u** = **w** is related to the buoyancy **b** in Fourier space through the usual SQG relation:

$$ \hat{\mathbf{u}} = i\mathbf{k}^{\perp} \frac{\hat{b}}{N \|\mathbf{k}\|}, \qquad (107) $$

where **k** is the horizontal wave-vector. The unresolved velocity, $\sigma \dot{\mathbf{B}}$, is also horizontal. Consequently, the variance tensor, **a**, is a 2 × 2 matrix and the SQG<sub>MU</sub> model is two-dimensional.

### 8.1.2 Our SQG simulations

High-resolution deterministic SQG simulations of test flows will provide references to which we will compare random simulations performed at a lower resolution. For this purpose, the high-resolution (512²) simulation outputs will be projected onto the space of low-resolution (128²) fields, i.e. it will be adequately filtered and subsampled to a low-resolution. We will refer to this projected fields as "observations" and denote them $b^o$. For all simulations, a standard hyperviscosity (HV) scheme has been introduced (Held et al., 1995):

$$ D_t b = \alpha \Delta^4 b \, dt, \qquad (108) $$

with a positive coefficient $\alpha$ proportional to $M_x^{-8}$ where $M_x$ denotes the grid size (i.e. 128 or 512). The domain size is a square box $L_x \times L_y = 1000 \, \text{km} \times 1000 \, \text{km}$ and the boundary conditions are doubly periodic.

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

$$ b_0(\mathbf{x}) \triangleq F\left(\mathbf{x} - \begin{pmatrix} 250 \, \text{km} \\ 250 \, \text{km} \end{pmatrix}\right) + F\left(\mathbf{x} - \begin{pmatrix} 750 \, \text{km} \\ 250 \, \text{km} \end{pmatrix}\right) \\ -F\left(\mathbf{x} - \begin{pmatrix} 250 \, \text{km} \\ 750 \, \text{km} \end{pmatrix}\right) - F\left(\mathbf{x} - \begin{pmatrix} 750 \, \text{km} \\ 750 \, \text{km} \end{pmatrix}\right), \qquad (109) $$

with

$$ F(\mathbf{x}) \triangleq B_0 \exp \left( -\frac{1}{2} \left( \frac{x^2}{\sigma_x^2} + \frac{y^2}{\sigma_y^2} \right) \right) \quad \text{and} \quad \begin{cases} \sigma_x = 67 \, \text{km}, \\ \sigma_y = 133 \, \text{km}. \end{cases} \qquad (110) $$

The amplitude of the buoyancy and the stratification are set to: $B_0 = 10^{-3} \, \text{m.s}^{-2}$ and $N = 3f_0$. The Coriolis frequency $f_0$ is fixed to $1.028 \times 10^{-4} \, \text{s}^{-1}$, which corresponds to a latitude of 45°.
A random version of this model (denoted SQG<sub>MU</sub>) can be derived from the location uncertainty principle (Resseguier et al., 2017b). It will keep the same structure except that the buoyancy is now transported in the stochastic sense (14). The horizontal velocity **u** = **w** is related to the buoyancy **b** in Fourier space through the usual SQG relation:

$$ \hat{\mathbf{u}} = i\mathbf{k}^{\perp} \frac{\hat{b}}{N\|\mathbf{k}\|}, \qquad (107) $$

where **k** is the horizontal wave-vector. The unresolved velocity, $\sigma \dot{\mathbf{B}}$, is also horizontal. Consequently, the variance tensor, **a**, is a 2 × 2 matrix and the SQG<sub>MU</sub> model is two-dimensional.

### 8.1.2 Our SQG simulations

High-resolution deterministic SQG simulations of test flows will provide references to which we will compare random simulations performed at a lower resolution. For this purpose, the high-resolution (512<sup>2</sup>) simulation outputs will be projected onto the space of low-resolution (128<sup>2</sup>) fields, i.e. it will be adequately filtered and subsampled to a low-resolution. We will refer to this projected fields as "observations" and denote them $b^o$. For all simulations, a standard hyperviscosity (HV) scheme has been introduced (Held et al., 1995):

$$ D_t b = \alpha \Delta^4 b \, dt, \qquad (108) $$

with a positive coefficient $\alpha$ proportional to $M_x^{-8}$ where $M_x$ denotes the grid size (i.e. 128 or 512). The domain size is a square box $L_x \times L_y = 1000 \, \text{km} \times 1000 \, \text{km}$ and the boundary conditions are doubly periodic.

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

$$ b_0(\mathbf{x}) \triangleq F\left(\mathbf{x} - \begin{pmatrix} 250 \, \text{km} \\ 250 \, \text{km} \end{pmatrix}\right) + F\left(\mathbf{x} - \begin{pmatrix} 750 \, \text{km} \\ 250 \, \text{km} \end{pmatrix}\right) \\ -F\left(\mathbf{x} - \begin{pmatrix} 250 \, \text{km} \\ 750 \, \text{km} \end{pmatrix}\right) - F\left(\mathbf{x} - \begin{pmatrix} 750 \, \text{km} \\ 750 \, \text{km} \end{pmatrix}\right), \qquad (109) $$

with

$$ F(\mathbf{x}) \triangleq B_0 \exp \left( -\frac{1}{2} \left( \frac{x^2}{\sigma_x^2} + \frac{y^2}{\sigma_y^2} \right) \right) \quad \text{and} \quad \begin{cases} \sigma_x = 67 \, \text{km}, \\ \sigma_y = 133 \, \text{km}. \end{cases} \qquad (110) $$

The amplitude of the buoyancy and the stratification are set to: $B_0 = 10^{-3} \, \text{m.s}^{-2}$ and $N = 3f_0$. The Coriolis frequency $f_0$ is fixed to $1.028 \times 10^{-4} \, \text{s}^{-1}$, which corresponds to a latitude of 45°.
Another type of initial buoyancy field is shown in figure 2(b), which is a homogeneous Gaussian random field generated from a spectrum of buoyancy with a prescribed slope equals to -5/3. This slope corresponds to the power law of a developped SQG turbulence.

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Two initial buoyancy (m.s⁻²) fields: (a) – A smooth field with four given vortices (cyclones in blue and anti-cyclones in red); (b) – A homogeneous Gaussian field generated from a -5/3-spectrum.</figcaption></figure>

We will now forecast the initial buoyancy fields of figure 2 with several ensembles of simulations and assess the performance of each of these ensembles. Short term and long term forecasts will be treated separately because they are very different in nature.

## 8.2 Short-term ensemble forecast

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

Several low-resolution ensembles have been forecast. Two ensembles rely on initial conditions random perturbations, three on dynamics under location uncertainty and five on the dissipation-adapted noise (see section 5.3.2). Specifically, the small-scale velocity in the dynamics under location uncertainty and the random perturbation of initial condition in $SQG_{PIC}$ have been both generated using three approaches: a homogeneous stationary model (denoted "Spectral"), a homogeneous non-stationary model (denoted "ADSD") and a heterogeneous non-stationary model (denoted "SVDpseudo"). For the
Another type of initial buoyancy field is shown in figure 2(b), which is a homogeneous Gaussian random field generated from a spectrum of buoyancy with a prescribed slope equals to -5/3. This slope corresponds to the power law of a developped SQG turbulence.

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Two initial buoyancy (m.s⁻²) fields: (a) – A smooth field with four given vortices (cyclones in blue and anti-cyclones in red); (b) – A homogeneous Gaussian field generated from a -5/3-spectrum.</figcaption></figure>

We will now forecast the initial buoyancy fields of figure 2 with several ensembles of simulations and assess the performance of each of these ensembles. Short term and long term forecasts will be treated separately because they are very different in nature.

## 8.2 Short-term ensemble forecast

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

Several low-resolution ensembles have been forecast. Two ensembles rely on initial conditions random perturbations, three on dynamics under location uncertainty and five on the dissipation-adapted noise (see section 5.3.2). Specifically, the small-scale velocity in the dynamics under location uncertainty and the random perturbation of initial condition in $SQG_{PIC}$ have been both generated using three approaches: a homogeneous stationary model (denoted "Spectral"), a homogeneous non-stationary model (denoted "ADSD") and a heterogeneous non-stationary model (denoted "SVDpseudo"). For the
High-resolution buoyancy

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Buoyancy ($m.s^{-2}$) at $t = 10, 15, 16, 18, 20, 30$ days of advection for the usual SQG model at resolution 512$^2$, based on the smooth initial field – figure 2(a).</figcaption></figure>

model under location uncertainty, it corresponds to parameterizations 6.3.2, 6.3.3 and 6.3.10 respectively, of section 6.3. For the “Spectral” method, the small-scale energy is specified by the diffusion coefficient $\frac{a}{2} = 9 m^2.s^{-1}$. We also need to choose the spectrum slope and the minimum and maximum wavelengths of the unresolved velocity. For more, information on the values of these parameters, the reader can refer to Resseguier et al. (2017b). In contrast, for the improved “ADSD” method, we do not need to choose any parameter. The parameterization must adapt itself using large-scale informations at each time step. Despite the difference of flow and resolution, the very same code – without any modification – is used here and in Resseguier et al. (2019) to
High-resolution buoyancy

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Buoyancy ($m.s^{-2}$) at $t = 10, 15, 16, 18, 20, 30$ days of advection for the usual SQG model at resolution 512$^2$, based on the smooth initial field – figure 2(a).</figcaption></figure>

model under location uncertainty, it corresponds to parameterizations 6.3.2, 6.3.3 and 6.3.10 respectively, of section 6.3. For the “Spectral” method, the small-scale energy is specified by the diffusion coefficient $\frac{a}{2} = 9 m^2.s^{-1}$. We also need to choose the spectrum slope and the minimum and maximum wavelengths of the unresolved velocity. For more, information on the values of these parameters, the reader can refer to Resseguier et al. (2017b). In contrast, for the improved “ADSD” method, we do not need to choose any parameter. The parameterization must adapt itself using large-scale informations at each time step. Despite the difference of flow and resolution, the very same code – without any modification – is used here and in Resseguier et al. (2019) to
High-resolution buoyancy

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Buoyancy (m.s⁻²) at t = 3, 6, 9, 12, 15, 18 days of advection for the usual SQG model at resolution 512², based on the random initial field – figure 2(b).</figcaption></figure>

define the “ADSD” small-scale velocity. For the ensembles with perturbed initial conditions, similar sampling methods are adopted. But the obtained “small-scale” random fields are used only at the initial time. As an example, figure 5 shows the initial perturbation spectrum used in the “Spectral” method for the first initial condition (109). Then, it is added to the initial condition b(x, t = 0). Besides, the dissipation-adapted noise 5.3.2 has been implemented in SQG_MUMATLAB code. For this method – denoted “WavHypervis” – we use a wavelet basis (e_k)_{k∈Z}. We have forecast 5 ensembles with 5 different values of the scaling factor ζ. For instance, figure 6 shows us that under the same low-resolution, when the scaling factor ζ increases (i.e. when the noise
High-resolution buoyancy

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Buoyancy (m.s⁻²) at t = 3, 6, 9, 12, 15, 18 days of advection for the usual SQG model at resolution 512², based on the random initial field – figure 2(b).</figcaption></figure>

define the “ADSD” small-scale velocity. For the ensembles with perturbed initial conditions, similar sampling methods are adopted. But the obtained “small-scale” random fields are used only at the initial time. As an example, figure 5 shows the initial perturbation spectrum used in the “Spectral” method for the first initial condition (109). Then, it is added to the initial condition b(x, t = 0). Besides, the dissipation-adapted noise 5.3.2 has been implemented in SQG_MUMATLAB code. For this method – denoted “WavHypervis” – we use a wavelet basis (e_k)_{k∈Z}. We have forecast 5 ensembles with 5 different values of the scaling factor ζ. For instance, figure 6 shows us that under the same low-resolution, when the scaling factor ζ increases (i.e. when the noise
<figure><img src="image_7.png" /><figcaption>Fig. 5 Kinetic energy spectrum ($m^2.s^{-4}/(rad.m^{-1})$), at the initial time, of the mean buoyancy, in blue, spectrum of its random perturbation, in red, and slope $-5/3$ in black. The initial perturbation is restricted to a narrow spectral band. This random initial condition has been used to simulate an ensemble with the deterministic SQG model.</figcaption></figure>

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

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

$$MSB \approx \frac{N_e + 1}{N_e} MEV, \qquad (111)$$

where

$$MSB \triangleq \frac{1}{M} \sum_{j=1}^{M} (b^{o} - \hat{\mathbb{E}}b)^{2} (t, \mathbf{x}_{j}), \qquad (112)$$

$$MEV \triangleq \frac{1}{M} \sum_{j=1}^{M} \widehat{Var}(b)(t, \mathbf{x}_j), \qquad (113)$$

$M$ the number of grid points $\mathbf{x}_j$, $\hat{\mathbb{E}}\{f\}$ and $\widehat{Var}(f)$ the empirical mean and the empirical variance of $f$, computed from the ensemble, respectively. At the same time, a classical error metric – the ensemble mean square error (MSE) is also considered :

$$MSE(t) = \hat{\mathbb{E}} \left( \frac{1}{M} \sum_{j=1}^{M} (b^{o} - b)^{2} (t, \mathbf{x}_{j}) \right). \qquad (114)$$
<figure><img src="image_7.png" /><figcaption>Fig. 5 Kinetic energy spectrum ($m^2.s^{-4}/(rad.m^{-1})$), at the initial time, of the mean buoyancy, in blue, spectrum of its random perturbation, in red, and slope $-5/3$ in black. The initial perturbation is restricted to a narrow spectral band. This random initial condition has been used to simulate an ensemble with the deterministic SQG model.</figcaption></figure>

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

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

$$MSB \approx \frac{N_e + 1}{N_e} MEV, \qquad (111)$$

where

$$MSB \triangleq \frac{1}{M} \sum_{j=1}^{M} (b^{o} - \hat{\mathbb{E}}b)^{2} (t, \mathbf{x}_{j}), \qquad (112)$$

$$MEV \triangleq \frac{1}{M} \sum_{j=1}^{M} \widehat{Var}(b)(t, \mathbf{x}_j), \qquad (113)$$

$M$ the number of grid points $\mathbf{x}_j$, $\hat{\mathbb{E}}\{f\}$ and $\widehat{Var}(f)$ the empirical mean and the empirical variance of $f$, computed from the ensemble, respectively. At the same time, a classical error metric – the ensemble mean square error (MSE) is also considered :

$$MSE(t) = \hat{\mathbb{E}} \left( \frac{1}{M} \sum_{j=1}^{M} (b^{o} - b)^{2} (t, \mathbf{x}_{j}) \right). \qquad (114)$$
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Comparison of one realisation of the buoyancy field after 17 days of advection, according to various energy scaling factor in "WavHypervis".</figcaption></figure>

In figures 7(a) and 7(b), we compare the criterions (111) and (114) pro-
duced by our models. For the 4 vortices field 2(a), it seems that for all models
produce slight errors during the first 10 days. Their forecast skills tend to
a stationary state after about 17 days' advections. In this case, *WavHypervis*
and *SVDpseudo* under location uncertainty have a better ensemble spread, yet
conversely the homogeneous models – *Spectral* and *ADSD* have lower MSE.
As illustrated in figures 7(c) and 7(d), in the *WavHypervis* method, when the
factor $\zeta$ increases, both the spreading and the ensemble errors increase. In fig-
ures 7(e) and 7(f), we compare the results with the spectral initial field 2(b),
in which the ensembles models tend to a stationary state much more faster. In
this case, *WavHypervis* and *SVDpseudo* under location uncertainty have still
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Comparison of one realisation of the buoyancy field after 17 days of advection, according to various energy scaling factor in "WavHypervis".</figcaption></figure>

In figures 7(a) and 7(b), we compare the criterions (111) and (114) pro-
duced by our models. For the 4 vortices field 2(a), it seems that for all models
produce slight errors during the first 10 days. Their forecast skills tend to
a stationary state after about 17 days' advections. In this case, *WavHypervis*
and *SVDpseudo* under location uncertainty have a better ensemble spread, yet
conversely the homogeneous models – *Spectral* and *ADSD* have lower MSE.
As illustrated in figures 7(c) and 7(d), in the *WavHypervis* method, when the
factor $\zeta$ increases, both the spreading and the ensemble errors increase. In fig-
ures 7(e) and 7(f), we compare the results with the spectral initial field 2(b),
in which the ensembles models tend to a stationary state much more faster. In
this case, *WavHypervis* and *SVDpseudo* under location uncertainty have still
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Comparisons of ensemble forecasts for different stochastic models based on the two tested initial conditions : (a) – The mean square error produced by different models under the background of 4 vortices; (b) – The corresponding difference between mean square bias and mean ensemble variance; (c), (d) – Using the empirical noise (5.3.2) with wavelet basis and various energy scaling factors, under the background of 4 vortices; (e), (f) – Results under the spectral background. Notice that all these results are normalized by the amplitude $-B_0$- of the referent initial buoyancy fields.</figcaption></figure>

better ensemble spread. The forecasts by *Spectral* and *ADSD* consist of less errors in ensemble.

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

As explained above, distinct models yield distinct spreading based criterions (111), (114) and Talagrand diagram. Nonetheless, these two criterions
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Comparisons of ensemble forecasts for different stochastic models based on the two tested initial conditions : (a) – The mean square error produced by different models under the background of 4 vortices; (b) – The corresponding difference between mean square bias and mean ensemble variance; (c), (d) – Using the empirical noise (5.3.2) with wavelet basis and various energy scaling factors, under the background of 4 vortices; (e), (f) – Results under the spectral background. Notice that all these results are normalized by the amplitude $-B_0$-of the referent initial buoyancy fields.</figcaption></figure>

better ensemble spread. The forecasts by *Spectral* and *ADSD* consist of less errors in ensemble.

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

As explained above, distinct models yield distinct spreading based criterions (111), (114) and Talagrand diagram. Nonetheless, these two criterions
<figure><img src="image_10.png" /><figcaption>**Table 1** Talagrand diagram on some forecast days under the background of 4 vortices : (a) – MU Spectral; (b) – MU ADSD; (c) – WavHypervis 50%; (d) – MU SVDpseudo; (e) – PIC Spectral; (f) – PIC SVDpseudo.</figcaption></figure>

type provide opposite conclusions. An objective analysis need hence others metrics like proper score (see section 7).

To begin, we will compare the numerical results with the proper score CRPS (see the section 7.3). The CRPS is defined as a pointwise score. Applied on ensemble of spatio-temporal fields, maps of CRPS can be represented at each fixed time step as in the figure 8. As expected on all maps of the figure 8, the normalized CRPS is relatively high on small scale structures and low on the center of each vortices. Indeed, turbulent structures are the hardest things to reproduced. Moreover, the figure 8 suggests that the model *Spectral* in *SQG<sub>PIC</sub>* is the worst model since it has the largest CRPS values. Model *WavHypervis* in *SQG<sub>MU</sub>* shows smaller CRPS. Thus reveals that this model makes fewer local errors at this advection time.
<figure><img src="image_10.png" /><figcaption>**Table 1** Talagrand diagram on some forecast days under the background of 4 vortices : (a) – MU Spectral; (b) – MU ADSD; (c) – WavHypervis 50%; (d) – MU SVDpseudo; (e) – PIC Spectral; (f) – PIC SVDpseudo.</figcaption></figure>

type provide opposite conclusions. An objective analysis need hence others metrics like proper score (see section 7).

To begin, we will compare the numerical results with the proper score CRPS (see the section 7.3). The CRPS is defined as a pointwise score. Applied on ensemble of spatio-temporal fields, maps of CRPS can be represented at each fixed time step as in the figure 8. As expected on all maps of the figure 8, the normalized CRPS is relatively high on small scale structures and low on the center of each vortices. Indeed, turbulent structures are the hardest things to reproduced. Moreover, the figure 8 suggests that the model *Spectral* in *SQG<sub>PIC</sub>* is the worst model since it has the largest CRPS values. Model *WavHypervis* in *SQG<sub>MU</sub>* shows smaller CRPS. Thus reveals that this model makes fewer local errors at this advection time.
<figure><img src="image_11.png" /><figcaption>Fig. 8 CRPS of each model after 20 days of advection with an initial condition based on a smooth field with four given vortices</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 8 CRPS of each model after 20 days of advection with an initial condition based on a smooth field with four given vortices</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Table 2** Talagrand diagram on some forecast days under the spectral background : (a) – MU Spectral; (b) – MU ADSD; (c) – WavHypervis 50%; (d) – MU SVDpseudo; (e) – PIC Spectral; (f) – PIC SVDpseudo.</figcaption></figure>

The CRPS has only been computed at one fixed time. So, again, it is difficult to conclude on which model is better.

To push further the analysis, we maps at different times. The figure 9 reveals a part of the spatial evolution of those CRPS maps for the model *Spectral* in *SQG<sub>MU</sub>*. With the random gaussian field as initial condition, we can see on the figure 9 that the CRPS structures are first small and spread over the spatial domain. Then, these small scale CRPS structures merge and create larger structures of larger intensity. The merging is due to the inverse energy cascade of SQG turbulence (Lapeyre, 2017) which merges the badly resolved small-scale turbulence structures. On the top of that, the new large vortices have chaotic trajectories. The difficulty for the models to track these trajectories yield large CRPS values in the centers of those vortices.
<figure><img src="image_12.png" /><figcaption>**Table 2** Talagrand diagram on some forecast days under the spectral background : (a) – MU Spectral; (b) – MU ADSD; (c) – WavHypervis 50%; (d) – MU SVDpseudo; (e) – PIC Spectral; (f) – PIC SVDpseudo.</figcaption></figure>

The CRPS has only been computed at one fixed time. So, again, it is difficult to conclude on which model is better.

To push further the analysis, we maps at different times. The figure 9 reveals a part of the spatial evolution of those CRPS maps for the model *Spectral* in *SQG*<sub>MU</sub>. With the random gaussian field as initial condition, we can see on the figure 9 that the CRPS structures are first small and spread over the spatial domain. Then, these small scale CRPS structures merge and create larger structures of larger intensity. The merging is due to the inverse energy cascade of SQG turbulence (Lapeyre, 2017) which merges the badly resolved small-scale turbulence structures. On the top of that, the new large vortices have chaotic trajectories. The difficulty for the models to track these trajectories yield large CRPS values in the centers of those vortices.
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Normalized CRPS of the model mu spectrum during 5,10 and 20 days of advection with an initial condition based on an initial random gaussian field.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Normalized CRPS of the model mu spectrum during 5,10 and 20 days of advection with an initial condition based on an initial random gaussian field.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 10** Box plot of the normalized CRPS computed on the ensemble at each spatial points of each model at different time step of advection for two different initial condition (left: Random gaussian field, right: Four given vortices).</figcaption></figure>

Our analysis describes the variability of the CRPS in time and space. But, the difficulty to select a model remains to summarize the CRPS information.

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

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

The model *SVDpseudo* in *SQG<sub>PIC</sub>* has a low energy score after two weeks of advection for the spectral background and after three weeks for the four-vortice initial condition. With the spectral background, after two weeks large-scale structures has formed and the model *SVDpseudo* in *SQG<sub>PIC</sub>* makes less error with the reference. Nevertheless, during the first week of advection,
<figure><img src="image_14.png" /><figcaption>**Fig. 10** Box plot of the normalized CRPS computed on the ensemble at each spatial points of each model at different time step of advection for two different initial condition (left: Random gaussian field, right: Four given vortices).</figcaption></figure>

Our analysis describes the variability of the CRPS in time and space. But, the difficulty to select a model remains to summarize the CRPS information.

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

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

The model *SVDpseudo* in *SQG<sub>PIC</sub>* has a low energy score after two weeks of advection for the spectral background and after three weeks for the four-vortice initial condition. With the spectral background, after two weeks large-scale structures has formed and the model *SVDpseudo* in *SQG<sub>PIC</sub>* makes less error with the reference. Nevertheless, during the first week of advection,
<figure><img src="image_15.png" /><figcaption>**Fig. 11** Normalized energy score for each model and spatial multivariate ensemble on all time steps and both initial condition (left: Random gaussian field, right: Four given vortices).</figcaption></figure>

the model *SVDpseudo* in *SQG<sub>PIC</sub>* has a really large energy score because, small-scale structures are present and are not well reproduced by *SVDpseudo* in *SQG<sub>PIC</sub>*. So, in this case, the models *SQG<sub>PIC</sub>* are less efficient than the models *SQG<sub>MU</sub>* and *WavHypervis*.

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

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

We will now consider proper scores and variograms where the temporal structures are considered instead of the spatial ones. The figure 12 reflects that the model *Spectral* in *SQG<sub>PIC</sub>* has the largest normalized energy score and variogram-0.5 score (see section 7.5). On the edges of the centers of the vortices, this model encompasses the worst temporal structures. The models *SVDpseudo* in *SQG<sub>MU</sub>* and *WavHypervis* seem to perform better. Nonetheless, in the model *Spectral* in *SQG<sub>MU</sub>* minimum values are located in the center of the vortices which indicates that the temporal evolution of the vortices is well reproduced by this model.
<figure><img src="image_15.png" /><figcaption>**Fig. 11** Normalized energy score for each model and spatial multivariate ensemble on all time steps and both initial condition (left: Random gaussian field, right: Four given vortices).</figcaption></figure>

the model *SVDpseudo* in *SQG<sub>PIC</sub>* has a really large energy score because, small-scale structures are present and are not well reproduced by *SVDpseudo* in *SQG<sub>PIC</sub>*. So, in this case, the models *SQG<sub>PIC</sub>* are less efficient than the models *SQG<sub>MU</sub>* and *WavHypervis*.

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

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

We will now consider proper scores and variograms where the temporal structures are considered instead of the spatial ones. The figure 12 reflects that the model *Spectral* in *SQG<sub>PIC</sub>* has the largest normalized energy score and variogram-0.5 score (see section 7.5). On the edges of the centers of the vortices, this model encompasses the worst temporal structures. The models *SVDpseudo* in *SQG<sub>MU</sub>* and *WavHypervis* seem to perform better. Nonetheless, in the model *Spectral* in *SQG<sub>MU</sub>* minimum values are located in the center of the vortices which indicates that the temporal evolution of the vortices is well reproduced by this model.
<figure><img src="image_16.png" /><figcaption>**Fig. 12** Normalized energy score and variogram-0.5 score of the temporal multivariate ensemble of each models at each spatial points with an initial condition based on the four vortices (left column: normalized energy score of each model, right: normalized variogram-0.5 score of each model).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 12** Normalized energy score and variogram-0.5 score of the temporal multivariate ensemble of each models at each spatial points with an initial condition based on the four vortices (left column: normalized energy score of each model, right: normalized variogram-0.5 score of each model).</figcaption></figure>
<table><caption>**Table 3** Model performance by score, MSE: Mean squared error (pointwise), MSB: Mean squared bias (pointwise), CRPS: Continuous ranked proper score (pointwise), ES: Energy score (spatial and time).</caption><thead><tr><th>Model</th><th>MSE</th><th>MSB</th><th>CRPS</th><th>ES</th><th>VS</th></tr></thead><tbody><tr><td>SQG<sub>MUSPEC</sub></td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td></tr><tr><td>SQG<sub>MUADSD</sub></td><td>++</td><td>+</td><td>+</td><td>+</td><td>+</td></tr><tr><td>SQG<sub>MUSVD</sub></td><td>-</td><td>++</td><td>++</td><td>++</td><td>++</td></tr><tr><td>WavHypervis<sub>50</sub></td><td>-</td><td>+</td><td>++</td><td>++</td><td>++</td></tr><tr><td>SQG<sub>PICSPEC</sub></td><td>-</td><td>--</td><td>--</td><td>--</td><td>-</td></tr><tr><td>SQG<sub>PICSVD</sub></td><td>--</td><td>-</td><td>+</td><td>+</td><td>+</td></tr></tbody></table>

## 8.2.1 Conclusion of the short-term forecasts analysis

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

In the table 3, $SQG_{PIC}$ methods obtained the lowest overall verification scores. In contrast, ADSD in $SQG_{MU}$, SVDpseudo in $SQG_{MU}$ and WavHypervis models present the best performances. Therefore, we recommend one of this method for UQ tasks, and we strongly advice to avoid relying only on initial conditions randomization.

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

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

## 8.3 Long-term ensemble forecast

After the first month of advection, the flow initiated with (109) (figure 3) breaks its symmetries. This leads to a chaotic behavior and a complete loss
<table><caption>**Table 3** Model performance by score, MSE: Mean squared error (pointwise), MSB: Mean squared bias (pointwise), CRPS: Continuous ranked proper score (pointwise), ES: Energy score (spatial and time).</caption><thead><tr><th>Model</th><th>MSE</th><th>MSB</th><th>CRPS</th><th>ES</th><th>VS</th></tr></thead><tbody><tr><td>SQG<sub>MUSPEC</sub></td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td></tr><tr><td>SQG<sub>MUADSD</sub></td><td>++</td><td>+</td><td>+</td><td>+</td><td>+</td></tr><tr><td>SQG<sub>MUSVD</sub></td><td>-</td><td>++</td><td>++</td><td>++</td><td>++</td></tr><tr><td>WavHypervis<sub>50</sub></td><td>-</td><td>+</td><td>++</td><td>++</td><td>++</td></tr><tr><td>SQG<sub>PICSPEC</sub></td><td>-</td><td>--</td><td>--</td><td>--</td><td>-</td></tr><tr><td>SQG<sub>PICSVD</sub></td><td>--</td><td>-</td><td>+</td><td>+</td><td>+</td></tr></tbody></table>

### 8.2.1 Conclusion of the short-term forecasts analysis

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

In the table 3, $SQG_{PIC}$ methods obtained the lowest overall verification scores. In contrast, ADSD in $SQG_{MU}$, SVDpseudo in $SQG_{MU}$ and WavHypervis models present the best performances. Therefore, we recommend one of this method for UQ tasks, and we strongly advice to avoid relying only on initial conditions randomization.

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

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

## 8.3 Long-term ensemble forecast

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

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

### 8.3.1 Chaotic test flow and resolution issues

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

To trigger a particular transition, we introduce two infinitesimal modifications in the initial condition. For the sake of simplicity, we focus on the bifurcation associated with the warm vortices. To do so, the merging of cold vortices in $x = 500$ km will be forced, by adding an infinitesimal small-scale cold eddy in $(x, y) = (480, 750)$ (in km). The barycentre of northern structures becomes slightly closer to $x = 500$ km. This gives rise to the desired transition, as shown in Figures 13 and 14. To trigger the bifurcation associated with the southern warm vortices, an infinitesimal small-scale eddy has been added in $(x, y) = (20, 250)$ (in km). If the eddy is warm, the southern barycentre is moved closer to $x = 0$, and the two warm vortices merge near $x = 0$ (scenario 1), as shown in Figure 13. If this eddy is cold, the southern barycentre is moved closer to $x = 500$ km, and the two warm vortices merge near $x = 500$ km (scenario 2), as shown in Figure 14. The exact expression of the initial
of (deterministic) predictability. In this section, we will focus on this regime. As demonstrated here, the bifurcation associated with the symmetry breaking is efficiently tracked using the model under location uncertainty. Using few realizations, the probability density functions of each subsequent scenarios is well characterized, whereas methods based on random initial conditions do not converge. The identification of several scenarios, done at each time step in a reduced subspace, is obtained by Principal Component Analysis (PCA), also termed Empirical Orthogonal Functions method (EOF).

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

### 8.3.1 Chaotic test flow and resolution issues

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

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

$$b(x, t = 0) = b_0(x) + 0.3s_w F \left( \frac{1}{40} \left( x - \begin{pmatrix} 20 \text{ km} \\ 250 \text{ km} \end{pmatrix} \right) \right) - 0.3F \left( \frac{1}{40} \left( x - \begin{pmatrix} 480 \text{ km} \\ 750 \text{ km} \end{pmatrix} \right) \right), \quad (115)$$

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

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

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

### 8.3.2 Stochastic analysis

Unlike the deterministic SQG model, the SQG dynamics under location uncertainty, with the exact same initialization $s_w = 1$, yields several likely transitions. Here, for simplicity, only the variant Spectral of SQG$_{MU}$ is considered (neither ADSD nor SVDpseudo). The calibration Spectral of SQG$_{MU}$ is detailed in section 6.3.2. In Figure 16, we show two realizations of the SQG$_{MU}$ dynamics. One of those realizations corresponds to the reference scenario (scenario 1), the other does not. The model encodes several likely transitions, and
condition is the following:

$$b(\boldsymbol{x}, t = 0) = b_0(\boldsymbol{x}) + 0.3s_w F \left( \frac{1}{40} \left( \boldsymbol{x} - \begin{pmatrix} 20 \text{ km} \\ 250 \text{ km} \end{pmatrix} \right) \right) - 0.3F \left( \frac{1}{40} \left( \boldsymbol{x} - \begin{pmatrix} 480 \text{ km} \\ 750 \text{ km} \end{pmatrix} \right) \right), \quad (115)$$

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

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

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

### 8.3.2 Stochastic analysis

Unlike the deterministic SQG model, the SQG dynamics under location uncertainty, with the exact same initialization $s_w = 1$, yields several likely transitions. Here, for simplicity, only the variant Spectral of SQG$_{MU}$ is considered (neither ADSD nor SVDpseudo). The calibration Spectral of SQG$_{MU}$ is detailed in section 6.3.2. In Figure 16, we show two realizations of the SQG$_{MU}$ dynamics. One of those realizations corresponds to the reference scenario (scenario 1), the other does not. The model encodes several likely transitions, and
Scenario 1

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

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

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

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

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

line) and the mean spectrum of the tracer random component, $\hat{\mathbb{E}}\{\Gamma_{b-\hat{\mathbb{E}}\{b\}}\}$,
Scenario 2

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

line) and the mean spectrum of the tracer random component, $\hat{\mathbb{E}}\{\Gamma_{b-\hat{\mathbb{E}}\{b\}}\}$,
<figure><img src="image_19.png" /><figcaption>Fig. 15 Buoyancy (m.s⁻²) after 70 days of advection for the SQG model at resolution 512² (left) and at resolution 128² (right). The small-scale perturbation in the initial condition activating the symmetry breaking in the reference simulation (512²), is dissipated in few days in the low-resolution simulation. This makes this single low-resolution simulation erroneous.</figcaption></figure>

(shaded grey):

$$ \hat{\mathbb{E}}\{\Gamma_b\} = \hat{\mathbb{E}}\left\{\frac{1}{\mu(\Omega)} \oint \|\mathbf{k}\|_2 |\hat{b}|^2 d\theta_k\right\}, \quad (116) $$

$$ = \frac{1}{\mu(\Omega)} \oint \|\mathbf{k}\|_2 \underbrace{\left( |\hat{\mathbb{E}}\{\hat{b}\}|^2 + \hat{\mathbb{E}}\{ |\hat{b} - \hat{\mathbb{E}}\{\hat{b}\}|^2 \} \right)}_{\text{Energy of the mean}} d\theta_k, \quad (117) $$

$$ = \Gamma_{\hat{\mathbb{E}}\{b\}} + \hat{\mathbb{E}}\{\Gamma_{b-\hat{\mathbb{E}}\{b\}}\}, \quad (118) $$

where $\mu(\Omega)$ denotes the area of the domain $\Omega$, $\hat{\mathbb{E}}\{f\}$ the empirical mean of $f$, computed from the ensemble, and $\hat{f}$ the Fourier transform of $f$.

<figure><img src="image_20.png" /><figcaption>Fig. 16 Two realizations of buoyancy (m.s⁻²) after 70 days of advection for the SQG under location uncertainty at resolution 128². Event though the small-scale perturbation in the initial condition is dissipated in few days, the small-scale component of the dynamics under location uncertainty randomly triggers the symmetry breaking. Therefore, some realizations follow the right transitions (left) and some do not (right).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 15 Buoyancy (m.s⁻²) after 70 days of advection for the SQG model at resolution 512² (left) and at resolution 128² (right). The small-scale perturbation in the initial condition activating the symmetry breaking in the reference simulation (512²), is dissipated in few days in the low-resolution simulation. This makes this single low-resolution simulation erroneous.</figcaption></figure>

(shaded grey):

$$ \hat{\mathbb{E}}\{\Gamma_b\} = \hat{\mathbb{E}}\left\{\frac{1}{\mu(\Omega)} \oint \|\mathbf{k}\|_2 |\hat{b}|^2 d\theta_k\right\}, \quad (116) $$

$$ = \frac{1}{\mu(\Omega)} \oint \|\mathbf{k}\|_2 \left( \underbrace{|\hat{\mathbb{E}}\{\hat{b}\}|^2}_{\text{Energy of the mean}} + \underbrace{\hat{\mathbb{E}}\{|\hat{b} - \hat{\mathbb{E}}\{\hat{b}\}|^2\}}_{\text{Variance}} \right) d\theta_k, \quad (117) $$

$$ = \Gamma_{\hat{\mathbb{E}}\{b\}} + \hat{\mathbb{E}}\{\Gamma_{b-\hat{\mathbb{E}}\{b\}}\}, \quad (118) $$

where $\mu(\Omega)$ denotes the area of the domain $\Omega$, $\hat{\mathbb{E}}\{f\}$ the empirical mean of $f$, computed from the ensemble, and $\hat{f}$ the Fourier transform of $f$.

<figure><img src="image_20.png" /><figcaption>Fig. 16 Two realizations of buoyancy (m.s⁻²) after 70 days of advection for the SQG under location uncertainty at resolution 128². Event though the small-scale perturbation in the initial condition is dissipated in few days, the small-scale component of the dynamics under location uncertainty randomly triggers the symmetry breaking. Therefore, some realizations follow the right transitions (left) and some do not (right).</figcaption></figure>
Uncertainty backscatter

<figure><img src="image_21.png" /><figcaption>**Fig. 17** Tracer spectrum ($m^2 \cdot s^{-4} / (\text{rad.m}^{-1})$) after 30 days (top) and 70 days (bottom) of advection for SQG model at resolution 512² (green), one realization of SQG_MU model, $\Gamma_{b(1)}$, at resolution 128² (red dashed line), the spectrum of the mean, $\Gamma_{\hat{E}\{b\}}$, (blue line) and the mean spectrum of the tracer random component, $\hat{E}\{\Gamma_{b-\hat{E}\{b\}}\}$, (shaded grey). After being stacked, the two last plots represent the mean spectrum: $\hat{E}\{\Gamma_b\} = \Gamma_{\hat{E}\{b\}} + \hat{E}\{\Gamma_{b-\hat{E}\{b\}}\}$. The more thick the shaded grey area is, the more variance is contained at this scale. At $t=30$ days, the variance remains at small scales but this small-scale uncertainty activates the symmetry breaking. This results in a variance backscattering with a thickening of the spectrum of the random component at large scales (visualized at $t=70$). Since large scales influence strongly the small scales the small-scale variance is also enforced.</figcaption></figure>

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

<figure><img src="image_21.png" /><figcaption>**Fig. 17** Tracer spectrum ($m^2 \cdot s^{-4} / (\text{rad.m}^{-1})$) after 30 days (top) and 70 days (bottom) of advection for SQG model at resolution 512² (green), one realization of SQG_MU model, $\Gamma_b^{(1)}$, at resolution 128² (red dashed line), the spectrum of the mean, $\Gamma_{\hat{E}\{b\}}$, (blue line) and the mean spectrum of the tracer random component, $\hat{E}\{\Gamma_{b-\hat{E}\{b\}}\}$, (shaded grey). After being stacked, the two last plots represent the mean spectrum: $\hat{E}\{\Gamma_b\} = \Gamma_{\hat{E}\{b\}} + \hat{E}\{\Gamma_{b-\hat{E}\{b\}}\}$. The more thick the shaded grey area is, the more variance is contained at this scale. At $t=30$ days, the variance remains at small scales but this small-scale uncertainty activates the symmetry breaking. This results in a variance backscattering with a thickening of the spectrum of the random component at large scales (visualized at $t=70$). Since large scales influence strongly the small scales the small-scale variance is also enforced.</figcaption></figure>

One can then wonder about the number of realizations following the right scenario. This necessitates the analysis of the 4-dimensional spatio-temporal random field $(b^{(i)}(x_j, y_k, t_l))_{ijkl}$. The superscript (i) designates the i-th realization of the ensemble. To reduce the associated dimension, a Principal Component Analysis (PCA) – also termed Empirical Orthogonal Function (EOF) representation – is performed over the realizations, at a fixed time $t$. Within this analysis, unlike usual EOF representation, the time axis is replaced by the realization index. At a given time, it helps to represent the whole ensemble (200 realizations of the random field) by the ensemble mean field and few other
EOF spatial modes. In the present case, the buoyancy is approximated as:

$$b^{(i)}(\mathbf{x}, t) = \hat{\mathbb{E}}(b)(\mathbf{x}, t) + \sum_{n=1}^{N_e} c_n^{(i)}(t) \Psi(\mathbf{x}, t), \quad (119)$$

$$\approx \hat{\mathbb{E}}(b)(\mathbf{x}, t) + \sum_{n=1}^{N_{EOF}} c_n^{(i)}(t) \Psi_n(\mathbf{x}, t), \quad (120)$$

where $N_e$ is the size of the ensemble, $N_{EOF} \ll N_e$ is the number of EOF modes chosen to described the whole ensemble and the $(\Psi_n)_{1\le n\le N_{EOF}}$ denote the EOF spatial modes. Those spatial fields are orthogonal:

$$\frac{1}{M} \sum_{j=1}^{M} \Psi_n(\mathbf{x}_j, t) \Psi_m(\mathbf{x}_j, t) = \delta_{nm}, \quad (121)$$

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

$$\underbrace{\hat{\mathbb{E}}\left\{\frac{1}{M} \sum_{j=1}^{M} b^2(\mathbf{x}_j, t)\right\}}_{\text{Mean energy}} = \frac{1}{M} \sum_{j=1}^{M} \left(\hat{\mathbb{E}}(b)\right)^2(\mathbf{x}_j, t) + \frac{1}{M} \sum_{j=1}^{M} \widehat{Var}(b)(\mathbf{x}_j, t), \quad (122)$$

$$= \frac{1}{M} \sum_{j=1}^{M} \left( \hat{\mathbb{E}}(b) \right)^2 (\mathbf{x}_j, t) + \underbrace{\sum_{n=1}^{N_e} \hat{\mathbb{E}}\{c_n^2(t)\}}_{\text{Mean energy of the } n\text{-th EOF}}. \quad (123)$$

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

$$b^{(i)}(\mathbf{x}, t) = \hat{\mathbb{E}}(b)(\mathbf{x}, t) + \sum_{n=1}^{N_e} c_n^{(i)}(t) \Psi(\mathbf{x}, t), \quad (119)$$

$$\approx \hat{\mathbb{E}}(b)(\mathbf{x}, t) + \sum_{n=1}^{N_{EOF}} c_n^{(i)}(t) \Psi_n(\mathbf{x}, t), \quad (120)$$

where $N_e$ is the size of the ensemble, $N_{EOF} \ll N_e$ is the number of EOF modes chosen to described the whole ensemble and the $(\Psi_n)_{1\le n\le N_{EOF}}$ denote the EOF spatial modes. Those spatial fields are orthogonal:

$$\frac{1}{M} \sum_{j=1}^{M} \Psi_n(\mathbf{x}_j, t) \Psi_m(\mathbf{x}_j, t) = \delta_{nm}, \quad (121)$$

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

$$\underbrace{\hat{\mathbb{E}}\left\{\frac{1}{M} \sum_{j=1}^{M} b^2(\mathbf{x}_j, t)\right\}}_{\text{Mean energy}} = \frac{1}{M} \sum_{j=1}^{M} \left(\hat{\mathbb{E}}(b)\right)^2 (\mathbf{x}_j, t) + \frac{1}{M} \sum_{j=1}^{M} \widehat{\mathrm{Var}}(b)(\mathbf{x}_j, t), \quad (122)$$

$$= \frac{1}{M} \sum_{j=1}^{M} \left(\hat{\mathbb{E}}(b)\right)^2 (\mathbf{x}_j, t) + \underbrace{\sum_{n=1}^{N_e} \hat{\mathbb{E}}\{c_n^2(t)\}}_{\text{Mean energy of the } n\text{-th EOF}}. \quad (123)$$

In Figure 18, left part, the energy associated with the different spatial modes is displayed. At $t = 30$ days, the energy of the mean field, denoted as a 0-th order EOF (index+1 = 1), is much larger than the variance field, described by the other EOFs. On the contrary, at $t = 70$ days, after the breaking symmetry, the energy of the mean and the variance have the same order of magnitude. The variance is mainly explained by the first EOF. Thus, at the first order, the randomness of the tracer is approximately encoded by this first EOF coefficient, $c_1$. Its probability density function and the joint probability density function for the two first EOF coefficients presented in Figure 18 are unimodal at $t = 30$ days and bimodal at $t = 70$ days. Note that the tracer is clearly non-Gaussian. The symmetry breaking has created two likely scenarios in the ensemble. The scenario corresponding to negative values of the first EOF coefficient (probability of 47%) is called scenario A and the scenario corresponding to positive values (probability of 53%) is called scenario
<figure><img src="image_22.png" /><figcaption>**Fig. 18** Energy of the EOFs ($m^2.s^{-4}$) (left), probability density function of the first EOF coefficient (middle) and joint probability density function for the two first EOF coefficients (right) after 30 days (top) and 70 days (bottom). The ensemble is simulated according to the SQG model under moderate uncertainty at resolution $128^2$. At $t = 30$ days, the energy of the mean field, denoted as a 0-th order EOF (index+1 = 1), is much higher than the variance field, described by the other EOFs. On the contrary, at $t = 70$ days (after the symmetry breaking), the energy of the mean and variance have the same order of magnitude. The variance is mainly explained by the first EOF. Therefore, at the first order, the random component of the tracer can be approximated by this first EOF. The probability density function of the first EOF coefficient and the joint probability density function for the two first EOF coefficients are unimodal at $t = 30$ days and bimodal at $t = 70$ days. The breaking symmetry has created two likely scenarios, which are very different from one another. The scenario A corresponds to negative values of the first EOF coefficient (probability of 47%) whereas the scenario B corresponds to positive value (probability of 53%). The red line separates the probability density function between the two scenarios.</figcaption></figure>

B. The red line at zero separates the probability density function between the two scenarios. In Figure 19, the same probability density function along time is plotted. The bifurcation is clearly visible. Also shown, the mean buoyancy ($m.s^{-2}$) of the two likely scenarios are represented after 70 days of advection.
<figure><img src="image_22.png" /><figcaption>**Fig. 18** Energy of the EOFs ($m^2.s^{-4}$) (left), probability density function of the first EOF coefficient (middle) and joint probability density function for the two first EOF coefficients (right) after 30 days (top) and 70 days (bottom). The ensemble is simulated according to the SQG model under moderate uncertainty at resolution $128^2$. At $t = 30$ days, the energy of the mean field, denoted as a 0-th order EOF (index+1 = 1), is much higher than the variance field, described by the other EOFs. On the contrary, at $t = 70$ days (after the symmetry breaking), the energy of the mean and variance have the same order of magnitude. The variance is mainly explained by the first EOF. Therefore, at the first order, the random component of the tracer can be approximated by this first EOF. The probability density function of the first EOF coefficient and the joint probability density function for the two first EOF coefficients are unimodal at $t = 30$ days and bimodal at $t = 70$ days. The breaking symmetry has created two likely scenarios, which are very different from one another. The scenario A corresponds to negative values of the first EOF coefficient (probability of 47%) whereas the scenario B corresponds to positive value (probability of 53%). The red line separates the probability density function between the two scenarios.</figcaption></figure>

B. The red line at zero separates the probability density function between the two scenarios. In Figure 19, the same probability density function along time is plotted. The bifurcation is clearly visible. Also shown, the mean buoyancy ($m.s^{-2}$) of the two likely scenarios are represented after 70 days of advection.
The two mean fields, $\hat{\mathbb{E}}(b|1)$ and $\hat{\mathbb{E}}(b|2)$, are calculated as:

$$ \hat{\mathbb{E}}(b|1) = \hat{\mathbb{E}}(b) + \hat{\mathbb{E}}(c_1|1)\Psi_1 = \hat{\mathbb{E}}(b) + \left( \frac{1}{\#\{i|c_1^{(i)} < 0\}} \sum_{c_1^{(i)} < 0} c_1^{(i)} \right) \Psi_1, \quad (124) $$

$$ \hat{\mathbb{E}}(b|2) = \hat{\mathbb{E}}(b) + \hat{\mathbb{E}}(c_1|2)\Psi_1 = \hat{\mathbb{E}}(b) + \left( \frac{1}{\#\{i|c_1^{(i)} \ge 0\}} \sum_{c_1^{(i)} \ge 0} c_1^{(i)} \right) \Psi_1, \quad (125) $$

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

For sake of comparison, we also show results obtained using the deterministic SQG model, with the “Spectral” perturbation of the initial conditions. This method was already used in the previous section for short-term forecast. We first briefly recall how these perturbations are generated. Initial small-scale buoyancy perturbations are assumed Gaussian and sampled from a $(-\frac{5}{3})$ spectrum, as shown in Figure 5. These perturbations should not change the large-scale flow before the predictability time (about one month). Accordingly, we require those perturbations to be of small amplitudes and restrict them to small scales. The same analysis is performed, including the EOF decomposition and the distinction between two likely scenarios. Figure 20 gathers the results. After 70 days of advection, the scenario A is hardly visible in the probability density function of the first EOF. The ensemble estimates a probability of only 39% for the reference scenario. The $SQG_{MU}$ ensemble estimated a probability of 47%. Moreover, the probability density function is very noisy. This suggests that the ensemble may not be converged, i.e. the empirical statistics of the ensemble will change if the ensemble size grows. The scenario A completely disappears in the joint probability density function for the two first EOF coefficients. The probability density function of the first EOF along time, before and after the symmetry breaking, exhibits very narrow branches associated with high probabilities, compared to Figure 19. It indicates that randomized initial conditions may lead to underdispersive ensemble. Furthermore, in Figure 20, trajectories of some realizations are still visible after the bifurcation. This confirms that the probability density function did not converge. Indeed, to estimate this density, we use the well-known Parzen-Rosenblatt estimator (Rosenblatt, 1956; Parzen, 1962): each realization is associated with a kernel
The two mean fields, $\hat{\mathbb{E}}(b|1)$ and $\hat{\mathbb{E}}(b|2)$, are calculated as:

$$ \hat{\mathbb{E}}(b|1) = \hat{\mathbb{E}}(b) + \hat{\mathbb{E}}(c_1|1)\Psi_1 = \hat{\mathbb{E}}(b) + \left( \frac{1}{\#\{i|c_1^{(i)} < 0\}} \sum_{c_1^{(i)} < 0} c_1^{(i)} \right) \Psi_1, \quad (124) $$

$$ \hat{\mathbb{E}}(b|2) = \hat{\mathbb{E}}(b) + \hat{\mathbb{E}}(c_1|2)\Psi_1 = \hat{\mathbb{E}}(b) + \left( \frac{1}{\#\{i|c_1^{(i)} \ge 0\}} \sum_{c_1^{(i)} \ge 0} c_1^{(i)} \right) \Psi_1, \quad (125) $$

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

For sake of comparison, we also show results obtained using the deterministic SQG model, with the “Spectral” perturbation of the initial conditions. This method was already used in the previous section for short-term forecast. We first briefly recall how these perturbations are generated. Initial small-scale buoyancy perturbations are assumed Gaussian and sampled from a $(-\frac{5}{3})$ spectrum, as shown in Figure 5. These perturbations should not change the large-scale flow before the predictability time (about one month). Accordingly, we require those perturbations to be of small amplitudes and restrict them to small scales. The same analysis is performed, including the EOF decomposition and the distinction between two likely scenarios. Figure 20 gathers the results. After 70 days of advection, the scenario A is hardly visible in the probability density function of the first EOF. The ensemble estimates a probability of only 39% for the reference scenario. The $SQG_{MU}$ ensemble estimated a probability of 47%. Moreover, the probability density function is very noisy. This suggests that the ensemble may not be converged, i.e. the empirical statistics of the ensemble will change if the ensemble size grows. The scenario A completely disappears in the joint probability density function for the two first EOF coefficients. The probability density function of the first EOF along time, before and after the symmetry breaking, exhibits very narrow branches associated with high probabilities, compared to Figure 19. It indicates that randomized initial conditions may lead to underdispersive ensemble. Furthermore, in Figure 20, trajectories of some realizations are still visible after the bifurcation. This confirms that the probability density function did not converge. Indeed, to estimate this density, we use the well-known Parzen-Rosenblatt estimator (Rosenblatt, 1956; Parzen, 1962): each realization is associated with a kernel
Stochastic tracking of the bifurcation

<figure><img src="image_23.png" /><figcaption>**Fig. 19** Probability density function of the first EOF coefficient along time (top), buoyancy ($m.s^{-2}$) after 70 days of advection for the mean of the two likely scenarios A and B of the SQG<sub>MU</sub> model at resolution 128<sup>2</sup> (respectively middle left and middle right) and the reference scenarios 1 and 2 at high resolution 512<sup>2</sup> (respectively bottom left and bottom right). The bifurcation is clearly visible on the top plot. The two likely scenarios differ from the sign of the first EOF coefficient. They are almost associated with an equal probability of occurrence: a probability of 47% for scenario A and a probability of 53% for scenario B. The scenario A (respectively B) is similar to the scenario 1 (respectively 2).</figcaption></figure>
Stochastic tracking of the bifurcation

<figure><img src="image_23.png" /><figcaption>**Fig. 19** Probability density function of the first EOF coefficient along time (top), buoyancy ($m.s^{-2}$) after 70 days of advection for the mean of the two likely scenarios A and B of the SQG<sub>MU</sub> model at resolution 128<sup>2</sup> (respectively middle left and middle right) and the reference scenarios 1 and 2 at high resolution 512<sup>2</sup> (respectively bottom left and bottom right). The bifurcation is clearly visible on the top plot. The two likely scenarios differ from the sign of the first EOF coefficient. They are almost associated with an equal probability of occurrence: a probability of 47% for scenario A and a probability of 53% for scenario B. The scenario A (respectively B) is similar to the scenario 1 (respectively 2).</figcaption></figure>
and the estimator is the sum of those kernels. Here, some realizations or set of few realizations are isolated and create spikes in the estimator. More realizations would be needed to have almost continuously distributed realizations. In other words, the ensemble is not converged. This drawback could be expected for at least two reasons. First, the structure of the initial noise contains little physical information, while the dimension of the state space is huge. Without phase information, covering all the possibilities requires a very large number of realizations. Furthermore, the subgrid tensor diffuses the small-scales components of the tracer where the ensemble variability is encoded. This is a known feature of ensemble forecasts: ensemble members tend to align with most unstable directions of the dynamics (Le Dimet and Talagrand, 1986; Ng et al., 2011; Gottwald and Harlim, 2013). Since small scales are stabilized by the subgrid tensor, the ensemble shrinks to span a smaller large-scale unstable subspace (Sapsis, 2013). On the contrary, the stochastic model associates phase and intermittency with the noise and continuously injects it into the dynamics. The phase information or inhomogeneity as well as the non-Gaussianity come from the multiplicative structure. Even though the uncorrelated velocity is only prescribed by a spectrum, the tracer gradients have phase and dynamically constraint the regions of application of the noise. This process makes the stochastic forcing much more efficient. Hence, a smaller number of realizations are needed. In Figure 20, the convergence of the probability density function of scenario B (positive values of the first EOF coefficient) seems slightly better than the density of scenario A. Unfortunately, the bottom Figures shows that the scenario B is not the one followed by the high-resolution simulation. Let us note that the reference is deterministic. Accordingly, the reference probability density is a dirac measure. Indeed, the deterministic reference initial condition is assumed to be known and is used in all large-scale simulations. Moreover, the reference dynamic defined by the high-resolution SQG model is deterministic as the real ocean dynamics is. The bad description of the scenario A tends to confirm that the SQG model with randomized initial conditions fails to describe the bifurcation.

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

### 8.3.3 Conclusion of the long-term forecasts analysis

In this study, long-time forecasts of two different SQG models have been compared. The first one corresponds to the classical SQG equations with a random initial condition. The second one is the SALT-LU SQG, derived from a
and the estimator is the sum of those kernels. Here, some realizations or set of few realizations are isolated and create spikes in the estimator. More realizations would be needed to have almost continuously distributed realizations. In other words, the ensemble is not converged. This drawback could be expected for at least two reasons. First, the structure of the initial noise contains little physical information, while the dimension of the state space is huge. Without phase information, covering all the possibilities requires a very large number of realizations. Furthermore, the subgrid tensor diffuses the small-scales components of the tracer where the ensemble variability is encoded. This is a known feature of ensemble forecasts: ensemble members tend to align with most unstable directions of the dynamics (Le Dimet and Talagrand, 1986; Ng et al., 2011; Gottwald and Harlim, 2013). Since small scales are stabilized by the subgrid tensor, the ensemble shrinks to span a smaller large-scale unstable subspace (Sapsis, 2013). On the contrary, the stochastic model associates phase and intermittency with the noise and continuously injects it into the dynamics. The phase information or inhomogeneity as well as the non-Gaussianity come from the multiplicative structure. Even though the uncorrelated velocity is only prescribed by a spectrum, the tracer gradients have phase and dynamically constraint the regions of application of the noise. This process makes the stochastic forcing much more efficient. Hence, a smaller number of realizations are needed. In Figure 20, the convergence of the probability density function of scenario B (positive values of the first EOF coefficient) seems slightly better than the density of scenario A. Unfortunately, the bottom Figures shows that the scenario B is not the one followed by the high-resolution simulation. Let us note that the reference is deterministic. Accordingly, the reference probability density is a dirac measure. Indeed, the deterministic reference initial condition is assumed to be known and is used in all large-scale simulations. Moreover, the reference dynamic defined by the high-resolution SQG model is deterministic as the real ocean dynamics is. The bad description of the scenario A tends to confirm that the SQG model with randomized initial conditions fails to describe the bifurcation.

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

### 8.3.3 Conclusion of the long-term forecasts analysis

In this study, long-time forecasts of two different SQG models have been compared. The first one corresponds to the classical SQG equations with a random initial condition. The second one is the SALT-LU SQG, derived from a
Same diagnostic based on a deterministic model with random initial conditions

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

<figure><img src="image_24.png" /><figcaption>**Fig. 20** Energy of the EOFs ($m^2.s^{-4}$) (top left), probability density function of the first EOF coefficient (top middle), joint probability density function for the two first EOF coefficients (top right) after 70 days, probability density function of the first EOF coefficient along time (middle), buoyancy ($m.s^{-2}$) after 70 days of advection for the means of the two likely scenarios A and B (bottom left and bottom right respectively). The ensemble is simulated according to the usual deterministic SQG model with random initial conditions at resolution 128². The joint probability density function for the two first EOF coefficients suggests only one likely scenario. The central Figure confirms that this model fails to correctly describe the bifurcation. Indeed, the probability density function appears to be not converged in this case. The bottom Figures show that the worst resolved scenario (the scenario A) is the one similar to the true reference scenario (the scenario 1).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. 21** Minimum normalized error along time in the ensemble with random initial conditions (red) and in the ensemble with random dynamics (blue). The square error was integrated over the space and divided by the energy of the reference.</figcaption></figure>

stochastic expression of the transport equations (see section 6). Both models are compared to a high-resolution simulation reference.

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

To identify and separate the scenarios, a simple threshold is sufficient regarding the high energy distributed along the first PCA axis. For more realistic flow, the number of scenarios is likely to be larger. Hence, this classification may be inadequate and more advanced clustering methods (e.g. k-means and expectation-maximization algorithms) are probably required. In this paper, this ensemble aggregation was used as a diagnosis. But, it could also improve filtering methods. Indeed, for geophysical flows, the state-space dimension being large, the ensemble size is generally too small to encode all possibilities. So, when an observation is assimilated, only few realizations of the ensemble are close enough, and all the others are considered useless. This often leads to
<figure><img src="image_25.png" /><figcaption>**Fig. 21** Minimum normalized error along time in the ensemble with random initial conditions (red) and in the ensemble with random dynamics (blue). The square error was integrated over the space and divided by the energy of the reference.</figcaption></figure>

stochastic expression of the transport equations (see section 6). Both models are compared to a high-resolution simulation reference.

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

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

While the numerical analysis of this section 8 is applied to a random-
ized version of the SQG model with toy initial conditions, the multiplicative
structure of the noise and the balance with diffusion hold for any fluid dynam-
ics models under location uncertainty. This suggests that similar conclusions
could be expected in more complex problems. Hence, this SALT-LU method
opens for new ensemble forecasts methods, for both short-term forecasts and
climate projections where uncertainty quantification is a main issue (Allen and
Stainforth, 2002).

9 Conclusion

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

After presenting a new UQ method (WavHypervis) adaptable to most de-
terministic subgrid dissipation, we have focused on a recent family of stochas-
tic subgrid models: the dynamics under location uncertainty (LU) and the
stochastic advection by Lie transport (SALT). These frameworks rely on a
time scale separation of the velocity field, and a stochastic Navier-Stokes
model. In the latter, the large-scale velocity component is transported – up to
some forcings – by the small-scale and the large-scale random velocity com-
ponents. Rigorously derived from stochastic calculus theory, this stochastic
transport can be decomposed into a skew-symmetric multiplicative noise and
an eddy-viscosity-like diffusion term. Being a transport, it naturally ensures
the conservation of many physical invariants. LU and SALT models differ in
the interpretations of "transport" (classical fluid dynamics VS geometric in-
terpretation) and of "large-scale velocity". The first interpretation difference
implies in particular that the dynamics under location uncertainty conserves
kinetic energy whereas stochastic advection by Lie transport conserves helicity
filter degeneracy. Considering distance to likely scenarios rather than distance to realizations may help preventing this deficiency.

While the numerical analysis of this section 8 is applied to a random-
ized version of the SQG model with toy initial conditions, the multiplicative
structure of the noise and the balance with diffusion hold for any fluid dynam-
ics models under location uncertainty. This suggests that similar conclusions
could be expected in more complex problems. Hence, this SALT-LU method
opens for new ensemble forecasts methods, for both short-term forecasts and
climate projections where uncertainty quantification is a main issue (Allen and
Stainforth, 2002).

9 Conclusion

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

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

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

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

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

Another possible research focus is the improvement of SALT-LU parameterizations. *ADSD* method is a promising one since it is parameter-free, fast, accurate and could in theory adapt itself to any turbulent dynamics. Nevertheless, this method is currently defined in the Fourier space. Developing a physical-domain-based *ADSD* implementation would be very useful. Scaling estimations of relative diffusion (Keating et al., 2011) or velocity structure function (Gawędzki and Kupiainen, 1995) and Matérn covariances (Williams and Rasmussen, 2006; Lim and Teo, 2009; Lilly et al., 2017; Resseguier et al., 2019) could probably help in this task. An adaptation of this method to complex boundary conditions would also be needed. The Gaussian unresolved velocity moments conditionally to the unresolved velocity boundary conditions may give a path toward this direction. New spatial anisotropic and heterogeneous SALT-LU parameterizations – like *SVDpseudo* – would also be helpful. Current works based on data-driven methods (related e.g. machine learning, Koopman operator, Girsanov-based maximum likelihood estimations) are on going.

LU and SALT long-term purpose concerns data assimilation. Preminary studies have been published towards that direction. Using the pseudo-stochastic Navier-Stokes model under location uncertainty (see section 6.2.3) and a dynamics error model, Yang and Mémin (2017) have assimilated high-
and circulation. The second interpretation difference leads to a modification of the large-scale advection in LU models. After an extended theoretical description, we have reviewed existing parameterization choices – i.e. subgrid velocity statistics modeling – for the SALT-LU framework.

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

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

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

Another possible research focus is the improvement of SALT-LU parameterizations. *ADSD* method is a promising one since it is parameter-free, fast, accurate and could in theory adapt itself to any turbulent dynamics. Nevertheless, this method is currently defined in the Fourier space. Developing a physical-domain-based *ADSD* implementation would be very useful. Scaling estimations of relative diffusion (Keating et al., 2011) or velocity structure function (Gawędzki and Kupiainen, 1995) and Matérn covariances (Williams and Rasmussen, 2006; Lim and Teo, 2009; Lilly et al., 2017; Resseguier et al., 2019) could probably help in this task. An adaptation of this method to complex boundary conditions would also be needed. The Gaussian unresolved velocity moments conditionally to the unresolved velocity boundary conditions may give a path toward this direction. New spatial anisotropic and heterogeneous SALT-LU parameterizations – like *SVDpseudo* – would also be helpful. Current works based on data-driven methods (related e.g. machine learning, Koopman operator, Girsanov-based maximum likelihood estimations) are on going.

LU and SALT long-term purpose concerns data assimilation. Preminary studies have been published towards that direction. Using the pseudo-stochastic Navier-Stokes model under location uncertainty (see section 6.2.3) and a dynamics error model, Yang and Mémin (2017) have assimilated high-
resolution observations into simple shallow-water models. Yang and Mémin (2018) have applied the Navier-Stokes model under location uncertainty (see section 6.2.1) to ensemble filters. They have proposed estimations of both the covariance model error and the variance tensor. Using SALT and particle filters, Cotter et al. (2018b) assimilate data in a 2D Euler dynamics. Works on similar assimilation procedures with a quasi-geostrophic dynamics are on going. Yet, SALT-LU-based data assimilation is still in its infancy. Those first promising results need to be extended and assessed on more realistic fluid dynamic model. For flow control purposes, very fast data assimilation procedures could also be developed in stochastic reduced order frameworks (Resseguier, 2017, chapter 8).

## A Spectral energy flux in models under location uncertainty

In the homogeneous case, $\sigma d\mathbf{B}_t = \tilde{\sigma} * d\mathbf{B}_t$ and $a = a_0 \mathbb{I}_d$ does not depend on $x$. Therefore, a – possibly active – tracer $q$ is solution of:

$$0 = D_t q = d_t q + (w dt + \sigma d\mathbf{B}_t) \cdot \nabla q - \frac{a_0}{2} \Delta q dt. \quad (126)$$

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

$$d_t \hat{q}(k) = -\frac{1}{(2\pi)d} \int_{\mathbb{R}^d} (\dot{w}(k')dt + \hat{\sigma}(k')d\hat{B}_t(k')) \cdot i(k-k')\hat{q}(k-k')dk' \\ -\frac{a_0}{2} \|k\|^2 \hat{q}(k)dt. \quad (127)$$

To find a PDE associated to the Fourier modes of the tracer (Fokker-Planck equation), we need an SDE in $\mathbb{R}$. The previous one is in $\mathbb{C}$. Either we consider the closed coupled equations on the real part, $\Re(\hat{q}(k))$, and imaginary part, $\Im(\hat{q}(k))$, either we decompose trough modulus and phase. The equation on the modulus is not closed but still very instructive. From now on, the dimension of the spatial space, $d$, will be 2. However, we expect that the following results are still correct for $d = 3$. Since $\sigma d\mathbf{B}_t$ is divergence free, we set $\tilde{\sigma} = \nabla^\perp \psi_\sigma$. We can notice that:

$$d \langle \mathbf{B}_t(x_1), \mathbf{B}_t(x_2) \rangle = \delta(x_1 - x_2) dt, \quad (128)$$

and then,

$$\Re \left\{ d \left\langle \hat{B}_t(k_1), \overline{\hat{B}_t(k_2)} \right\rangle \right\} = \Re \left\{ \iint_{(\mathbb{R}^d)^2} \delta(x-y) e^{-ik_1 \cdot x} e^{+ik_2 \cdot y} dxdy \right\} dt \quad (129)$$

$$= (2\pi)^d \delta(k_1 - k_2) dt. \quad (130)$$
resolution observations into simple shallow-water models. Yang and Mémin (2018) have applied the Navier-Stokes model under location uncertainty (see section 6.2.1) to ensemble filters. They have proposed estimations of both the covariance model error and the variance tensor. Using SALT and particle filters, Cotter et al. (2018b) assimilate data in a 2D Euler dynamics. Works on similar assimilation procedures with a quasi-geostrophic dynamics are on going. Yet, SALT-LU-based data assimilation is still in its infancy. Those first promising results need to be extended and assessed on more realistic fluid dynamic model. For flow control purposes, very fast data assimilation procedures could also be developed in stochastic reduced order frameworks (Resseguier, 2017, chapter 8).

## A Spectral energy flux in models under location uncertainty

In the homogeneous case, $\sigma d\mathbf{B}_t = \tilde{\sigma} * d\mathbf{B}_t$ and $a = a_0 \mathbb{I}_d$ does not depend on $x$. Therefore, a – possibly active – tracer $q$ is solution of:

$$0 = D_t q = d_t q + (w dt + \sigma d\mathbf{B}_t) \cdot \nabla q - \frac{a_0}{2} \Delta q dt. \quad (126)$$

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

$$d_t \hat{q}(k) = -\frac{1}{(2\pi)d} \int_{\mathbb{R}^d} (\dot{w}(k')dt + \hat{\sigma}(k')d\hat{B}_t(k')) \cdot i(k-k')\hat{q}(k-k')dk' \\ -\frac{a_0}{2} \|k\|^2 \hat{q}(k)dt. \quad (127)$$

To find a PDE associated to the Fourier modes of the tracer (Fokker-Planck equation), we need an SDE in $\mathbb{R}$. The previous one is in $\mathbb{C}$. Either we consider the closed coupled equations on the real part, $\Re(\hat{q}(k))$, and imaginary part, $\Im(\hat{q}(k))$, either we decompose trough modulus and phase. The equation on the modulus is not closed but still very instructive. From now on, the dimension of the spatial space, $d$, will be 2. However, we expect that the following results are still correct for $d = 3$. Since $\sigma d\mathbf{B}_t$ is divergence free, we set $\tilde{\sigma} = \nabla^\perp \psi_\sigma$. We can notice that:

$$d \langle \mathbf{B}_t(x_1), \mathbf{B}_t(x_2) \rangle = \delta(x_1 - x_2)dt, \quad (128)$$

and then,

$$\Re \left\{ d \left\langle \hat{B}_t(k_1), \overline{\hat{B}_t(k_2)} \right\rangle \right\} = \Re \left\{ \iint_{(\mathbb{R}^d)^2} \delta(x-y) e^{-ik_1 \cdot x} e^{+ik_2 \cdot y} dxdy \right\} dt \quad (129)$$

$$= (2\pi)^d \delta(k_1 - k_2)dt. \quad (130)$$
The classical Ito formula in $\mathbb{C}$, equation (127) and the above result lead to:

$$d_t |\hat{q}(k)|^2 = 2\Re\{\hat{q}(k)d\overline{\hat{q}(k)}\} + \Re\{d\langle \hat{q}(k), \overline{\hat{q}(k)} \rangle\}, \quad (131)$$

$$= -a_0 \|k\|_2^2 \hat{q}(k) dt + \underbrace{\frac{2}{(2\pi)^d} \Re \int_{\mathbb{R}^d} \hat{w}(k') \cdot i(k-k') \hat{q}(k-k') \overline{\hat{q}(k)} dk' dt}_{=F_1}$$

$$+ \underbrace{\frac{2}{(2\pi)^d} \Re \int_{\mathbb{R}^d} \hat{\sigma}(k') d\bar{B}_t(k') \cdot i(k-k') \hat{q}(k-k') \overline{\hat{q}(k)} dk'}_{=dF_2} \quad (132)$$

$$+ \frac{(2\pi)^d}{(2\pi)^{2d}} \int_{\mathbb{R}^d} |(k-k')^T \hat{\sigma}(k') \hat{q}(k-k')|^2 dk'dt. \quad (133)$$

$dF_2$ does not modify the mean spectral energy budget since $\mathbb{E}(dF_2) = 0$, even if it modifies the law of $|\hat{q}(k)|^2$. Note that the mean spectral energy budget is the study of the evolution of $\mathbb{E}|\hat{q}(k)|^2$ whereas the spectral energy budget of the mean is the much more trivial study of $|\mathbb{E}(\hat{q}(k))|^2$. We can explicit the Laplacian term using the following expression of the variance tensor, based on the homogeneity assumption:

$$a_{ij} = a_{ij}(0) = \int_{\mathbb{R}^d} \check{\sigma}_i(z) \check{\sigma}_j(z) dz = \int_{\mathbb{R}^d} (\nabla^\perp \psi_\sigma(z))_i (\nabla^\perp \psi_\sigma(z))_j dz, \quad (134)$$

$$= \frac{1}{(2\pi)^d} \Re \int_{\mathbb{R}^d} \overline{(\nabla^\perp \psi_\sigma(k'))_i (\nabla^\perp \psi_\sigma(k'))_j} d\bar{k}', \quad (135)$$

$$= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} \underbrace{(k'^\perp)_i (k'^\perp)_j}_{\text{even function if } i = j, \text{ odd function if } i \neq j} \overline{\psi_\sigma(\|k'\|)}^2 dk', \quad (136)$$

$$= \underbrace{\frac{1}{d} \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |\|k'\| \overline{\psi_\sigma(\|k'\|)}|^2 dk'}_{=a_0} \delta_{ij}. \quad (137)$$

Injecting this into the spectral energy budget yields:

$$d_t |\hat{q}(k)|^2 = F_1 dt + dF_2 - \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |\|k\|\|k'\| \hat{\psi}_\sigma(k') \hat{q}(k)|^2 dk'dt \quad (138)$$

$$+ \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |k \cdot k'^\perp \hat{\psi}_\sigma(k') \hat{q}(k-k')|^2 dk'dt, \quad (138)$$

$$= F_1 dt + dF_2 + \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |\|k\|\|k'\| \hat{\psi}_\sigma(k')|^2 \left( \sin^2(\langle k, k' \rangle) |\hat{q}(k-k')|^2 - \frac{1}{d} |\hat{q}(k)|^2 \right) dk'dt. \quad (139)$$

Taking the expectation of the previous equation highlights the terms which contribute to the mean spectral energy budget.

$$\partial_t \mathbb{E} |\hat{q}(k)|^2 = \mathbb{E}(F_1) + \frac{1}{(2\pi)^d} \mathbb{E} \int_{\mathbb{R}^d} |\|k\|\|k'\| \hat{\psi}_\sigma(k')|^2 \left( \sin^2(\langle k, k' \rangle) |\hat{q}(k-k')|^2 - \frac{1}{d} |\hat{q}(k)|^2 \right) dk'. \quad (140)$$
The classical Ito formula in $\mathbb{C}$, equation (127) and the above result lead to:

$$d_t |\hat{q}(k)|^2 = 2 \Re \{ \hat{q}(k) d\overline{\hat{q}(k)} \} + \Re \{ d \langle \hat{q}(k), \overline{\hat{q}(k)} \rangle \}, \quad (131)$$

$$= -a_0 \|k\|^2_2 \hat{q}(k) dt + \underbrace{\frac{2}{(2\pi)^d} \Re \int_{\mathbb{R}^d} \hat{w}(k') \cdot i(k-k') \hat{q}(k-k') \overline{\hat{q}(k)} dk' dt}_{=F_1}$$

$$+ \underbrace{\frac{2}{(2\pi)^d} \Re \int_{\mathbb{R}^d} \hat{\sigma}(k') d\bar{B}_t(k') \cdot i(k-k') \hat{q}(k-k') \overline{\hat{q}(k)} dk'}_{=dF_2} \quad (132)$$

$$+ \frac{(2\pi)^d}{(2\pi)^{2d}} \int_{\mathbb{R}^d} |(k-k')^T \hat{\sigma}(k') \hat{q}(k-k')|^2 dk'dt. \quad (133)$$

$dF_2$ does not modify the mean spectral energy budget since $\mathbb{E}(dF_2) = 0$, even if it modifies the law of $|\hat{q}(k)|^2$. Note that the mean spectral energy budget is the study of the evolution of $\mathbb{E}|\hat{q}(k)|^2$ whereas the spectral energy budget of the mean is the much more trivial study of $|\mathbb{E}(\hat{q}(k))|^2$. We can explicit the Laplacian term using the following expression of the variance tensor, based on the homogeneity assumption:

$$a_{ij} = a_{ij}(0) = \int_{\mathbb{R}^d} \check{\sigma}_i(z) \check{\sigma}_j(z) dz = \int_{\mathbb{R}^d} (\nabla^\perp \psi_\sigma(z))_i (\nabla^\perp \psi_\sigma(z))_j dz, \quad (134)$$

$$= \frac{1}{(2\pi)^d} \Re \int_{\mathbb{R}^d} \overline{(\nabla^\perp \psi_\sigma(k'))_i (\nabla^\perp \psi_\sigma(k'))_j} d k', \quad (135)$$

$$= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} \underbrace{(k'^\perp)_i (k'^\perp)_j}_{\text{even function if } i = j, \text{ odd function if } i \neq j} \overline{\psi_\sigma(\|k'\|)}^2 dk', \quad (136)$$

$$= \underbrace{\frac{1}{d} \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |\|k'\| \overline{\psi_\sigma(\|k'\|)}|^2 dk'}_{=a_0} \delta_{ij}. \quad (137)$$

Injecting this into the spectral energy budget yields:

$$d_t |\hat{q}(k)|^2 = F_1 dt + dF_2 - \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |\|k\| \|k'\| \|\hat{\psi}_\sigma(k') \hat{q}(k)\|^2 dk'dt \quad (138)$$

$$+ \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |k \cdot k'^\perp \hat{\psi}_\sigma(k') \hat{q}(k-k')|^2 dk'dt, \quad (138)$$

$$= F_1 dt + dF_2 + \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |\|k\| \|k'\| \|\hat{\psi}_\sigma(k')\|^2 (\sin^2(\langle k, k' \rangle) |\hat{q}(k-k')|^2 - \frac{1}{d} |\hat{q}(k)|^2) dk'dt. \quad (139)$$

Taking the expectation of the previous equation highlights the terms which contribute to the mean spectral energy budget.

$$\partial_t \mathbb{E} |\hat{q}(k)|^2 = \mathbb{E}(F_1) + \frac{1}{(2\pi)^d} \mathbb{E} \int_{\mathbb{R}^d} |\|k\| \|k'\| \|\hat{\psi}_\sigma(k')\|^2 (\sin^2(\langle k, k' \rangle) |\hat{q}(k-k')|^2 - \frac{1}{d} |\hat{q}(k)|^2) dk'. \quad (140)$$
## B Stratonovich material derivative

Here, we express the Stratonovich material derivative:

$$D_t^S q(\mathbf{x}, t) \triangleq \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) \bigg|_{\mathbf{X}_t=\mathbf{x}} \quad (141)$$

as a function of the Ito material derivative

$$D_t q(\mathbf{x}, t) \triangleq (q(\mathbf{X}_{t+dt}, t+dt) - q(\mathbf{X}_{t}, t))|_{\mathbf{X}_t=\mathbf{x}} \quad (142)$$

and its forcing

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

To do derive these results, we will apply the Stratonovich-Ito-notation-change formula:

$$r \circ ds = rds + \frac{1}{2} d \langle r, s \rangle, \quad (144)$$

which can be found in Kunita (1997), theorem 3.2.5 page 60. First, let us rewrite equation (143) in Lagrangian coordinates:

$$d(q(\mathbf{X}_t, t)) = q(\mathbf{X}_{t+dt}, t+dt) - q(\mathbf{X}_t, t) = F(\mathbf{X}_t, t)dt + \mathbf{H}^T(\mathbf{X}_t, t)d\mathbf{B}_t. \quad (145)$$

Then, we apply the formula (144) to obtain the corresponding Stratonovich equation:

$$d(q(\mathbf{X}_t, t)) = 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), \quad (146)$$

$$= F(\mathbf{X}_t, t)dt - d\langle \mathbf{H}^T(\mathbf{X}_t, t), \mathbf{B}_t \rangle + \mathbf{H}^T(\mathbf{X}_t, t) \circ d\mathbf{B}_t, \quad (147)$$

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

since, by application of the Ito-Wentzell-formula,

$$d(H^T(X_t, t)) = (dH)^T(X_t, t) + (\sigma dB_t \cdot \nabla) H^T + (\text{others terms in } dt). \quad (149)$$

Finally, rewriting everything in the Eulerian grid $X_t = \mathbf{x}$ gives the result:

$$D_t^S q(\mathbf{x}, t) = \left( F - \frac{1}{2} \frac{d}{dt} \langle \mathbf{H}^T, \mathbf{B}_t \rangle - \frac{1}{2} \mathrm{tr} \left( (\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T \right) \right) (\mathbf{x}, t)dt + \mathbf{H}^T(\mathbf{x}, t) \circ dB_t. \quad (150)$$

After this, we can reapply formula (144) in the Eulerian grid to get:

$$D_t^S q(\mathbf{x}, t) = \left( F - \frac{1}{2} \mathrm{tr} \left( (\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T \right) \right) (\mathbf{x}, t)dt + \mathbf{H}^T(\mathbf{x}, t)dB_t, \quad (151)$$

$$= D_t q(\mathbf{x}, t) - \frac{1}{2} \mathrm{tr} \left( (\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T \right) (\mathbf{x}, t)dt. \quad (152)$$

by identification with the balance (142).

## C Effective resolution and inertial range

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

$$F(\|\mathbf{k}\|) = \exp(-\nu t \|\mathbf{k}\|^{2p}). \quad (153)$$
## B Stratonovich material derivative

Here, we express the Stratonovich material derivative:

$$D_t^S q(\mathbf{x}, t) \triangleq \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) \bigg|_{\mathbf{X}_t=\mathbf{x}}, \quad (141)$$

as a function of the Ito material derivative

$$D_t q(\mathbf{x}, t) \triangleq (q(\mathbf{X}_{t+dt}, t+dt) - q(\mathbf{X}_{t}, t))|_{\mathbf{X}_t=\mathbf{x}} \quad (142)$$

and its forcing

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

To do derive these results, we will apply the Stratonovich-Ito-notation-change formula:

$$r \circ ds = rds + \frac{1}{2} d \langle r, s \rangle, \quad (144)$$

which can be found in Kunita (1997), theorem 3.2.5 page 60. First, let us rewrite equation (143) in Lagrangian coordinates:

$$d(q(\mathbf{X}_t, t)) = q(\mathbf{X}_{t+dt}, t+dt) - q(\mathbf{X}_t, t) = F(\mathbf{X}_t, t)dt + \mathbf{H}^T(\mathbf{X}_t, t)d\mathbf{B}_t. \quad (145)$$

Then, we apply the formula (144) to obtain the corresponding Stratonovich equation:

$$d(q(\mathbf{X}_t, t)) = 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), \quad (146)$$

$$= F(\mathbf{X}_t, t)dt - d\langle \mathbf{H}^T(\mathbf{X}_t, t), \mathbf{B}_t \rangle + \mathbf{H}^T(\mathbf{X}_t, t) \circ d\mathbf{B}_t, \quad (147)$$

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

since, by application of the Ito-Wentzell-formula,

$$d(\mathbf{H}^T(\mathbf{X}_t, t)) = (d\mathbf{H})^T(\mathbf{X}_t, t) + (\boldsymbol{\sigma} d\mathbf{B}_t \cdot \nabla) \mathbf{H}^T + (\text{others terms in } dt). \quad (149)$$

Finally, rewriting everything in the Eulerian grid $\mathbf{X}_t = \mathbf{x}$ gives the result:

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

After this, we can reapply formula (144) in the Eulerian grid to get:

$$D_t^S q(\mathbf{x}, t) = \left( F - \frac{1}{2} \mathrm{tr} \left( (\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T \right) \right) (\mathbf{x}, t) dt + \mathbf{H}^T(\mathbf{x}, t) d\mathbf{B}_t, \quad (151)$$

$$= D_t q(\mathbf{x}, t) - \frac{1}{2} \mathrm{tr} \left( (\boldsymbol{\sigma}^T \nabla) \mathbf{H}^T \right) (\mathbf{x}, t) dt. \quad (152)$$

by identification with the balance (142).

## C Effective resolution and inertial range

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

$$F(\|\mathbf{k}\|) = \exp(-\nu t \|\mathbf{k}\|^{2p}). \quad (153)$$
If the hyperviscosity $\nu$ is well chosen, we may expect that at the Shanon resolution $\pi/\Delta x = \kappa_M$, only 10% of the energy is left by the filter, i.e.

$$F(\kappa_M) = 1/10. \qquad (154)$$

A ratio smaller than 10% may lead to an over-damped simulation. Moreover, the precise value of this ratio does not influence much our final estimate.

We may define the effective resolution as the scale $\kappa = \kappa_m$ where the deterministic subgrid model influence is negligible. There, we may expect the filter to be equal to 95%, i.e.:

$$F(\kappa_m) = 95/100. \qquad (155)$$

The ratio $\kappa_m/\kappa_M$ can then be derived from formulas (153), (154) and (155).

## Acknowledgments

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

## Conflict of Interest Statement

Two authors are currently employed by a private company.

## References

* Allen M, Stainforth D (2002) Towards objective probabilistic climate forecasting. Nature 419(6903):228-228

* Anderson J, Anderson S (1999) A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts. Monthly Weather Review 127(12):2741-2758

* Anderson JL (1996) A method for producing and evaluating probabilistic forecasts from ensemble model integrations. Journal of Climate 9(7):1518-1530

* Berloff P (2005) Random-forcing model of the mesoscale oceanic eddies. Journal of Fluid Mechanics 529:71-95

* Berloff P, McWilliams J (2002) Material transport in oceanic gyres. part II: Hierarchy of stochastic models. Journal of Physical Oceanography 32(3):797-830

* Berner J, Shutts G, Leutbecher M, Palmer T (2009) A spectral stochastic kinetic energy backscatter scheme and its impact on flow-dependent predictability in the ECMWF ensemble prediction system. Journal of the Atmospheric Sciences 66(3):603-626

* Berner J, Ha SY, Hacker J, Fournier A, Snyder C (2011) Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations. Monthly Weather Review 139(6):1972-1995

* Berner J, Achatz U, Batte L, Camara ADL, Crommelin D, Christensen H, Colangeli M, Dolaptchiev S, Franzke C, Friederichs P, Imkeller P, Jarvinen H, Juricke S, Kitsios V, Lott F, Lucarini V, Mahajan S, Palmer T, Penland C, Storch JSV, Sakradzija M, Weniger M, Weisheimer A, Williams P, Yano JI (2015) Stochastic parameterization: towards a new view of weather and climate models. Tech. rep., arXiv:1510.08682 [physics.ao-ph]

* Blumen W (1978) Uniform potential vorticity flow: part I. theory of wave interactions and two-dimensional turbulence. Journal of the Atmospheric Sciences 35(5):774-783
If the hyperviscosity $\nu$ is well chosen, we may expect that at the Shanon resolution $\pi/\Delta x = \kappa_M$, only 10% of the energy is left by the filter, i.e.

$$F(\kappa_M) = 1/10. \qquad (154)$$

A ratio smaller than 10% may lead to an over-damped simulation. Moreover, the precise value of this ratio does not influence much our final estimate.

We may define the effective resolution as the scale $\kappa = \kappa_m$ where the deterministic subgrid model influence is negligible. There, we may expect the filter to be equal to 95%, i.e.:

$$F(\kappa_m) = 95/100. \qquad (155)$$

The ratio $\kappa_m/\kappa_M$ can then be derived from formulas (153), (154) and (155).

## Acknowledgments

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

## Conflict of Interest Statement

Two authors are currently employed by a private company.

## References

* Allen M, Stainforth D (2002) Towards objective probabilistic climate forecasting. Nature 419(6903):228-228

* Anderson J, Anderson S (1999) A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts. Monthly Weather Review 127(12):2741-2758

* Anderson JL (1996) A method for producing and evaluating probabilistic forecasts from ensemble model integrations. Journal of Climate 9(7):1518-1530

* Berloff P (2005) Random-forcing model of the mesoscale oceanic eddies. Journal of Fluid Mechanics 529:71-95

* Berloff P, McWilliams J (2002) Material transport in oceanic gyres. part II: Hierarchy of stochastic models. Journal of Physical Oceanography 32(3):797-830

* Berner J, Shutts G, Leutbecher M, Palmer T (2009) A spectral stochastic kinetic energy backscatter scheme and its impact on flow-dependent predictability in the ECMWF ensemble prediction system. Journal of the Atmospheric Sciences 66(3):603-626

* Berner J, Ha SY, Hacker J, Fournier A, Snyder C (2011) Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations. Monthly Weather Review 139(6):1972-1995

* Berner J, Achatz U, Batte L, Camara ADL, Crommelin D, Christensen H, Colangeli M, Dolaptchiev S, Franzke C, Friederichs P, Imkeller P, Jarvinen H, Juricke S, Kitsios V, Lott F, Lucarini V, Mahajan S, Palmer T, Penland C, Storch JSV, Sakradzija M, Weniger M, Weisheimer A, Williams P, Yano JI (2015) Stochastic parameterization: towards a new view of weather and climate models. Tech. rep., arXiv:1510.08682 [physics.ao-ph]

* Blumen W (1978) Uniform potential vorticity flow: part I. theory of wave interactions and two-dimensional turbulence. Journal of the Atmospheric Sciences 35(5):774-783
Bocquet M, Sakov P (2014) An iterative ensemble kalman smoother. Quarterly Journal of the Royal Meteorological Society 140(682):1521–1535

Boussinesq J (1877) Essai sur la théorie des eaux courantes. Mémoires présentés par divers savants à l'Académie des Sciences, 23 (1): 1–680

Brankart J (2013) Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modeling. Ocean Modelling 66:64–76

Brzeźniak Z, Capiński M, Flandoli F (1991) Stochastic partial differential equations and turbulence. Mathematical Models and Methods in Applied Sciences 1(01):41–59

Buehner M (2005) Ensemble-derived stationary and flow-dependent background-error covariance: Evaluation in a quasi-operational nwp setting. Quarterly Journal of the Royal Meteorological Society 131(607):1013–1043

Buizza R (2016a) Representing model uncertainty for climate forecasts. In: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading

Buizza R (2016b) Weather prediction in a world of uncertainties: should ensembles simulate the effect of model approximations? In: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading

Buizza R, Miller M, Palmer T (1999) Stochastic representation of model uncertainties in the ECMWF ensemble prediction system. Quarterly Journal Royal Meteorological Society 125:2887–2908

Cai S, Mémin E, Dérian P, Xu C (2018) Motion estimation under location uncertainty for turbulent fluid flows. Experiments in Fluids 59(1):8

Candy J (2011) Bayesian signal processing: Classical, modern and particle filtering methods, vol 54. John Wiley & Sons

Capet X, Klein P, Hua B, Lapeyre G, McWilliams J (2008) Surface kinetic energy transfer in surface quasi-geostrophic flows. J Fluid Mech 604:165–174

Chandramouli P, Heitz D, Laizet S, Mémin E (2018) Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty. Computers & Fluids 168:170–189

Chandramouli P, Mémin E, Chapron B, nd D Heitz SL (2019) Deciphering the role of small-scale inhomogeneity on flow structuration: a stochastic approach, submitted

Chapron B, Dérian P, Mémin E, Resseguier V (2018) Large-scale flows under location uncertainty: a consistent stochastic framework. Quarterly Journal of the Royal Meteorological Society 144(710):251–260

Chasnov J (1991) Simulation of the Kolmogorov inertial subrange using an improved subgrid model. Physics of Fluids A: Fluid Dynamics (1989-1993) 3(1):188–200

Chekroun MD, Kondrashov D, Ghil M (2011) Predicting stochastic systems by noise sampling, and application to the el niño-southern oscillation. Proceedings of the National Academy of Sciences 108(29):11766–11771

Chekroun MD, Liu H, McWilliams JC (2017) The emergence of fast oscillations in a reduced primitive equation model and its implications for closure theories. Computers & Fluids 151:3–22

Chow PL (2014) Stochastic partial differential equations. CRC Press

Cintolesi C, Mémin E (2019a) Pseudo-stochastic simulation of turbulent channel flows with near-wall modelling, URL https://hal.archives-ouvertes.fr/hal-02044818, working paper or preprint

Cintolesi C, Mémin E (2019b) Stochastic modelling of turbulent flows for numerical simulations, URL https://hal.archives-ouvertes.fr/hal-02044809

Cotter C, Gottwald G, Holm D (2017) Stochastic partial differential fluid equations as a diffusive limit of deterministic lagrangian multi-time dynamics. arXiv preprint arXiv:170600287

Cotter C, Crisan D, Holm DD, Pan W, Shevchenko I (2018a) Modelling uncertainty using circulation-preserving stochastic transport noise in a 2-layer quasi-geostrophic model. arXiv preprint arXiv:180205711

Cotter C, Crisan D, Holm DD, Pan W, Shevchenko I (2018b) Sequential monte carlo for stochastic advection by lie transport. In preparation

Cotter CJ, Crisan D, Holm DD, Pan W, Shevchenko I (2018c) Numerically modelling stochastic lie transport in fluid dynamics. arXiv preprint arXiv:180109729
Bocquet M, Sakov P (2014) An iterative ensemble kalman smoother. Quarterly Journal of the Royal Meteorological Society 140(682):1521–1535

Boussinesq J (1877) Essai sur la théorie des eaux courantes. Mémoires présentés par divers savants à l'Académie des Sciences, 23 (1): 1–680

Brankart J (2013) Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modeling. Ocean Modelling 66:64–76

Brzeźniak Z, Capiński M, Flandoli F (1991) Stochastic partial differential equations and turbulence. Mathematical Models and Methods in Applied Sciences 1(01):41–59

Buehner M (2005) Ensemble-derived stationary and flow-dependent background-error covariance: Evaluation in a quasi-operational nwp setting. Quarterly Journal of the Royal Meteorological Society 131(607):1013–1043

Buizza R (2016a) Representing model uncertainty for climate forecasts. In: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading

Buizza R (2016b) Weather prediction in a world of uncertainties: should ensembles simulate the effect of model approximations? In: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading

Buizza R, Miller M, Palmer T (1999) Stochastic representation of model uncertainties in the ECMWF ensemble prediction system. Quarterly Journal Royal Meteorological Society 125:2887–2908

Cai S, Mémin E, Dérian P, Xu C (2018) Motion estimation under location uncertainty for turbulent fluid flows. Experiments in Fluids 59(1):8

Candy J (2011) Bayesian signal processing: Classical, modern and particle filtering methods, vol 54. John Wiley & Sons

Capet X, Klein P, Hua B, Lapeyre G, McWilliams J (2008) Surface kinetic energy transfer in surface quasi-geostrophic flows. J Fluid Mech 604:165–174

Chandramouli P, Heitz D, Laizet S, Mémin E (2018) Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty. Computers & Fluids 168:170–189

Chandramouli P, Mémin E, Chapron B, nd D Heitz SL (2019) Deciphering the role of small-scale inhomogeneity on flow structuration: a stochastic approach, submitted

Chapron B, Dérian P, Mémin E, Resseguier V (2018) Large-scale flows under location uncertainty: a consistent stochastic framework. Quarterly Journal of the Royal Meteorological Society 144(710):251–260

Chasnov J (1991) Simulation of the Kolmogorov inertial subrange using an improved subgrid model. Physics of Fluids A: Fluid Dynamics (1989-1993) 3(1):188–200

Chekroun MD, Kondrashov D, Ghil M (2011) Predicting stochastic systems by noise sampling, and application to the el niño-southern oscillation. Proceedings of the National Academy of Sciences 108(29):11766–11771

Chekroun MD, Liu H, McWilliams JC (2017) The emergence of fast oscillations in a reduced primitive equation model and its implications for closure theories. Computers & Fluids 151:3–22

Chow PL (2014) Stochastic partial differential equations. CRC Press

Cintolesi C, Mémin E (2019a) Pseudo-stochastic simulation of turbulent channel flows with near-wall modelling, URL https://hal.archives-ouvertes.fr/hal-02044818, working paper or preprint

Cintolesi C, Mémin E (2019b) Stochastic modelling of turbulent flows for numerical simulations, URL https://hal.archives-ouvertes.fr/hal-02044809

Cotter C, Gottwald G, Holm D (2017) Stochastic partial differential fluid equations as a diffusive limit of deterministic lagrangian multi-time dynamics. arXiv preprint arXiv:170600287

Cotter C, Crisan D, Holm DD, Pan W, Shevchenko I (2018a) Modelling uncertainty using circulation-preserving stochastic transport noise in a 2-layer quasi-geostrophic model. arXiv preprint arXiv:180205711

Cotter C, Crisan D, Holm DD, Pan W, Shevchenko I (2018b) Sequential monte carlo for stochastic advection by lie transport. In preparation

Cotter CJ, Crisan D, Holm DD, Pan W, Shevchenko I (2018c) Numerically modelling stochastic lie transport in fluid dynamics. arXiv preprint arXiv:180109729
Craig G, Forbes R, Abdalla S, Balsamo G, Bechtold P, Berner J, Buizza R, Pallares AC, Meutter PD, Dueben P, Frogner IL, Gagnon N, Hodyss D, Holm D, Lock SJ, Machulskaya E, Najafi H, Ollinaho P, Selz T, Separovic L, Subramanian A, Tseng J, Weisheimer A (2016) Reports from working group 1: What are the sources of model error and how can we improve the physical basis of model uncertainty representation? In: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading

Crisan D, Flandoli F, Holm DD (2017) Solution properties of a 3d stochastic euler fluid equation. arXiv preprint arXiv:170406989

Da Prato G, Zabczyk J (1992) Stochastic Equations in Infinite Dimensions. Encyclopedia of Mathematics and its Applications, Cambridge University Press

Desroziers G, Berre L, Chapnik B, Poli P (2005) Diagnosis of observation, background and analysis-error statistics in observation space. Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 131(613):3385-3396

Doucet A, Johansen A (2009) A tutorial on particle filtering and smoothing: Fifteen years later. Handbook of Nonlinear Filtering 12:656-704

Doucet A, De Freitas N, Gordon N (2001) Sequential Monte Carlo methods in practice. Springer

Dukowicz JK, Smith RD (1997) Stochastic theory of compressible turbulent fluid transport. Physics of Fluids 9(11):3523-3529

Dwivedi S, Franzke CL, Lunkeit F (2019) Energetically consistent scale adaptive stochastic and deterministic energy backscatter schemes for an atmospheric model. Quarterly Journal of the Royal Meteorological Society

Ehrendorfer M (1997) Predicting the uncertainty of numerical weather forecasts: A review. METEOROLOGISCHE ZEITSCHRIFT-BERLIN- 6:147-183

Falkovich G, Gawędzki K, Vergassola M (2001) Particles and fields in fluid turbulence. Reviews of modern Physics 73(4):913

Farrell B, Ioannou P (2014) Statistical state dynamics: a new perspective on turbulence in shear flow. arXiv preprint arXiv:14128290

Ferrari R, Nikurashin M (2010) Suppression of eddy diffusivity across jets in the Southern Ocean. Journal of Physical Oceanography 40(7):1501-1519

Flandoli F (2011) The interaction between noise and transport mechanisms in PDEs. Milan Journal of Mathematics 79(2):543-560

Frank JE, Gottwald GA (2013) Stochastic homogenization for an energy conserving multi-scale toy model of the atmosphere. Physica D: Nonlinear Phenomena 254:46-56

Franzke C, Majda A, Vanden-Eijnden E (2005) Low-order stochastic mode reduction for a realistic barotropic model climate. Journal of the atmospheric sciences 62(6):1722-1745

Franzke C, O'Kane T, Berner J, Williams P, Lucarini V (2015) Stochastic climate theory and modeling. Wiley Interdisciplinary Reviews: Climate Change 6(1):63-78

Franzke CL (2017) Extremes in dynamic-stochastic systems. Chaos: An Interdisciplinary Journal of Nonlinear Science 27(1):012101

Frederiksen JS (1999) Subgrid-scale parameterizations of eddy-topographic force, eddy viscosity, and stochastic backscatter for flow over topography. Journal of the atmospheric sciences 56(11):1481-1494

Frisch U (1995) Turbulence: the legacy of AN Kolmogorov. Cambridge university press

Gawędzki K, Kupiainen A (1995) Anomalous scaling of the passive scalar. Physical review letters 75(21):3834

Gay-Balmaz F, Holm DD (2018) Stochastic geometric models with non-stationary spatial correlations in lagrangian fluid flows. Journal of nonlinear science 28(3):873-904

Gent P, Mcwilliams J (1990) Isopycnal mixing in ocean circulation models. Journal of Physical Oceanography 20(1):150-155

Germano M, Piomelli U, Moin P, Cabot W (1991) A dynamic subgrid-scale eddy viscosity model. Phys of Fluids 3:1760-1765

Givon D, Kupferman R, Stuart A (2004) Extracting macroscopic dynamics: model problems and algorithms. Nonlinearity 17(6):R55

Gneiting T, Raftery AE (2007) Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association 102(477):359-378
Craig G, Forbes R, Abdalla S, Balsamo G, Bechtold P, Berner J, Buizza R, Pallares AC, Meutter PD, Dueben P, Frogner IL, Gagnon N, Hodyss D, Holm D, Lock SJ, Machulskaya E, Najafi H, Ollinaho P, Selz T, Separovic L, Subramanian A, Tseng J, Weisheimer A (2016) Reports from working group 1: What are the sources of model error and how can we improve the physical basis of model uncertainty representation? In: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading

Crisan D, Flandoli F, Holm DD (2017) Solution properties of a 3d stochastic euler fluid equation. arXiv preprint arXiv:170406989

Da Prato G, Zabczyk J (1992) Stochastic Equations in Infinite Dimensions. Encyclopedia of Mathematics and its Applications, Cambridge University Press

Desroziers G, Berre L, Chapnik B, Poli P (2005) Diagnosis of observation, background and analysis-error statistics in observation space. Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 131(613):3385-3396

Doucet A, Johansen A (2009) A tutorial on particle filtering and smoothing: Fifteen years later. Handbook of Nonlinear Filtering 12:656-704

Doucet A, De Freitas N, Gordon N (2001) Sequential Monte Carlo methods in practice. Springer

Dukowicz JK, Smith RD (1997) Stochastic theory of compressible turbulent fluid transport. Physics of Fluids 9(11):3523-3529

Dwivedi S, Franzke CL, Lunkeit F (2019) Energetically consistent scale adaptive stochastic and deterministic energy backscatter schemes for an atmospheric model. Quarterly Journal of the Royal Meteorological Society

Ehrendorfer M (1997) Predicting the uncertainty of numerical weather forecasts: A review. METEOROLOGISCHE ZEITSCHRIFT-BERLIN- 6:147-183

Falkovich G, Gawędzki K, Vergassola M (2001) Particles and fields in fluid turbulence. Reviews of modern Physics 73(4):913

Farrell B, Ioannou P (2014) Statistical state dynamics: a new perspective on turbulence in shear flow. arXiv preprint arXiv:14128290

Ferrari R, Nikurashin M (2010) Suppression of eddy diffusivity across jets in the Southern Ocean. Journal of Physical Oceanography 40(7):1501-1519

Flandoli F (2011) The interaction between noise and transport mechanisms in PDEs. Milan Journal of Mathematics 79(2):543-560

Frank JE, Gottwald GA (2013) Stochastic homogenization for an energy conserving multi-scale toy model of the atmosphere. Physica D: Nonlinear Phenomena 254:46-56

Franzke C, Majda A, Vanden-Eijnden E (2005) Low-order stochastic mode reduction for a realistic barotropic model climate. Journal of the atmospheric sciences 62(6):1722-1745

Franzke C, O'Kane T, Berner J, Williams P, Lucarini V (2015) Stochastic climate theory and modeling. Wiley Interdisciplinary Reviews: Climate Change 6(1):63-78

Franzke CL (2017) Extremes in dynamic-stochastic systems. Chaos: An Interdisciplinary Journal of Nonlinear Science 27(1):012101

Frederiksen JS (1999) Subgrid-scale parameterizations of eddy-topographic force, eddy viscosity, and stochastic backscatter for flow over topography. Journal of the atmospheric sciences 56(11):1481-1494

Frisch U (1995) Turbulence: the legacy of AN Kolmogorov. Cambridge university press

Gawędzki K, Kupiainen A (1995) Anomalous scaling of the passive scalar. Physical review letters 75(21):3834

Gay-Balmaz F, Holm DD (2018) Stochastic geometric models with non-stationary spatial correlations in lagrangian fluid flows. Journal of nonlinear science 28(3):873-904

Gent P, Mcwilliams J (1990) Isopycnal mixing in ocean circulation models. Journal of Physical Oceanography 20(1):150-155

Germano M, Piomelli U, Moin P, Cabot W (1991) A dynamic subgrid-scale eddy viscosity model. Phys of Fluids 3:1760-1765

Givon D, Kupferman R, Stuart A (2004) Extracting macroscopic dynamics: model problems and algorithms. Nonlinearity 17(6):R55

Gneiting T, Raftery AE (2007) Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association 102(477):359-378
Gneiting T, Raftery AE, Westveld III AH, Goldman T (2005) Calibrated probabilistic forecasting using ensemble model output statistics and minimum crps estimation. Monthly Weather Review 133(5):1098-1118

Gottwald G, Harlim J (2013) The role of additive and multiplicative noise in filtering complex dynamical systems. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 469(2155):20130096

Gottwald G, Melbourne I (2013) Homogenization for deterministic maps and multiplicative noise. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 469(2156)

Gottwald G, Crommelin D, Franzke C (2015) Stochastic climate theory. In: Nonlinear and Stochastic Climate Dynamics, Cambridge University Press

Grimit EP, Gneiting T, Berrocal V, Johnson NA (2006) The continuous ranked probability score for circular variables and its application to mesoscale forecast ensemble verification. Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 132(621C):2925-2942

Grooms I, Majda A (2014) Stochastic superparameterization in quasigeostrophic turbulence. Journal of Computational Physics 271:78-98

Gugole F, Franzke C (219) Numerical development and evaluation of an energy conserving conceptual stochastic climate model. Mathematics of climate and weather forecasting In press

Gula J, Molemaker J, McWilliams J (2015) Gulf stream dynamics along the southeastern us seaboard. Journal of Physical Oceanography 45(3):690-715

Hamill TM (2001) Interpretation of rank histograms for verifying ensemble forecasts. Monthly Weather Review 129(3):550-560

Hamill TM, Colucci SJ (1997) Verification of eta-rsm short-range ensemble forecasts. Monthly Weather Review 125(6):1312-1327

Harouna SK, Mémin E (2017) Stochastic representation of the reynolds transport theorem: revisiting large-scale modeling. Computers & Fluids 156:456-469

Hasselmann K (1976) Stochastic climate models. part I: theory. Tellus 28:473-485

Held I, Pierrehumbert R, Garner S, Swanson K (1995) Surface quasi-geostrophic dynamics. Journal of Fluid Mechanics 282:1-20

Hersbach H (2000) Decomposition of the continuous ranked probability score for ensemble prediction systems. Weather and Forecasting 15(5):559-570

Holm D (2015) Variational principles for stochastic fluid dynamics. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 471(2176)

Holmes P, Lumley J, Berkooz G (1996) Turbulence, coherence structures, dynamical systems and symmetry. Cambridge university press

Jain A, Timofeyev I, Vanden-Eijnden E (2014) Stochastic mode-reduction in models with conservative fast sub-systems. arXiv preprint arXiv:14103004

Jansen MF, Held IM (2014) Parameterizing subgrid-scale eddy effects using energetically consistent backscatter. Ocean Modelling 80:36-48

Jordan A (2016) Facets of forecast evaluation

Keating S, Smith S, Kramer P (2011) Diagnosing lateral mixing in the upper ocean with virtual tracers: Spatial and temporal resolution dependence. Journal of Physical Oceanography 41(8):1512-1534

Keating S, Majda A, Smith S (2012) New methods for estimating ocean eddy heat transport using satellite altimetry. Monthly Weather Review 140(5):1703-1722

Klein P, Hua B, Lapeyre G, Capet X, Le Gentil S, Sasaki H (2008) Upper ocean turbulence from high-resolution 3D simulations. Journal of Physical Oceanography 38(8):1748-1763

Klyatskin V (2005) Stochastic equations through the eye of the physicist: Basic concepts, exact results and asymptotic approximations. Elsevier

Klyatskin V, Woyczynski W, Gurarie D (1996) Short-time correlation approximations for diffusing tracers in random velocity fields: A functional approach. In: Stochastic modelling in physical oceanography, Springer, pp 221-269

Kondrashov D, Berloff P (2015) Stochastic modeling of decadal variability in ocean gyres. Geophysical Research Letters 42(5):1543-1553

Kondrashov D, Chekroun MD, Ghil M (2015) Data-driven non-markovian closure models. Physica D: Nonlinear Phenomena 297:33-55
Gneiting T, Raftery AE, Westveld III AH, Goldman T (2005) Calibrated probabilistic forecasting using ensemble model output statistics and minimum crps estimation. Monthly Weather Review 133(5):1098-1118

Gottwald G, Harlim J (2013) The role of additive and multiplicative noise in filtering complex dynamical systems. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 469(2155):20130096

Gottwald G, Melbourne I (2013) Homogenization for deterministic maps and multiplicative noise. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 469(2156)

Gottwald G, Crommelin D, Franzke C (2015) Stochastic climate theory. In: Nonlinear and Stochastic Climate Dynamics, Cambridge University Press

Grimit EP, Gneiting T, Berrocal V, Johnson NA (2006) The continuous ranked probability score for circular variables and its application to mesoscale forecast ensemble verification. Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 132(621C):2925-2942

Grooms I, Majda A (2014) Stochastic superparameterization in quasigeostrophic turbulence. Journal of Computational Physics 271:78-98

Gugole F, Franzke C (219) Numerical development and evaluation of an energy conserving conceptual stochastic climate model. Mathematics of climate and weather forecasting In press

Gula J, Molemaker J, McWilliams J (2015) Gulf stream dynamics along the southeastern us seaboard. Journal of Physical Oceanography 45(3):690-715

Hamill TM (2001) Interpretation of rank histograms for verifying ensemble forecasts. Monthly Weather Review 129(3):550-560

Hamill TM, Colucci SJ (1997) Verification of eta-rsm short-range ensemble forecasts. Monthly Weather Review 125(6):1312-1327

Harouna SK, Mémin E (2017) Stochastic representation of the reynolds transport theorem: revisiting large-scale modeling. Computers & Fluids 156:456-469

Hasselmann K (1976) Stochastic climate models. part I: theory. Tellus 28:473-485

Held I, Pierrehumbert R, Garner S, Swanson K (1995) Surface quasi-geostrophic dynamics. Journal of Fluid Mechanics 282:1-20

Hersbach H (2000) Decomposition of the continuous ranked probability score for ensemble prediction systems. Weather and Forecasting 15(5):559-570

Holm D (2015) Variational principles for stochastic fluid dynamics. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 471(2176)

Holmes P, Lumley J, Berkooz G (1996) Turbulence, coherence structures, dynamical systems and symmetry. Cambridge university press

Jain A, Timofeyev I, Vanden-Eijnden E (2014) Stochastic mode-reduction in models with conservative fast sub-systems. arXiv preprint arXiv:14103004

Jansen MF, Held IM (2014) Parameterizing subgrid-scale eddy effects using energetically consistent backscatter. Ocean Modelling 80:36-48

Jordan A (2016) Facets of forecast evaluation

Keating S, Smith S, Kramer P (2011) Diagnosing lateral mixing in the upper ocean with virtual tracers: Spatial and temporal resolution dependence. Journal of Physical Oceanography 41(8):1512-1534

Keating S, Majda A, Smith S (2012) New methods for estimating ocean eddy heat transport using satellite altimetry. Monthly Weather Review 140(5):1703-1722

Klein P, Hua B, Lapeyre G, Capet X, Le Gentil S, Sasaki H (2008) Upper ocean turbulence from high-resolution 3D simulations. Journal of Physical Oceanography 38(8):1748-1763

Klyatskin V (2005) Stochastic equations through the eye of the physicist: Basic concepts, exact results and asymptotic approximations. Elsevier

Klyatskin V, Woyczynski W, Gurarie D (1996) Short-time correlation approximations for diffusing tracers in random velocity fields: A functional approach. In: Stochastic modelling in physical oceanography, Springer, pp 221-269

Kondrashov D, Berloff P (2015) Stochastic modeling of decadal variability in ocean gyres. Geophysical Research Letters 42(5):1543-1553

Kondrashov D, Chekroun MD, Ghil M (2015) Data-driven non-markovian closure models. Physica D: Nonlinear Phenomena 297:33-55
Kraichnan R (1968) Small-scale structure of a scalar field convected by turbulence. Physics of Fluids (1958-1988) 11(5):945-953

Kraichnan R (1994) Anomalous scaling of a randomly advected passive scalar. Physical review letters 72(7):1016

Kunita H (1997) Stochastic flows and stochastic differential equations, vol 24. Cambridge university press

Kurtz T (1973) A limit theorem for perturbed operator semigroups with applications to random evolutions. Journal of Functional Analysis 12(1):55-67

Laio F, Tamea S (2007) Verification tools for probabilistic forecasts of continuous hydrological variables. Hydrology and Earth System Sciences 11(4):1267-1277

Lapeyre G (2017) Surface quasi-geostrophy. Fluids 2(1):7

Le Dimet FX, Talagrand O (1986) Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. Tellus A 38(2):97-110

Le Maitre O, Reagan M, Najm H, Ghanem R, Knio O (2002) A stochastic projection method for fluid flow. II. random process. Journal of Computational Physics 181:9-44

Leith C (1971) Atmospheric predictability and two-dimensional turbulence. Journal of the Atmospheric Sciences 28(2):145-161

Leith C (1990) Stochastic backscatter in a subgrid-scale model: Plane shear mixing layer. Physics of Fluids A: Fluid Dynamics (1989-1993) 2(3):297-299

Leutbechner M, Ollinaha P, Lock SJ, Lang S, Bechtold P, Beljaars A, Bozzo A, Forbes R, Haiden R T Hogan, Sandu I (2016) Stochastic representations of model uncertainties in the ifs. In: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading

Lilly D (1966) On the application of the eddy viscosity concept in the inertial subrange of turbulence. Tech. Rep. 123, NCAR

Lilly J, Sykulski A, Early J, Olhede S (2017) Fractional brownian motion, the Matérn process, and stochastic modeling of turbulent dispersion. Nonlinear Processes in Geophysics 24(1):481-514

Lim S, Teo L (2009) Generalized whittle-matérn random field as a model of correlated fluctuations. Journal of Physics A: Mathematical and Theoretical 42(10):105202

Liu C, Xiao Q, Wang B (2009) An ensemble-based four-dimensional variational data assimilation scheme. part ii: Observing system simulation experiments with advanced research wrf (arw). Monthly Weather Review 137(5):1687-1704

Lorenz E (1969) The predictability of a flow which possesses many scales of motion. Tellus 21:289-307

Lu F, Lin KK, Chorin AJ (2017) Data-based stochastic model reduction for the kuramoto-sivashinsky equation. Physica D: Nonlinear Phenomena 340:46-57

Lucarini V, Blender R, Herbert C, Ragone F, Pascale S, Wouters J (2014) Mathematical and physical ideas for climate science. Reviews of Geophysics 52(4):809-859

Majda A (2015) Statistical energy conservation principle for inhomogeneous turbulent dynamical systems. Proceedings of the National Academy of Sciences 112(29):8937-8941

Majda A, Timofeyev I, Eijnden E (1999) Models for stochastic climate prediction. Proceedings of the National Academy of Sciences 96(26):14687-14691

Majda A, Timofeyev I, Vanden Eijnden E (2001) A mathematical framework for stochastic climate models. Communications on Pure and Applied Mathematics 54(8):891-974

Majda A, Franzke C, Khouider B (2008) An applied mathematics perspective on stochastic modelling for climate. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 366(1875):2427-2453

Mana P, Zanna L (2014) Toward a stochastic parameterization of ocean mesoscale eddies. Ocean Modelling 79:1-20

Matheson JE, Winkler RL (1976) Scoring rules for continuous probability distributions. Management science 22(10):1087-1096

McCabe A, Swinbank R, Tennant W, Lock A (2016) Representing model error in the met office convection permitting ensemble prediction system. In: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading

Mémin E (2014) Fluid flow dynamics under location uncertainty. Geophysical & Astrophysical Fluid Dynamics 108(2):119-146, DOI 10.1080/03091929.2013.836190

Métais O, Lesieur M (1986) Statistical predictability of decaying turbulence. Journal of the atmospheric sciences 43(9):857-870
Kraichnan R (1968) Small-scale structure of a scalar field convected by turbulence. Physics of Fluids (1958-1988) 11(5):945-953

Kraichnan R (1994) Anomalous scaling of a randomly advected passive scalar. Physical review letters 72(7):1016

Kunita H (1997) Stochastic flows and stochastic differential equations, vol 24. Cambridge university press

Kurtz T (1973) A limit theorem for perturbed operator semigroups with applications to random evolutions. Journal of Functional Analysis 12(1):55-67

Laio F, Tamea S (2007) Verification tools for probabilistic forecasts of continuous hydrological variables. Hydrology and Earth System Sciences 11(4):1267-1277

Lapeyre G (2017) Surface quasi-geostrophy. Fluids 2(1):7

Le Dimet FX, Talagrand O (1986) Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. Tellus A 38(2):97-110

Le Maitre O, Reagan M, Najm H, Ghanem R, Knio O (2002) A stochastic projection method for fluid flow. II. random process. Journal of Computational Physics 181:9-44

Leith C (1971) Atmospheric predictability and two-dimensional turbulence. Journal of the Atmospheric Sciences 28(2):145-161

Leith C (1990) Stochastic backscatter in a subgrid-scale model: Plane shear mixing layer. Physics of Fluids A: Fluid Dynamics (1989-1993) 2(3):297-299

Leutbechner M, Ollinaha P, Lock SJ, Lang S, Bechtold P, Beljaars A, Bozzo A, Forbes R, Haiden R T Hogan, Sandu I (2016) Stochastic representations of model uncertainties in the ifs. In: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading

Lilly D (1966) On the application of the eddy viscosity concept in the inertial subrange of turbulence. Tech. Rep. 123, NCAR

Lilly J, Sykulski A, Early J, Olhede S (2017) Fractional brownian motion, the Matérn process, and stochastic modeling of turbulent dispersion. Nonlinear Processes in Geophysics 24(1):481-514

Lim S, Teo L (2009) Generalized whittle-matérn random field as a model of correlated fluctuations. Journal of Physics A: Mathematical and Theoretical 42(10):105202

Liu C, Xiao Q, Wang B (2009) An ensemble-based four-dimensional variational data assimilation scheme. part ii: Observing system simulation experiments with advanced research wrf (arw). Monthly Weather Review 137(5):1687-1704

Lorenz E (1969) The predictability of a flow which possesses many scales of motion. Tellus 21:289-307

Lu F, Lin KK, Chorin AJ (2017) Data-based stochastic model reduction for the kuramoto-sivashinsky equation. Physica D: Nonlinear Phenomena 340:46-57

Lucarini V, Blender R, Herbert C, Ragone F, Pascale S, Wouters J (2014) Mathematical and physical ideas for climate science. Reviews of Geophysics 52(4):809-859

Majda A (2015) Statistical energy conservation principle for inhomogeneous turbulent dynamical systems. Proceedings of the National Academy of Sciences 112(29):8937-8941

Majda A, Timofeyev I, Eijnden E (1999) Models for stochastic climate prediction. Proceedings of the National Academy of Sciences 96(26):14687-14691

Majda A, Timofeyev I, Vanden Eijnden E (2001) A mathematical framework for stochastic climate models. Communications on Pure and Applied Mathematics 54(8):891-974

Majda A, Franzke C, Khouider B (2008) An applied mathematics perspective on stochastic modelling for climate. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 366(1875):2427-2453

Mana P, Zanna L (2014) Toward a stochastic parameterization of ocean mesoscale eddies. Ocean Modelling 79:1-20

Matheson JE, Winkler RL (1976) Scoring rules for continuous probability distributions. Management science 22(10):1087-1096

McCabe A, Swinbank R, Tennant W, Lock A (2016) Representing model error in the met office convection permitting ensemble prediction system. In: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading

Mémin E (2014) Fluid flow dynamics under location uncertainty. Geophysical & Astrophysical Fluid Dynamics 108(2):119-146, DOI 10.1080/03091929.2013.836190

Métais O, Lesieur M (1986) Statistical predictability of decaying turbulence. Journal of the atmospheric sciences 43(9):857-870