# Data-driven versus self-similar parameterizations for Stochastic Advection by Lie Transport and Location Uncertainty

Valentin Resseguier, Wei Pan, Baylor Fox-Kemper

► To cite this version:

Valentin Resseguier, Wei Pan, Baylor Fox-Kemper. Data-driven versus self-similar parameterizations for Stochastic Advection by Lie Transport and Location Uncertainty. Nonlinear Processes in Geo-physics, 2020, 27 (2), pp.209-234. 10.5194/npg-27-209-2020. hal-02371789v3

HAL Id: hal-02371789

https://hal.science/hal-02371789v3

Submitted on 29 May 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Data-driven versus self-similar parameterizations for stochastic advection by Lie transport and location uncertainty

Valentin Resseguier¹, Wei Pan², and Baylor Fox-Kemper³,⁴

¹Lab, SCALIAN DS, Espace Nobel, 2 Allée de Becquerel, 35700 Rennes, France

²Department of Mathematics, Imperial College London, London SW7 2AZ, UK

³Department of Earth, Environmental and Planetary Sciences (DEEPS), Brown University, Providence, RI 02912, USA

⁴Institute at Brown for Environment and Society (IBES), Brown University, Providence, RI 02912, USA

**Correspondence:** Valentin Resseguier (valentin.resseguier@scalian.com)

Received: 13 October 2019 – Discussion started: 16 October 2019

Revised: 23 January 2020 – Accepted: 9 February 2020 – Published: 16 April 2020

**Abstract.** Stochastic subgrid parameterizations enable en-
semble forecasts of fluid dynamic systems and ultimately
accurate data assimilation (DA). Stochastic advection by
Lie transport (SALT) and models under location uncer-
tainty (LU) are recent and similar physically based stochas-
tic schemes. SALT dynamics conserve helicity, whereas
LU models conserve kinetic energy (KE). After highlight-
ing general similarities between LU and SALT frameworks,
this paper focuses on their common challenge: the param-
eterization choice. We compare uncertainty quantification
skills of a stationary heterogeneous data-driven parameter-
ization and a non-stationary homogeneous self-similar pa-
rameterization. For stationary, homogeneous surface quasi-
geostrophic (SQG; QG) turbulence, both parameterizations
lead to high-quality ensemble forecasts. This paper also
discusses a heterogeneous adaptation of the homogeneous
parameterization targeted at a better simulation of strong
straight buoyancy fronts.

# 1 Introduction

Geophysical fluid dynamics and other turbulent flows cover
a wide range of spatial and temporal scales (see, e.g. Fox-
Kemper, 2018). However, numerical simulations and obser-
vations are limited in the range of scales that can be simu-
lated or observed. Accordingly, deterministic subgrid models
or parameterizations (e.g. Margolin et al., 2006; San et al.,
2013; Bachman et al., 2017; Voosen, 2018) and regulariza-
tions (e.g. nugget in optimal interpolation, Tandeo et al.,

2014) are incorporated into numerical simulations and mea-
surement processes. The limited resolution of these numer-
ical simulations introduces unknown errors which are con-
tinuously amplified during the course of the simulation and
which compound with the measurement errors to the extent
that the simulations model the true chaotic nature of fluid
dynamics. Improvements in accuracy of these simulations
can sometimes be made by judicious use of data assimilation
(DA) techniques. In the data assimilation process, observa-
tions (data) are integrated into the model state of the numer-
ical simulation limited by imperfect physics and imperfect
initial conditions. DA proceeds by alternating between fore-
cast and analysis cycles. In each analysis cycle, observations
are combined with the results from a prediction model (the
forecast) to produce a so-called analysis. The analysis is ex-
pected to be more accurate than either predictions based on
physical models without incorporating observations or pre-
dictions from physics-free interpolations of the observations.
In order to optimize the analysis cycle, confidence in both
the observations and the prediction needs to be quantified.
Observational accuracy depends on the instrumental preci-
sion and sampling; and well-known uncertainty quantifica-
tion (UQ) techniques are available to estimate the uncertainty
associated with these limitations. In contrast, the uncertainty
due to imperfect model physics is more difficult to quan-
tify because often a “perfect model”, exact solution, or di-
rect numerical simulation is not available or feasible for the
intended application.

Historically introduced to address the limited range of
scales simulated (e.g. Orszag, 1970; Leith, 1971), stochas-
tic subgrid parameterization is now a widely used UQ tool for assessing the impact of imperfect physics in numerical simulations. In particular, it is generally thought that building such stochastic parameterizations from physical principles promises accuracy and robustness (Berner et al., 2017). Several authors have addressed the issue of stochastic subgrid parameterization evaluation through averaging or homogenization procedures (e.g. Majda-Timofeyev-Vandeneijnden - MTV - method; see Majda et al., 1999). Correlated additive and multiplicative noises usually play a key role in these approaches (e.g. Gottwald and Melbourne, 2013). Sardeshmukh and Sura (2009), Sardeshmukh et al. (2015), and Pearson and Fox-Kemper (2018) highlight the relevance of skew-symmetric and multiplicative noises in geophysical fluid dynamic applications. Mimetic symmetries built into the fluid dynamics and stochastic subgrid model formulations may preserve some important physical invariants - a desired properties of many UQ models (e.g. Frank and Gottwald, 2013; Sapsis and Majda, 2013; Mitchell and Gottwald, 2012; Majda, 2015) and DA algorithms (Janjić et al., 2014). Other stochastic subgrid parameterizations involve memory terms, justified by the Mori-Zwanzig equations (Chekroun et al., 2011; Lu et al., 2017).

The dynamics under location uncertainty (LU) are a type of systematic stochastic subgrid parameterization introduced by Mikulevicius and Rozovskii (2004b) and Flandoli (2011) for the theory of well-posedness of stochastic partial differential equations and by Mémin (2014) for more applied purposes. More theoretical results about models under location uncertainty are discussed in Resseguier et al. (2017a). This framework has been successfully applied to uncertainty quantification (Resseguier et al., 2017b), geophysics (Resseguier et al., 2017c), attractor exploration (Chapron et al., 2018), and optical flow (Cai et al., 2018).

A similar class of models called stochastic advection by Lie transport (SALT) was first derived in Holm (2015) for the time-homogeneous subgrid parameterization and extended to the inhomogeneous case in Gay-Balmaz and Holm (2018). In the following, we will refer to this class of models as SALT. The approach used is Hamilton's variation principle with the constraint that fluid parcels move according to a Eulerian stochastic subgrid model; see Holm (2015) and Cotter et al. (2017). A local well-posedness result for the 3D SALT Euler equations is proved in Crisan et al. (2019), and a global well-posedness result for the 2D case is proved in Crisan and Lang (2019). Numerical studies for the implementation of the SALT subgrid parameterization are discussed in Cotter et al. (2019a, 2018) for a 2D Euler model and two-layer 2D QG dynamics. A first study on employing SALT dynamics for data assimilation is provided in Cotter et al. (2019b).

In both LU and SALT frameworks, the velocity is decomposed into a random large-scale component, $w$, and a time-uncorrelated component, $\sigma \dot{B} = \sigma dB_t/dt$. The latter is Gaussian (conditionally on some large-scale quantities), correlated in space, with possible heterogeneities and anisotropy.

Hereafter, this unresolved velocity component will further be assumed to be solenoidal. To parameterize those spatial correlations, we apply an infinite-dimensional linear operator, $\sigma$, to a $d$-dimensional space-time white noise, $\dot{B}$. Holm and coauthors rely on a different but equivalent notation. It directly explicates the spectral decomposition of the time-uncorrelated velocity spatial covariance. Specifically, the small-scale velocity reads $\sigma(\mathbf{x})\dot{B} = \sum_p \xi_p(\mathbf{x})dW_p(t)/dt$, where $(W_p)_p$ is a set of independent 1D Brownian motions and $(\xi_p)_p$ is a set of orthogonal functions.

That unresolved velocity is the central object of both SALT and LU approaches. The statistics of this velocity - defined either by $\sigma$ or $(\xi_p)_p$ - constitute the parameterization of the random model. In this paper, we discuss several choices for this parameterization and dedicate significant discussion to the physical principles and assumptions about the flow regime underlying each choice. To better motivate this work, we will highlight the strong links and some formal differences between the LU and SALT approaches. However, a full - theoretical and numerical - comparison of LU and SALT is beyond the scope of this paper.

The paper is structured as follows.

- Section 2.1 focuses on parametric non-data-driven choices. Specifically, we propose a non-stationary and tuning-free improvement of the work of Resseguier et al. (2017b) based on self-similar assumptions and possible heterogeneous modulation. Links with the Smagorinsky subgrid parameterization and non-linear energy fluxes are also discussed.

- Section 2.2 discusses a data-driven stationary but possibly heterogeneous parameterization implemented by Cotter et al. (2019a, 2018). The method relies on decomposition of empirical orthogonal functions (EOF) also known as proper orthogonal decomposition (POD) or principal component analysis (PCA).

- Section 3 begins by presenting the deterministic and stochastic surface quasi-geostrophic (SQG) models, their numerical setup, and simulation parameters. These are followed by detailed discussions of uncertainty quantification skills in Sect. 3.4.

- Sections 4 and 5 concludes this work.

Additionally, in Appendix A we include a mathematical summary of LU and SALT formulations for the interested readers. The aim here is to highlight the similarities and differences between the two approaches. Their differing sets of invariants are listed.

The choice of surface quasi-geostrophic flow is convenient, since the SALT and LU versions of the SQG dynamics coincide; thus the results presented apply to both SALT and LU models. Furthermore, present interest in oceanic submesoscale motions for which SQG dynamics are often a reasonable approximation (e.g. Lapeyre and Klein, 2006b; Cal-
lies et al., 2016) and the fact that SQG offers the possibility of reconstructing subsurface flow properties from surface satellite observations (e.g. LaCasce and Mahadevan, 2006; Isern-Fontanet et al., 2006, 2008; Klein et al., 2009; Wang et al., 2013) have resulted in a recent surge in the use of the SQG equations for geophysics. Accordingly, the discussion can focus on the parameterization choices in this timely application.

## 2 Parameterization of SALT and LU models: the statistics of the unresolved velocity

### 2.1 Parametric and self-similar model for the unresolved velocity

In this section, we propose parametric non-data-driven models for the unresolved velocity statistics. Similar to the Kraichnan model (Kraichnan, 1968, 1994; Gawedzky and Kupiainen, 1995; Majda and Kramer, 1999), Resseguier et al. (2017b) introduce an unresolved velocity which is homogeneous and isotropic. Even though this method leads to good numerical results, some parameters need to be tuned. To overcome this drawback, we propose here a parameter-free improvement based on what we call the absolute diffusivity spectral density (ADSD). New subgrid statistics will be defined at each time step from the resolved velocity kinetic energy (KE) spectrum. Finally, we also propose a heterogeneous modulation of the previous method based on the local energy flux and discuss the link with Smagorinsky-like subgrid parameterizations.

#### 2.1.1 Parameter-free and non-stationary homogeneous model

We first define the statistics of the small-scale velocity through what we call the ADSD. Then, we explicate how to sample the unresolved velocity from this ADSD.

In Resseguier et al. (2017b), the absolute diffusivity (i.e. KE times correlation time; Falkovich et al., 2001; Penland, 2003; Klyatskin, 2005; Vallis, 2006; Keating et al., 2011) of the unresolved velocity is twice the variance tensor trace, $tr(a) = E||σdB_t||^2/dt$, whereas the unresolved kinetic energy is $tr(a)/dt$. Clearly, this kinetic energy has no physical meaning. Indeed, it depends on the simulation time step, and one should have the possibility to choose the time step as close as possible to zero. Thus, the unresolved velocity amplitude is specified through an absolute diffusivity rather than through a KE. In the mathematics literature of homogenization, Kubo-type formulas may be seen as what physicists call absolute diffusivities. More generally, since the variance of a time-continuous white noise is infinite, it is more relevant to deal with absolute diffusivity rather than kinetic energy in order to describe the statistics of the time-uncorrelated velocity. Thus, keeping a spectral approach, we define – for any spatiotemporal field – an absolute diffusivity spectral den-

sity denoted $A(κ)$ at the wavenumber $κ$. We will rely on this ADSD rather than on the KE spectrum, $E(κ)$. Since the absolute diffusivity is the variance multiplied by the correlation time, it is naturally to define the ADSD as follows:

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

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

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

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

$$E(\kappa) = C^2 \kappa^{-s}, \quad (3)$$

and we obtain

$$A(\kappa) = C\kappa^{-r}, \quad (4)$$

where $r = (s + 3)/2$.

We aim at defining the unresolved velocity ADSD from the large-scale velocity. For this purpose, we will assume the self-similar model (Eq. 4) is valid at all spatial scales. At each time step, we compute the ADSD of the large-scale velocity, $A_w$, from Eq. (2). Then, we fit the coefficients $C$ and $r$ of Eq. (4), as illustrated by Fig. 1. Let us denote these coefficients with $C_w$ and $r_w$. Note that they are time-dependent because they depend on $w$ which is. More precisely, we estimate the coefficients $C_w$ and $r_w$ in a wavenumber interval which approximately represents an inertial range of fully resolved scales (i.e. before the spectrum roll-off). In the left panel of Fig. 1, this interval is delimited by two vertical lines.

From there, we can define the unresolved velocity ADSD in such a way that the total velocity – resolved plus unresolved – meets Eq. (4) at small spatial scales with the same coefficients ($C_w$ and $r_w$). Since the two velocity components are not correlated, the total ADSD is the sum of the ADSD of each velocity component. In the inertial range, this reads

$$A(\kappa) = \underbrace{A_w(\kappa)}_{\text{Leading at large scales}} + \underbrace{A_{\sigma\dot{B}}(\kappa)}_{\text{Leading at small scales}} \approx C_w \kappa^{-r_w}. \quad (5)$$

Therefore, the unresolved ADSD is chosen to compensate the resolved ADSD roll-off – introduced by the deterministic subgrid parameterization – at small scales. Specifically, the unresolved ADSD is set to

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

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

$$\kappa_M = \frac{\pi}{\Delta x}, \quad (7)$$
<figure><img src="image_2.png" /><figcaption>Figure 1. Buoyancy field (m s⁻²) (a, b), KE spectrum (m² s⁻²/(rad m⁻¹)) (c, d), and ADSD (m² s⁻¹/(rad m⁻¹)) (e, f) of the resolved velocity, at t = 0 d (a, c) and t = 15 d (b, d, f), in blue; ADSD of the unresolved velocity (without a multiplicative constant), in green; and slope −5/3 (c, d) and corresponding ADSD (e, f) in a black solid line. The two dashed vertical lines define the interval where coefficients C<sub>w</sub> and r<sub>w</sub> are fitted. The right-hand-side dashed vertical line is fixed (see Appendix B for its specification). The left-hand-side dashed line corresponds to the energy-injecting scale. This scale is estimated – at each time step – as corresponding to the maximum of the compensated KE spectrum (i.e. KE spectrum divided by a theoretical spectrum). The dashed oblique line is the resulting fit (it is set to meet A<sub>w</sub> in the left bound of the fitting interval). The unresolved velocity is still restricted to a narrow spectral band, specifically between the dotted dashed line and the right end of the plot. Nevertheless, the form of its spectrum varies in time and is determined by the spectrum of the large-scale velocity. This particular initial condition corresponds to case 2 studied by Constantin et al. (1994, 1999, 2012) at a resolution of 128².</figcaption></figure>

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

$$ \kappa_m = \left( \frac{\ln(0.95)}{\ln(0.1)} \right)^{1/p} \kappa_M, \qquad (8) $$

where p is the order of the Laplacian used as a deterministic subgrid tensor (i.e. $\frac{Db}{Dt} = -\nu(-\Delta)^p b$). The justification of the above formula is left in Appendix B. Compared to the work of Resseguier et al. (2017b), the value of $\kappa_m$ is less critical. Indeed, Eq. (6) implies a weaker unresolved ADSD at larger scales where the resolved ADSD, $A_w$, is stronger. This softens the threshold effect introduced by the band-pass filter $f_{BP}$.

In practice, we set an upper bound for the estimation of $r_w$. Without this upper bound, a concentration of energy at relatively large wavenumbers – scales smaller than $\kappa_m$ – in the resolved fields can become unstable. Indeed, this localized energy concentration would decrease the $r_w$ estimation, and hence increase the unresolved ADSD $A_\sigma \dot{B}$ through Eq. (6) at large wavenumbers – larger than $\kappa_m$. This implies a larger

noise intake, which can induce a larger concentration of energy at relatively large wavenumbers in the resolved fields, resulting in a positive feedback loop. To prevent these unphysical instabilities, the slope $r_w$ is bounded.

We have proposed a way to compute the unresolved ADSD $A_\sigma \dot{B}$ from the large-scale velocity statistics. From that ADSD, it is possible to sample the unresolved velocity, as explained in the following.

A 2D homogeneous divergence-free small-scale velocity can be constructed by filtering a 1D white noise $\dot{B}$ with a 2D divergence-free filter $\tilde{\sigma} = \nabla^\perp \tilde{\psi}_\sigma$ (Resseguier et al., 2017b).

$$ \sigma \dot{B} = \tilde{\sigma} \star \dot{B} = \nabla^\perp \tilde{\psi}_\sigma \star \dot{B}, \qquad (9) $$

where $\star$ denotes a spatial convolution. This velocity field can be easily sampled in Fourier space.

$$ \widehat{\sigma \dot{B}}(k) = ik^\perp \widehat{\psi}_\sigma(k) \widehat{B}(k) = \frac{1}{\sqrt{\Delta t}} ik^\perp \widehat{\psi}_\sigma (\|k\|) \frac{\widehat{dB_t}}{\sqrt{\Delta t}}(k), \qquad (10) $$

where $\widehat{dB_t}$ is the spatial Fourier transform of $dB_t$ and $\frac{dB_t}{\sqrt{\Delta t}}$ is a discrete scalar white-noise process of unit variance in
space and time. Thus, the small-scale velocity is defined by the Fourier transform of the streamfunction kernel, $\hat{\psi}_{\sigma}$.

In order to link the unresolved ADSD to the kernel $\hat{\sigma}$ which defines the unresolved velocity, we note that

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

where $\mu(\Omega)$ is the surface of the spatial domain $\Omega$ and $\theta_{\kappa}$ is the angle of the wave vector $\mathbf{k}$. From Eqs. (6)–(11), we can finally express the unresolved velocity as follows:

$$\hat{\sigma} \hat{B}(t, \mathbf{k}) = \frac{1}{\sqrt{\Delta t}} i \mathbf{k}^{\perp} \hat{\psi}_{\sigma}(t, \|\mathbf{k}\|) \frac{\widehat{dB_t}}{\sqrt{\Delta t}}(\mathbf{k}), \qquad (12)$$

$$= \frac{1}{\sqrt{\Delta t}} i \mathbf{k}^{\perp} \sqrt{\frac{\max(0, C_w \|\mathbf{k}\|^r_w - A_w (\|\mathbf{k}\|))}{2\pi \|\mathbf{k}\|^3}}$$

$$f_{BP}(\|\mathbf{k}\|) = \frac{\widehat{dB_t}}{\sqrt{\Delta t}}(\mathbf{k}). \qquad (13)$$

Again the simulated unresolved velocity ADSD is physically relevant, while the KE spectrum is not. Indeed, the simulated unresolved velocity ADSD is expected to match the true (time-correlated) unresolved velocity ADSD, whereas the KE spectra of the simulated and true unresolved velocities differ. Indeed, the true unresolved velocity correlation time spectral distribution $\tau(\kappa)$ is not restricted to the time step $dt$.

The tuning-free Eq. (12) is attractive because it is not a priori restricted to stationary flows or to a particular type of 2D incompressible dynamics. In order to appreciate that, we can numerically compare it with the previous method of Resseguier et al. (2017b). For stationary, fully developed turbulence and after including a tuning stage to optimize the match, both parameterizations give approximately the same results. Therefore, in order to perceive improvements between the two simulations, we must either work with non-stationary flows or with erroneous tuning. This latter scenario is of course the meaningful one for the application of such models to the real world, where turbulence is non-stationary and heterogeneous, and so region-by-region and season-by-season tuning is impossible to do with quality checks in place.

We believe that non-stationarity may yield a fair, yet idealized, comparison. So, we here compare the two parameterizations in a non-stationary case. Figure 2 below shows simulated buoyancy fields initialized with “case 2” of Constantin et al. (1994) and corresponding errors. After 2 d of advection, there is no turbulence yet, and a 128 × 128 resolution is sufficient to correctly resolve every scale. Therefore,

no stochastic subgrid parameterization is needed. The self-similar method automatically adapts to the situation, whereas the method of Resseguier et al. (2017b) introduces spurious buoyancy isoline roughness by randomly folding the isolines. Accordingly, the method of Resseguier et al. (2017b) introduces more errors.

## 2.1.2 Heterogeneity

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

In order to obtain a heterogeneous model of the unresolved velocity, we need a heterogeneous version of the ADSD (Eq. 4). Since the wavenumber, $\kappa$, cannot depend on the position, $\mathbf{x}$, the constant, $C_w$, and/or the spectrum slope, $r$, should do. A spatially varying spectrum slope is probably difficult to estimate. Hence, we restrict ourselves to a spatially varying constant, $C_w$, and a spatially homogeneous spectrum slope. The constant may also varies with the time and the wavenumber. According to the Kolmogorov theory (e.g. Frisch, 1995) and Eq. (3):

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

where $\epsilon_F$ is the energy flux through the spatial scales and $p = 1/3$ for a SQG flow (cst. is an abbreviation for constant). More specifically, the energy flux describes the energy moving from scales larger than $1/\kappa$ toward smaller scales and can be computed as follows:

$$\epsilon_F(\mathbf{x}, t, \kappa) \triangleq \overline{q_{\kappa}^{<}((\mathbf{w} \cdot \nabla) q)_{\kappa}^{<}}, \qquad (15)$$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Buoyancy field (m s⁻²) (**a**, **b**, **c**) and corresponding normalized error (dimensionless) (**d**, **e**) after 1 d of advection for the 1024 × 1024 reference deterministic simulation (**a**), the 128 × 128 LU simulation with the self-similar parameterization (**b**, **d**), and the 128 × 128 LU simulation with the method of Resseguier et al. (2017b) (**c**, **e**).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Presented from left to right are the buoyancy field (m s⁻²), KE spectrum (m² s⁻²/(rad m⁻¹)), ADSD (m² s⁻¹/(rad m⁻¹)), and variance tensor at *t* = 15 d for (from top to bottom) a 1024²-resolution deterministic dynamics and a 128²-resolution deterministic dynamics. This particular initial condition corresponds to case 2 studied by Constantin et al. (1994, 1999, 2012) at a resolution of 128². The low-resolution deterministic dynamics shows too many filaments and not enough eddies.</figcaption></figure>

where *q* is the transported (up to possible source terms) quantity, $g_K^<$ is the low-pass filtered version of *g* (setting to zero the Fourier modes of *g* which have frequencies larger than $κ$), and ⋙ stands for the spatial average (Frisch, 1995). For a SQG flow, *q* corresponds to the buoyancy normalized by the stratification: $q = b/N$. The energy flux is essentially a third-

order moment. It is very important because it describes the cascade of the flow by non-linear energy transfers (Frisch, 1995).

If the energy flux through scale is understood locally in space (as indeed Smagorinsky, 1963, also assumes), Eq. (14) provides a natural parameterization of the unresolved ve-
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Presented from left to right are the buoyancy field (m s⁻²), KE spectrum (m² s⁻²/(rad m⁻¹)), ADSD (m² s⁻¹/(rad m⁻¹)), and variance tensor at t = 15 d for 128²-resolution stochastic dynamics with (from top to bottom) homogeneous unresolved velocity, unresolved velocity modulated by ||∇b||¹/⁴, and unresolved velocity modulated by ϵ_F¹/⁶. This particular initial condition corresponds to case 2 studied by Constantin et al. (1994, 1999, 2012) at a resolution of 128². The low-resolution stochastic simulations do not show too many filaments and show enough eddies. Among the stochastic simulations, the energy-flux modulation better preserves the sharp straight fronts (e.g. front from (0 m, 2 × 10⁵ m) to (2.5 × 10⁵ m, 5 × 10⁵ m)).</figcaption></figure>

locity heterogeneities. We simply modulate the unresolved ADSD (Eq. 6) by the heterogeneous ratio ϵ_F^p/ϵ_F^p, averaged over the resolved inertial range wavenumbers.

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

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

where $\mathcal{P} = I_d - \Delta^{-1}\nabla\nabla^T$ is the projector onto the space of free-divergence functions. This parameterization is physically meaningful, since, locally in space, a stronger direct cascade at large scales (larger $\epsilon_F$ and thus larger $C_w$) suggests that the unresolved velocity (large) should maintain this cascade by folding smaller-scale tracer structures. Furthermore, considering that the energy flux is a third-order structure makes this parameterization relevant to differentiate between strait fronts and curved structures (e.g. eddies). Indeed, at least three points are needed to define a curvature and differentiate between these structures. Figure 4 confirms that this modulation enables a more accurate spatial distribution of the stochastic folding.

In order to keep a divergence-free velocity and the ensuing properties (e.g. energy conservation), the modulated velocity is projected onto the space of free-divergence functions, using the operator $\mathcal{P}$. Because of that we do not consider the
advection correction $w^* - w = -\frac{1}{2}(\nabla \cdot a)^T$ of the LU formalism. Indeed, here the variance tensor has the simple form $a = \frac{1}{d} \text{tr}(a) I_d$. As such, the advection correction is a gradient field and is hence removed by the projection onto the space of free-divergence functions.

Although relevant, the comparison of Figs. 3 and 4 – about front dynamics – remains qualitative. For a quantitative demonstration, adapted metrics should be used. Indeed, simple isotropic, homogeneous second-order statistics (e.g. KE spectra) cannot distinguish between eddies and filaments. Bi-spectra may overcome this drawback, since they express three-point statistics. This quantitative analysis would necessitate studies of the metrics themselves. Therefore, these analyses will be addressed in future work.

Many other closures rely on the Kolmogorov (Kolmogorov, 1941) model (Eqs. 3–14): in particular, the famous Smagorinsky model (Smagorinsky, 1963) and its variants (e.g. Fox-Kemper and Menemenlis, 2008; Bachman et al., 2017) provide a path to developing deterministic and dissipative scale-aware subgrid models. Typically, these models result in a Laplacian dissipation which involves a heterogeneous eddy diffusivity or viscosity coefficient, $v_{\text{Sm}}$. Aside from their heuristic theoretical justification, these Smagorinsky-type subgrid terms are formally equivalent to the turbulent dissipation of our stochastic model: $\nabla \cdot (\frac{1}{2} a \nabla q)$ (Resseguier et al., 2017a), where $a = 2v_{\text{Sm}} I_d$ and $q$ is a transported quantity. This similarity suggests that a Smagorinsky-type model could provide a good estimate for our variance tensor and thus for the heterogeneity of the unresolved velocity.

The goal of the Smagorinsky model is to optimize the KE spectrum by targeting a specific turbulent diffusive scale adapted to the simulation resolution. A turbulent dissipation coefficient expression can be derived from the Kolmogorov model (Eqs. 3 and 14) and the closure,

$$ \epsilon_F = \epsilon_D, \qquad (18) $$

where $\epsilon_D$ is the dissipation, a second-order moment related to the molecular or turbulent diffusion. To develop a Smagorinsky-type model, the resolved flux of the energy-like conserved invariant is equated to the dissipation of the energy invariant such that

$$ \epsilon_D = v_{\text{Sm}} \|\nabla q\|^2. \qquad (19) $$

From there, one can obtain an eddy diffusivity or viscosity coefficient proportional to $\|\nabla q\|^h$. For an SQG flow, the exponent is $h = 1/2$, where $q = b$ is the buoyancy (Bachman et al., 2017).

In the Kolmogorov theory of homogeneous and stationary turbulence, the energy flux is a constant of the flow, and the closure (Eq. 19) is an exact result of a simple energy budget over an ensemble mean. Indeed, there is no accumulation of energy at any scales in a stationary regime. This closure is

very useful since the dissipation (Eq. 19) is generally much simpler to compute than the energy flux. Nevertheless, in every flow realization, the energy flux and the diffusion vary with space and time, and they do not match each other locally. For instance, a strong straight front of an SQG flow involves a large dissipation but no energy cascade because the velocity is aligned with front (see Eq. 15). Moreover, in any bounded, limited resolution situation, the inertial cascade range is limited so the energy flux through scale varies with the wavenumber, especially outside the inertial range.

The discrepancy between energy flux and dissipation is not so much of an issue for the Smagorinsky model because its aim is the optimization of a second-order statistics at small scales. Unfortunately, this closure cannot be used to simplify the modulation computation in our parametric random model. Indeed, the stochastic dynamics relies on processes – such as folding – associated with higher-order statistics. Figure 4 (middle right) illustrates this statement. The following unresolved velocity parameterization is used there:

$$ \begin{gathered} \widetilde{\sigma \dot{B}}(t, k) = \frac{1}{\sqrt{\Delta t}} ik^{\perp} \sqrt{\frac{\max(0, C_w \|k\|^{-r_w} - A_w (\|k\|))}{2\pi \|k\|^3}} \\ f_{\text{BP}}(\|k\|) = \frac{\widehat{dB_t}}{\sqrt{\Delta t}}(k), \end{gathered} \qquad (20) $$

$$ \sigma \dot{B}(t, x) = \mathcal{P} \left\{ \frac{\|\nabla b\|^{1/4}(t, x)}{\sqrt{\|\nabla b\|^{1/2}(t)}} \underbrace{\widetilde{\sigma \dot{B}}(t, x)}_{\text{Heterogeneous modulation}} \right\}. \qquad (21) $$

Along sharp straight fronts, the dissipation will be larger. Accordingly, the associated modulation $\|\nabla b\|^{1/4}$ enhances the stochastic folding where one would need it to weaken.

## 2.2 Data-driven model for the unresolved velocity

In this section, we detail a procedure proposed by Cotter et al. (2019a) for the estimation of the (weighted) EOFs, $(\xi_k)_k$, involved in the unresolved velocity definition,

$$ \sigma \dot{B} = \sum_k \xi_k dW_k/dt. \qquad (22) $$

If we denote $\lambda_k$ as the $L^2$ norm of those EOFs, their normalized versions $(\xi_k = \xi_k/\lambda_k)_k$ are the eigenvectors of the self-adjoint operator,

$$ f \mapsto \int_{\Omega} \left[ \sigma(\mathbf{x}) \sigma^T(\mathbf{y}) / dt \right] f(\mathbf{y}) dy, \qquad (23) $$

defined by the small-scale velocity covariance $[\sigma(\mathbf{x})\sigma^T(\mathbf{y})/dt]$. $(\lambda_k/dt)_k$ is the set of eigenvalues of this operator.

The data-driven methods of Cotter et al. (2019a) relies on Lagrangian paths defined at two “resolutions”. The first paragraph of this section defines these two types of Lagrangian
paths. Then, we propose several ways to identify the time increments of the infinite-dimensional Brownian motion, $\sigma B_t$, from these Lagrangian paths. Preprocessing of the increments are needed in order to meet some structural assumptions. After this, we relate the increment covariance to the EOFs.

### 2.2.1 Preliminary definitions

We introduce two types of velocity field:

- a high-resolution velocity, $v$, on a fine mesh-grid (512²),

- a low-resolution velocity, $\bar{v}$, on a coarse mesh-grid (64²). This velocity field is a spatially low-pass-filtered version of $f$.

Then, two types of Lagrangian path are defined:

- a “high-resolution flow”, $X_{ij}(t_0, t)$, defined by the high-resolution velocity, $u$:

$$ \frac{d\mathbf{X}_{ij}(t)}{dt} = \mathbf{u}(t, \mathbf{X}_{ij}(t_0, t)) \quad \text{and} \quad \mathbf{X}_{ij}(t_0, t_0) = \mathbf{x}_{ij}, \text{ and} \tag{24} $$

with $\mathbf{x}_{ij} = (i \Delta x, j \Delta y)$.

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

$$ \frac{d\bar{\mathbf{X}}_{ij}(t_0, t)}{dt} = \bar{\mathbf{u}}(t, \bar{\mathbf{X}}_{ij}(t_0, t)) \quad \text{and} \quad \bar{\mathbf{X}}_{ij}(t_0, t_0) = \mathbf{x}_{ij}. \tag{25} $$

### 2.2.2 Candidate for the increment realization

In order to estimate the EOFs, $(\xi_i)_i$, involved in Eq. (A5), we assume that we can observe increments

$$ \sigma \Delta \mathbf{B}_{m\Delta T} = \sigma (\mathbf{B}_{(m+1)\Delta T} - \mathbf{B}_{m\Delta T}). \tag{26} $$

We will interpret the following residual flow increments as a realization of the above:

$$ \begin{aligned} \Delta \widetilde{\mathbf{X}}_{ij}^m &\triangleq \mathbf{X}_{ij}(m\Delta T, (m+1)\Delta T) \\ &\quad - \bar{\mathbf{X}}_{ij}(m\Delta T, (m+1)\Delta T). \end{aligned} \tag{27} $$

### 2.2.3 Preprocessing

The increments are supposed to be centred and divergence free (as $\nabla \cdot \xi_i = 0, \forall i$). Therefore, after computing the residual flow increments, $(\Delta \widetilde{\mathbf{X}}_{ij}^m)_m$, they are centred,

$$ \Delta \mathbf{X}_{ij}'' \triangleq \Delta \widetilde{\mathbf{X}}_{ij}^m - \hat{\mathbb{E}}\{\Delta \widetilde{\mathbf{X}}_{ij}^m\}, \tag{28} $$

with the estimator

$$ \hat{\mathbb{E}}\{\Delta \widetilde{\mathbf{X}}_{ij}^m\} \triangleq \frac{1}{N} \sum_{m=0}^{N-1} \Delta \widetilde{\mathbf{X}}_{ij}^m, \tag{29} $$

and projected onto the space of divergence-free functions,

$$ \Delta \mathbf{X}_{ij}^m \triangleq \mathcal{P}\{\Delta \mathbf{X}_{ij}''\}, \text{ with } \mathcal{P} = \mathbf{I}_d - \Delta^{-1}\nabla\nabla^T. \tag{30} $$

### 2.2.4 Covariance, quadratic covariation, and EOF

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

$$ \begin{aligned} \sum_k \xi_k(\mathbf{x}_{ij}) \xi_k^T(\mathbf{x}_{pq}) \\ &= \sigma(\mathbf{x}_{ij}) \sigma^T(\mathbf{x}_{pq}) \end{aligned} \tag{31} $$

$$ \approx \frac{1}{(N-1)\Delta T} \sum_{m=0}^{N-1} \Delta \mathbf{X}_{ij}^m (\Delta \mathbf{X}_{pq}^m)^T \tag{32} $$

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

$$ \int_{\Omega} \xi_i(\mathbf{x})^T \xi_j(\mathbf{x}) d\mathbf{x} = 0 \quad \text{if } i \neq j. \tag{33} $$

Finally, after estimating the $(\xi_k)_k$ offline, the ensemble forecast can generated online with Eq. (22).

## 3 Numerical simulations and uncertainty quantification

### 3.1 Surface quasi-geostrophic model

#### 3.1.1 Deterministic SQG

This is a simplified model to describe the ocean surface dynamics at mesoscales (i.e. horizontal length scale of the order of 100 km) (Blumen, 1978; Held et al., 1995; Lapeyre and Klein, 2006a; Constantin et al., 1994, 1999, 2012; Lapeyre, 2017). The buoyancy, $b$, is transported at the ocean surface by a horizontal velocity field,

$$ \frac{Db}{Dt} = S, \tag{34} $$

where $D/Dt$ stands for the horizontal (deterministic or stochastic) material derivative and $S$ represents possible sources and sinks. As the potential vorticity is assumed to be zero in the fluid interior, the Fourier transform of the velocity streamfunction, $\hat{\psi}$, is related to the Fourier transform of the buoyancy, $\hat{b}$, by the following SQG relationship:

$$ \hat{\psi} = \frac{1}{N \|\mathbf{k}\|} \hat{b}, \tag{35} $$

where $N$ is the stratification and $\mathbf{k}$ is the wave vector.

#### 3.1.2 Stochastic SQG

The LU and SALT versions of the SQG dynamics are formally similar to the deterministic model. However, the buoyancy transport (Eq. 34) has to be understood in the stochastic
sense (Eq. A4). Furthermore, the SQG relationship (Eq. 35) must be interpreted with

$$ w^* = \nabla^{\perp} \psi = \nabla^{\perp} (-\Delta)^{-1/2} (b/N), \quad (36) $$

in the SALT context, and $w = \nabla^{\perp} \psi = \nabla^{\perp} (-\Delta)^{-1/2} (b/N)$, in the LU one. Besides, in the LU SQG, the advecting drift $w^*$ has to be divergence-free. We enforce this constraint by projecting the drift correction ($w^* - w$) onto the space of free-divergence functions. So, we have

$$ w^* = P(w^* - w) + w = -\frac{1}{2} P(\nabla \cdot a)^T + \nabla^{\perp} (-\Delta)^{-1/2} (b/N), \quad (37) $$

where $P = I_d - \Delta^{-1}\nabla\nabla^T$. Indeed, the SQG model is derived from a QG model with a transport of the buoyancy at the surface. Then, the potential vorticity (PV) is assumed to be zero inside the fluid. In the SALT framework, the PV reads $0 = \text{PV}_{\text{slt}} = \nabla^{\perp} \cdot w^* + f + (f_0/N)^2 \partial_z^2 \psi$, $w^* = \nabla^{\perp} \psi$ and $b = f_0 \partial_z \psi$. For the LU dynamics, $0 = \text{PV}_{\text{lu}} = \nabla^{\perp} \cdot w + f + (f_0/N)^2 \partial_z^2 \psi$, $w = \nabla^{\perp} \psi$ and $b = f_0 \partial_z \psi$.

Nevertheless, the slight difference between the SQG SALT and the SQG LU models are not considered in this section since we neglect the advection correction of the LU framework ($w^* - w$). So, we simulate the stochastic transport (Eq. 34) coupled with the SQG relation (Eq. 36). The unresolved velocity statistics encoded in $\sigma$ are specified either by the self-similar method of Sect. 2.1.1 or by the data-driven method of Sect. 2.2.

### 3.1.3 Flow simulation

We perform a high-resolution simulation (512² spatial grid) of the deterministic SQG model with the following initial condition:

$$ b(x, t=0) = 0.2 B_0 \cos \left( \frac{2\pi}{L} (x+y) \right) + F(x-x_{00}^0) + F(x-x_{10}^0) - F(x-x_{01}^0) - F(x-x_{11}^0), \quad (38) $$

$$ \text{with } F(x) = B_0 \exp \left( -\frac{1}{2} \left( \frac{15}{L} \right)^2 \left( x^2 + \left( \frac{y}{2} \right)^2 \right) \right), \quad (39) $$

$$ \text{and } x_{ij}^0 = \frac{L}{4} \begin{pmatrix} 1 \\ 1 \end{pmatrix} + \frac{L}{2} \begin{pmatrix} i \\ j \end{pmatrix}. \quad (40) $$

The source term, $S$, involves a hyperviscosity, a linear drag, and an additive stationary forcing such that

$$ S \triangleq -\nu_{\text{HV}} \Delta^4 b - \frac{1}{\tau_F} b + A_F \sin(k_F^x x) \sin(k_F^y y). \quad (41) $$

The parameters of the simulations are summed up in Table 1. The influence of the initial condition remains for about a

<table><caption>**Table 1.** Parameters of the simulation.</caption><thead><tr><th>Parameters</th><th>Value</th></tr></thead><tbody><tr><td><i>L</i> (domain length)</td><td>10<sup>6</sup> m</td></tr><tr><td><i>B</i><sub>0</sub>/<i>N</i></td><td>3.24 m s<sup>-1</sup></td></tr><tr><td>√<span style="text-decoration: overline;">b<sup>2</sup></span>(<i>t</i> = 0)/<i>N</i></td><td>1.16 m s<sup>-1</sup></td></tr><tr><td>√<span style="text-decoration: overline;">b<sup>2</sup></span>(<i>t</i> = 100 d)/<i>N</i></td><td>2.28 m s<sup>-1</sup></td></tr><tr><td><i>A</i><sub>F</sub>τ<sub>F</sub>/<i>N</i></td><td>32.4 m s<sup>-1</sup></td></tr><tr><td>1/τ<sub>F</sub></td><td>2.16 × 10<sup>-7</sup> s<sup>-1</sup> ≈ 1/(54 d)</td></tr><tr><td>(<i>k</i><sub>F</sub><sup>x</sup>, <i>k</i><sub>F</sub><sup>y</sup>)</td><td>(3(2π/L), 2(2π/L))</td></tr><tr><td>ν<sub>HV</sub>/dx<sup>8</sup></td><td>3.39 × 10<sup>-8</sup> s<sup>-1</sup></td></tr></tbody></table>

month. Then, the turbulence is maintained by the forcing as illustrated by Fig. 5.

From $t = 50$ d to $t = 100$ d, the EOFs of the data-driven method are learned. From $t = 100$ d to $t = 130$ d, the deterministic high-resolution simulation is used as a reference. In this time interval, several stochastic low-resolution (64²) simulations are performed and compared. These simulations are initialized at $t = 100$ d with the reference high-resolution simulation projected at low resolution (i.e. keeping only the Fourier modes associated with a coarse-resolution grid).

## 3.2 Learning and analysis of the EOF

Here, we describe the convergence of the EOF estimation. The accuracy of this estimation is in particular a function of the number of snapshots used and of the Lagrangian advection time, $\Delta T$. We first describe the convergence with the number of snapshots and then – as a remark – the convergence with the Lagrangian advection time, $\Delta T$.

The EOFs are learned between $t = 50$ d and $t = 100$ d from 12 465 spatial fields of residual flow increments. Even though a large number of snapshots are used, the correlation time of the residual flow increments is about 1 d. This is not negligible compared to the estimation time window: 50 d. Accordingly, the EOFs are not fully converged. But, we expect this convergence to be sufficient for the present work. Moreover, for real applications at this spatial scale, we expect the unresolved velocity statistics to be non-stationary on temporal scales larger 50 d. Thus, learning a unresolved velocity stationary statistical representation on a larger time window might be difficult in practice. Therefore, even if our EOF estimation is not converged, it represents a realistic test case.

**Remark 1.** The integration time of the flows, $\Delta T$, is a critical parameter for the definition of the EOF. Indeed, for small advection time, $\Delta T$, the length of an increment of a Lagrangian flow path is proportional to the velocity and to the advection time. This is the so-called ballistic regime (Falkovich et al., 2001; Vallis, 2006). In particular, residual
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Buoyancy (m s⁻²) (left), KE spectrum (m² s⁻²/(rad m⁻¹)) (middle), and ADSDs (m² s⁻¹/(rad m⁻¹)) (right) at *t* = 0, 30, 50 and 70 d of advection, for the deterministic SQG model at a resolution of 512².</figcaption></figure>

flow increments squared norm would be proportional to $\Delta T^2$; the variance tensor estimator would be proportional to $\Delta T$; and the estimated EOFs would be proportional to $\sqrt{\Delta T}$. With a larger advection time, $\Delta T$, the Lagrangian velocity decorrelates – along the flow path – from the initial Lagrangian velocity. When the advection time becomes larger than the correlation time of the Lagrangian velocity, the flow path begins to act as a Brownian motion (Falkovich et al., 2001; Penland, 2003; Klyatskin, 2005; Vallis, 2006; Keating et al., 2011). The length of a displacement scales

as $\sqrt{\Delta T}$; i.e. the Lagrangian velocity acts as a white noise in time. This is the so-called diffusive regime. Figure 6 illustrates this convergence with the average tensor $a_0 = \frac{1}{d} \text{tr}(a)$ with $\Delta T$. In this paper, the Lagrangian advection time $\Delta T$ is computed from a CFL (Courant–Friedrichs–Lewy condition) at the coarse resolution of 64² – about 300 s. Although it corresponds to the ballistic regime (i.e. the flow increments are correlated), this choice is coherent with the work of Cotter et al. (2019a) and gives very good UQ results. Moreover, the residual flow increments – and hence the EOFs – are
<figure><img src="image_7.png" /></figure>

**Figure 6.** Estimation of the spatial average of the trace of the variance tensor divided by $d = 2$, $\frac{1}{Md} \sum_{q=1}^{M} \text{tr}(\mathbf{a}(\mathbf{x}_q)) = \frac{1}{Md} \sum_{q,i} \|\xi_i(\mathbf{x}_q)\|^2$, as a function of the Lagrangian advection time $\Delta T$, for method 1 (described in Sect. 2.2) (blue plot) and for method 2 (the flow increments defined by the integration of the small-scale velocity along the high-resolution flow: $\Delta \bar{X}_{ij}^m = \int_{k\Delta T}^{(k+1)\Delta T} (\mathbf{u} - \bar{\mathbf{u}})(t, \mathbf{X}_{ij}(m\Delta T, t))dt$) (black plot). The red plot relies on another method not described here. The bottom plot is an enhancement of the top plot. For low $\Delta T$ values, the Lagrangian velocities remain highly correlated during the interval $[0, \Delta T]$, i.e. during the Lagrangian advection. It is the ballistic regime where the flow increments scale linearly in time: $\|\Delta \bar{X}_{ij}^m\| \propto \Delta T$. In that regime, methods 1 and 2 are very similar, since the high-resolution, $\mathbf{X}_{ij}(m\Delta T, t)$, and low-resolution, $\bar{\mathbf{X}}_{ij}(m\Delta T, t)$, flows remain close to each other. Above $\Delta T = 1$ d, the final and initial small-scale Lagrangian velocities, $(\mathbf{u} - \bar{\mathbf{u}})$, are decorrelated. It is the diffusive regime where the flow increments scale as the square root of time: $\|\Delta \bar{\mathbf{X}}_{ij}^m\| \propto \sqrt{\Delta T}$. This is the relevant regime for the estimation of the EOF, since we meet the fundamental assumption of our model: a Brownian behaviour for the small-scale flow. For method 1, the diffusive regime is not visible, at least in this advection time window, and the EOF estimation is theoretically not possible. In contrast, method 2 provides a converged estimator of the variance tensor for an advection time $\Delta T > 1$ d.

spatially aliased, since a large part of those increments lives at small spatial scales but are spatially sampled on a coarse spatial grid, $\mathbf{x}_{ij}$. Figure 7 confirms this idea. If all the $2 \times 64^2 - 1 = 8191$ EOFs are considered, the ADSD of the unresolved velocity reveals a strong spatial aliasing.

Both the data-driven and the non-data driven methods exhibit large unresolved velocities at the smallest scales of the

coarse resolution grid (see Fig. 7). Nevertheless, the two ADSDs are distinct. In particular, the data-driven method involves some large spatial scales. One could think that these large-scale components would disappear if the flow $\mathbf{X}$ and $\bar{\mathbf{X}}$ are integrated in a Eulerian way rather than in a Lagrangian way. However, with the advection time, $\Delta T$, being very small, few differences are expected between a Eulerian and a Lagrangian advection. New numeric experiments confirm this idea (not shown). A complete study of these effects is beyond the scope of this paper.

### 3.3 One realization

We now simulate the LU SQG dynamics (equivalent to SALT SQG) (Eqs. 34–36), at low resolution ($64^2$), with two possible parameterizations for the unresolved velocity: either the data-driven model (see Sect. 2.2) or the parametric and self-similar model (see Sect. 2.1). For the data-driven model, we keep 200 EOFs. This choice will be explained in the next section. For all simulations, there is a unique reference initial condition at $t = 100$ d. This initial condition is the low-resolution ($64^2$) projection of the reference deterministic high-resolution ($512^2$) simulation (i.e. the Fourier modes associated with large wave vectors are set to zero, and the obtained spatial field is subsampled at the low resolution). Spatial fields, KE spectra and ADSD are plotted in Figs. 8 and 10. For comparison purposes, we also show the low-resolution deterministic SQG simulations and the low-resolution ($64^2$) projection of the reference deterministic high-resolution simulation. As already pointed out by Resseguier et al. (2017b) for free-decaying turbulence, the realizations of the LU dynamics are no worse than a low-resolution deterministic simulation. For the short-term simulations, our stochastic subgrid parameterizations have often weak improvements on the low-resolution simulations, even though, sometimes, the stochastic subgrid parameterization can improve the simulation. Indeed, Resseguier et al. (2017b) show that the LU dynamics at a resolution of $128 \times 128$ can trigger filament instabilities by random destabilization and hence obtain a more realistic proportion of eddies and filaments. This is confirmed by Figs. 3 and 4, also at a resolution of $128 \times 128$. In Fig. 8, the resolution is coarser ($64 \times 64$). Therefore, the stabilizing deterministic subgrid tensor (hyper viscosity) is stronger. This may explain an inhibition of filament instabilities here and hence less difference between deterministic and stochastic coarse simulations. Nevertheless, our main goal is not improving a single simulation. Our main goal is improving the uncertainty quantification – as developed below – without deteriorating single simulations.

At $t = 110$ d, all the spatial fields of Figs. 8 and 9 are still similar, whereas at $t = 120$ d, the spatial fields strongly differ. Thus, the predictability timescale – i.e. the time over which initial conditions are forgotten – is between 10 and 20 d. The larger features evolve more slowly, leading to a longer predictability timescale. This effect can be seen in the similar lo-
cation across simulations of larger vortices in Fig. 10, while the filaments between the vortices have lost all coherence. In Fig. 11, the larger, more-coherent features persist in the ensemble mean of the coarse-resolution simulations, while the smaller scales cancel. It is the persistent features which are the basis for robust forecasts as will become clear in the next section, and it is the improvement of these robust features that is the goal of the parameterizations, not the improvement of individual filaments in individual ensemble members.

## 3.4 Uncertainty quantification

We now forecast two ensembles of 200 realizations following the stochastic SQG dynamics (Eqs. 34–36), with the same unique reference initial condition at $t = 100$ d. Again, the ensembles members evolve at low resolution (64²). The first ensemble is generated with the data-driven model for the unresolved velocity (see Sect. 2.2), while the second ensemble is generated with the parametric and self-similar model for the unresolved velocity (see Sect. 2.1). For the data-driven model, we keep only 200 EOFs, since the number of EOFs cannot be larger than the ensemble size without increasing the complexity of the algorithm.

With such ensemble forecasts, we aim at representing the variety of possible behaviours of the fluid dynamic system. In particular, the spreading (i.e. the variance increase) of an ensemble is expected to make an ensemble be closer to the reference. In such a case, the standard deviation – at each point and at each time – is expected to be of the order of the bias. Figure 12 shows that both the data-driven and the self-similar parameterizations achieve this goal everywhere and for every time. The pointwise biases of the data-driven method, of the self-similar method, and of a deterministic simulation (at the same low resolution) are very similar (not shown). Therefore, we only plot the pointwise bias of the data-driven method. In Fig. 12, the represented biases and error estimations (mean ±1.96 times the pointwise standard deviation<fnref n="1" />) are normalized by the squared energy of the reference solution. Note that those relative pointwise biases increase very quickly with time.

After $t = 105$ d and especially after $t = 115$ d, a bifurcation plays a large role in the simulation error and its estimation. This bifurcation is due to the chaotic trajectory of an eddy – centred on (525, 300 km) at $t = 105$ d. Due to the incorrect trajectory in the ensemble mean, a large yellow spot develops in the bias images. Both ensembles capture well this variability by creating similar spots. According to the doubling of those spots, there are probably only two likely trajectories for this eddy, at least until $t = 118$ d.

Similar UQ results have been obtained on a free-turbulence flow (Resseguier et al., 2017b). This close anal-

ysis has also compared the UQ potential of LU/SALT algorithms against that of a deterministic dynamics with random initial conditions. The latter has shown an underestimation of errors by one order of magnitude.

Figure 12 offers a visual validation of our methods’ UQ potential in the whole spatial domain. Nevertheless, the pointwise laws of the ensembles are not Gaussian, since we consider a non-linear evolution law with multiplicative noise. Therefore, the pointwise confidence interval is not in general a symmetric interval centred on the pointwise mean with width $2 \times 1.96$ times the pointwise standard deviation. Such an interval is only a first approximation. In order to be more accurate in our analysis, we now focus on few spatial points. There, we compute – at each time – the true ensemble-based confidence intervals at 95% and at 50% and compare them with the reference value. Figure 13 displays the results at two grid points. Again, the results are impressively good. Similar UQ proxies have been presented by Cotter et al. (2019a) with a SALT 2D Euler dynamics with Dirichlet boundary conditions. This work also highlights the excellent UQ skills of our frameworks.

To conclude this quantitative UQ analysis, we study the effect of the number of realizations needed. How many is a main issue in operational weather forecast centres, since each realization is costly. Accordingly, we start again the UQ analysis but with 20 ensemble members only. Figure 14 shows that the 97.5% and the 2.5% quantiles get closer, compared to the 200-ensemble-member case. It means that the probability density tail estimations – i.e. the representations of extremes – are slower to converge than the more likely values. Nevertheless, the ensembles still capture well the reference dynamics of the centre of the distribution.

## 4 Conclusion

This paper develops the SALT and LU models which coexist in a single family of stochastic schemes. In addition to their general theoretical properties, we have discussed and compared possible parameterizations for this new scheme family.

Appendix A highlights the strong theoretical similarities between SALT and LU stochastic subgrid tensors at the heart of their parameterizations. These frameworks assume that tracers are transported by the sum of a resolved large-scale velocity and an unresolved time-uncorrelated velocity. As already mentioned in the literature, these subgrid models can conserve some but not all invariants. The SALT framework imposes helicity and circulation conservation in two and three dimensions and entrophy conservation in two dimensions, whereas LU dynamics strictly conserves kinetic energy. Yet, for a homogeneous unresolved turbulence, we have proved that – on average – LU dynamics also conserve the helicity and circulation.

This paper mainly focuses on numerical parameterization. We have formulated and described several parameterizations

<footnote n="1">¹Here, the buoyancy is not Gaussian. However, it is reasonable to believe that mean ±1.96 times the pointwise standard deviation remains a simple and convenient approximate metric to define an acceptable bias.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Buoyancy field (m s⁻²) (a), kinetic energy spectrum (m² s⁻²/(rad m⁻¹)) (b), and ADSD (m² s⁻¹/(rad m⁻¹)) (c) of w, at t = 105 d, ADSD of σ·Ḃ from the non-data-driven method (without a multiplicative constant); ADSD of σ·Ḃ from the data-driven method, with 2, 20, 200, 2000, and 8000 EOFs, and a slope of -5/3 (b); and a corresponding ADSD (c) in a black solid line. The two dashed vertical lines define the interval where coefficients Cw and rw are fitted. The dashed oblique line is the resulting fit (it is set to match Aw in the left bound of the interval).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Buoyancy fields (m s⁻²) (left), kinetic energy spectra (m² s⁻²/(rad m⁻¹)) (middle), and ADSDs (m² s⁻¹/(rad m⁻¹)) (right) of w, in blue, ADSDs of σ·Ḃ, in green, at t = 110 d, for (from top to bottom) the low-resolution (64²) projection of the reference deterministic high-resolution (512²) SQG dynamics; the low-resolution (64²) deterministic SQG dynamic; one realization of the low-resolution (64²) LU SQG dynamics (equivalent to SALT SQG) with the self-similar parameterization; and one realization of the low-resolution (64²) LU SQG dynamics (equivalent to SALT SQG) with data-driven parameterization.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Buoyancy fields (m s⁻²) (**a**, **b**), kinetic energy spectra (m² s⁻²/(rad m⁻¹)) (**c**, **d**), and ADSDs (m² s⁻¹/(rad m⁻¹)) (**e**, **f**) of w, at t = 110 d, for ensemble mean of the low-resolution (64²) LU SQG dynamics (equivalent to SALT SQG) with self-similar parameterization (**a**, **c**, **e**) and with data-driven parameterization (**b**, **d**, **f**).</figcaption></figure>

which apply to both SALT and LU frameworks. Parame-
terization for SALT or LU means choosing a spatial co-
variance for the unresolved small-scale velocity (i.e. choos-
ing σ). The stationary homogeneous self-similar parameter-
ization of Resseguier et al. (2017b) has been improved to
make it non-stationary and tuning-free. Spectral properties
are learned “on the fly” from the resolved large-scale ve-
locity as is common in the practice of large-eddy simula-
tion. A heterogeneous modulation of that parameterization
– based on the energy flux through scales – is also proposed.
This modulation naturally comes into play when considering
Kolmogorov-like energy cascades and spectra. Because the
energy flux quantifies the energy cascade, which is assumed
to be constant across scales, the modulated unresolved ve-
locity acts only where there is a large-scale energy cascade.
This parameterization improvement enables a better simula-
tion of straight strong fronts in SQG dynamics. However, the
usual convenient approximation produced in the Smagorin-
sky diagnosis of the energy flux – the dissipation – cannot
be used accurately here. We also recall the stationary het-
erogeneous data-driven parameterization method of Cotter
et al. (2019a, 2018). Here, small Lagrangian displacement
increments are computed at distinct resolutions from high-
resolution simulation outputs. The spectral decomposition of
the Lagrangian displacements covariance leads to a light and
convenient representation of the unresolved velocity statis-
tics: the Lagrangian displacement EOFs.

Finally, two tuning-free parameterizations have been nu-
merically compared: the new non-stationary homogeneous
self-similar parameterization for LU and the stationary het-
erogeneous data-driven one of Cotter et al. (2019a, 2018) for
SALT. The test case is a homogeneous and stationary forced
SQG turbulence whose LU and SALT versions are equiva-

lent. Single realizations of the LU test case are found to be
– at least – as good as the result of the corresponding deter-
ministic simulation at the same resolution. For both param-
eterizations, the ensemble forecasts are found to predict the
right amplitudes and positions of numerical errors. The un-
certainty skills of the ensemble forecasts remain impressively
good even with only 20 realizations.

5 Discussion

A lot of theoretical and practical questions remain about
SALT and LU schemes.

The numerical explorations of this paper have focused on
SQG dynamics because – except for the neglected advection
correction – SALT and LU SQG models coincide. This sim-
plification has enabled a clearer parameterization compari-
son. An interesting dual study would be the comparison of
distinct SALT and LU models but with a fixed parameteriza-
tion (i.e. fixed σ model).

For this purpose, the simplest dynamics to consider would
be 2D incompressible Navier–Stokes equations. The 3D or
quasi-geostrophic (QG) dynamics would be appropriate, al-
though they possess a larger parameter space and more costly
computation. The SALT version will conserve helicity and
circulation, while the LU one will conserve kinetic energy.
It would be a very interesting exercise to examine the iner-
tial cascades of energy, enstrophy, and potential enstrophy in
these systems and how they are affected by the SALT ver-
sus LU assumptions. Nonetheless, even an understanding of
these cascades is probably not sufficient to objectively con-
clude which framework is more appropriate in general, as the
answer will likely depend on the application, model resolu-
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Buoyancy fields ($m s^{-2}$) (left); kinetic energy spectra ($m^2 s^{-2} / (\text{rad m}^{-1})$) (middle); and ADSDs ($m^2 s^{-1} / (\text{rad m}^{-1})$) (right) of $w$, in blue, and ADSDs of $\sigma \dot{B}$, in green, at $t = 120$ d, for (from top to bottom) the low-resolution (64$^2$) projection of the reference deterministic high-resolution (512$^2$) SQG dynamics, the low-resolution (64$^2$) deterministic SQG dynamic, one realization of the low-resolution (64$^2$) LU SQG dynamics with self-similar parameterization, and one realization of the low-resolution (64$^2$) LU SQG dynamics (equivalent to SALT SQG) with data-driven parameterization.</figcaption></figure>

tion, etc. A variety of idealized and applied numerical studies will probably be necessary to gain insight into how to optimize in a variety of settings. And even for one specific situation, the quality and skill metrics of choice are not obvious. Even in the simple cases studied here, the usual second-order statistics are probably not sufficient to discriminate between these dynamics. Furthermore, it is important to note that the exact conservation of vorticity by the SALT scheme and energy by the LU scheme may be counterproductive in certain heterogeneous settings where the small-scale features are required to exhibit systematic property transport. One interesting example is the oceanic wind-driven gyre, whose circulation magnitude depends critically on the gyre and its turbulence to relocate the source of vorticity and energy from the wind to the regions where dissipation occurs. An eddy vorticity

transport across large-scale streamlines (i.e. affecting the interpretation of the Kelvin circulation theorem) is needed (Fox-Kemper and Pedlosky, 2004). A cross-streamline energy transport is also part of the system equilibration (Scott and Straub, 1998). In the real ocean, the eddies shed in the Agulhas retroflection are a classic example of organized small-scale transport (Gordon, 1985) as are the “eddy canons” that fire mesoscale eddies across the Antarctic Circumpolar Current (Hallberg and Gnanadesikan, 2006).

In the presence of turbulence heterogeneity, another distinction between SALT and LU models is the advection correction. It is also the single distinction between LU dynamics and Mikulevicius and Rozovskii (2004a). So far, it is not clear what constitutes a “large-scale velocity” in practice. Again, this probably depends on the situation. Even though
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Buoyancy fields (m s⁻²) (**a**, **b**), kinetic energy spectra (m² s⁻²/(rad m⁻¹)) (**c**, **d**), and ADSDs (m² s⁻¹/(rad m⁻¹)) (**e**, **f**) of w, at t = 120 d, for the ensemble mean of the low-resolution (64²) LU SQG dynamics (equivalent to SALT SQG) with self-similar parameterization (**a**, **c**, **e**) and with data-driven parameterization (**b**, **d**, **f**).</figcaption></figure>

Cotter et al. (2017) have provided some first theoretical clues, further explorations are needed. Numerical and experimental work will probably help if one identifies and studies simple heterogeneous flows meeting the main assumption of SALT and LU models: the velocity timescale gap.

The parameterizations for SALT and LU presented in this study (e.g. the ADSD method) could be naturally extended to physical-domain-based implementation, as is an important step in evaluating schemes for operational use, e.g. Pearson et al. (2017). Indeed, Sect. 2.1 details this tuning-free parameterization in the Fourier space only. A convenient dual of self-similar spectra (i.e. spectra scaling as ||*k*||⁶) in the physical space is the widely used Matérn covariance (Williams and Rasmussen, 2006; Lim and Teo, 2009; Lilly et al., 2017). Thus, the spatial filter $\tilde{\sigma}$ – appearing in $\sigma dB_t = \tilde{\sigma}^* dB_t$ – of a possible future 2D physical-domain-based ADSD parameterization would probably be built from the curl of a Matérn covariance.

Nevertheless, the homogeneity of the ADSD parameterization is a drawback. Heterogeneous solutions need to be developed. Unfortunately, the ADSD heterogeneous modulation of this paper remains expensive to compute, even making use of the Fourier space. Moreover, even though its third-order physical meaning and its natural association with SALT and LU is very instructive, its additive value for UQ is not clear.

The heterogeneous parameterization method of Cotter et al. (2019a, 2018) is valuable in the presence of some turbulence heterogeneities (e.g. heterogeneities induced by fixed boundary conditions). However, the practical usefulness of this parameterization remains debatable, because its heterogeneity is assumed to be stationary. Other ways for calibrating

the subgrid unresolved statistics using machine learning is also on-going research.

Many parameterizations for SALT and LU are now available. Regardless of their differences, the resulting stochastic dynamics and associated UQ skill seems to be relatively independent of the particular parameterization choice. Notwithstanding, we think that learning-free heterogeneous and non-stationary parameterizations may further improve SALT and LU UQ skills and/or enable smaller ensemble size and thus more efficient computation.

Finally, the main goal of all these studies is SALT- and LU-based data assimilation. Cotter et al. (2019b) have opened the way, using a particle filter and the heterogeneous parameterization method of Cotter et al. (2019a, 2018). Observed data from a dynamical system of O(10⁶) degrees of freedom were successfully assimilated into the parameterized system of O(10⁴) degrees of freedom. We expect that many other data assimilation studies will follow.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Normalized buoyancy bias absolute value, $|\hat{E}(b) - b^{\text{ref}}|$, (dimensionless) of the stochastic SQG model (left) and its estimation (1.96 times the standard deviation of the ensemble) for the data-driven method (middle) and the self-similar one (right) at a resolution of 64² at (from top to bottom) *t* = 105, 110, 115, and 120 d. The reference is the usual SQG model at a resolution of 512², adequately filtered and subsampled. The stochastic simulations and the reference have a common initial condition at *t* = 100 d. Points used for the local UQ analysis of Figs. 13 and 14 are highlighted in green.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Confidence interval at 95% (light purple) and at 50% (dark purple) along time, computed from the low-resolution stochastic SQG 200-member ensemble with the data-driven (a, b) and with the self-similar (c, d) parameterization, at the points (500, 250 km) (a, c) and (250, 500 km) (b, d). The high-resolution reference is superimposed in red.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Confidence interval at 95% (light purple) and at 50% (dark purple) along time, computed from the low-resolution stochastic SQG 20-member ensemble with the data-driven (a, b) and with the self-similar (c, d) parameterization, at the points (500, 250 km) (a, c) and (250, 500 km) (b, d). The high-resolution reference is superimposed in red.</figcaption></figure>
## Appendix A: Theoretical motivations

Here, we briefly highlight the similarities and differences between the location uncertainty model and the stochastic Lie transport model. To simplify the comparison, we work in the Stratonovich representation and restrict ourselves to the incompressible case.

### A1 Lagrangian path

From a large-scale under-resolved point of view, SALT and LU methods assume that the fluid velocity is partially uncorrelated in time. Hence, the position of a Lagrangian particle $X_t$ evolves in time according to

$$ \frac{dX_t}{dt} = w^*(X_t, t) + (\sigma \dot{B})(X_t, t), \quad (A1) $$

where $\sigma \dot{B}$ is time-uncorrelated and $w^*$ is the Stratonovich large-scale drift term. Formally, for a spatial domain $\Omega \subset \mathbb{R}^d$, the process $(t \mapsto B_t)$ is a cylindrical $\mathbf{I}_d(E^2(\Omega))^d$-Wiener process (see Da Prato and Zabczyk, 1992, and Prévôt and Röckner, 2007, for more information on infinite-dimensional Wiener processes and cylindrical $\mathbf{I}_d$-Wiener processes). Recently, Cotter et al. (2017) have rigorously shown that such a decomposition corresponds to the limit of a deterministic flow when the correlation time of the small-scale velocity goes to zero.

### A2 Notation correspondences

The Lagrangian equation can also be written with Itō notations as follows:

$$ X_{t+dt} - X_t = w(X_t, t)dt + \sigma(X_t, t)(B_{t+dt} - B_t). \quad (A2) $$

Note the difference in the drift term when compared with Eq. (A1). Table A1 summarizes the notations differences between SALT and LU both for Itō and Stratonovich notations.

### A3 Scalar transport

For a – possibly active – scalar-tracer denoted $q$, the SALT prescribes the same type of evolution equation as the LU models.

$$ \frac{Dq}{Dt}(\mathbf{x}, t) \triangleq \frac{d}{dt}(q(X_t, t))|_{X_t=x} = 0, \quad (A3) $$

where the material derivative, $\frac{Dq}{Dt}$, (in Stratonovich notations) is simply

$$ \frac{Dq}{Dt} = \frac{\partial}{\partial t} + \dot{\mathbf{x}}_t \cdot \nabla, \quad (A4) $$

with

$$ \dot{\mathbf{x}}_t(\mathbf{x}) \triangleq w^*(\mathbf{x}, t) + (\sigma \dot{B})(\mathbf{x}, t). \quad (A5) $$

In Itō form, the transport Eq. (A3) makes explicit the turbulent diffusion and the centred anti-symmetric multiplicative noise (Resseguier et al., 2017a).

We again highlight that this analysis holds for both SALT and LU approaches.

### A4 Euler models

Nonetheless, the incompressible stochastic transports of velocity and vorticity differ between LU and SALT. First, the SALT approach considers the transport – through a specific Lagrangian choice and up to some forcings – of the Stratonovich large-scale linear momentum, $\rho w^*$, whereas the LU Euler derivation assumes the transport of the Itō large-scale linear momentum, $\rho w$. Furthermore, due to its geometrical approach, the SALT Euler equation involves an additional term.

Specifically, the LU Euler equation with neither viscosity nor Coriolis force reads

$$ \frac{D}{Dt} \underbrace{w}_{\text{transport}} = -\frac{1}{\rho} \nabla p, \quad (A6) $$

Due to the transport of $\rho w$

whereas the SALT version is

$$ \frac{D}{Dt} \underbrace{w^*}_{\text{transport}} + \underbrace{\nabla \dot{\mathbf{x}}_t^T w^*}_{\text{Additional term}} = -\frac{1}{\rho} \nabla p. \quad (A7) $$

Due to the transport of $\rho w^*$

Unlike in the SALT model, the large-scale transported velocity of the LU Euler, $w$, differs from the large-scale transporting velocity, $w^*$. The latter implicitly appears in the stochastic transport operator, $\frac{D}{Dt}$, of both the SALT and LU equations. Thus, in the LU equations, we can see the correction
<table><caption>**Table A1.** Notation equivalences between LU and SALT approaches.</caption><thead><tr><th></th><th>LU</th><th>SALT</th></tr></thead><tbody><tr><td>Unresolved velocity</td><td>(σ&#x0307;)(x, t) = ∫<sub>Ω</sub>dz&#x0307;(x, z) dB<sub>t</sub>(z)/dt</td><td>∑<sub>p</sub>ξ<sub>p</sub>(x)dW<sub>p</sub>(t)/dt</td></tr><tr><td>Covariance</td><td>σ(x)σ<sup>T</sup>(y)/dt = ∫<sub>Ω</sub>dz&#x0307;(x, z)σ<sup>T</sup>(y, z)/dt</td><td>∑<sub>p</sub>ξ<sub>p</sub>(x)ξ<sub>p</sub><sup>T</sup>(y)/dt</td></tr><tr><td>Variance tensor</td><td>a(x) = σ(x)σ<sup>T</sup>(x) = ∫<sub>Ω</sub>dz&#x0307;(x, z)σ<sup>T</sup>(x, z)</td><td>∑<sub>p</sub>ξ<sub>p</sub>(x)ξ<sub>p</sub><sup>T</sup>(x)</td></tr><tr><td>Stratonovich drift</td><td>w*</td><td>u</td></tr><tr><td>Itō drift</td><td>w = w* + 1/2(∇·a)<sup>T</sup></td><td>u + 1/2(∇·(∑<sub>p</sub>ξ<sub>p</sub>ξ<sub>p</sub><sup>T</sup>))<sup>T</sup></td></tr></tbody></table>

$w^* - w = -\frac{1}{2}(\nabla \cdot a)^T$ as a modification of the large-scale advection. Note that this modification cancels for a homogeneous small-scale velocity. Otherwise, whether the Itō drift,

$$w(x, t) \stackrel{\Delta}{=} \mathbb{E} (X_{t+dt} - X_t | X_t = x) / dt, \quad (A8)$$

or the Stratonovich drift,

$$w^*(x, t) \stackrel{\Delta}{=} \mathbb{E} (X_{t+dt/2} - X_{t-dt/2} | X_t = x) / dt, \quad (A9)$$

should be (randomly) transported is still an open question. A related question is how to interpret the large-scale velocity derived from observations or from numerical simulations. Depending on the situation, this velocity may better correspond to the Itō drift or to the Stratonovich drift. Cotter et al. (2017) may help answer this question from a theoretical perspective.

On top of this large-scale advection difference between SALT and LU modelling, the SALT additional term, $\nabla \dot{x}_t^T w^*$, has major consequences on the dynamics invariants, as explained in Sect. A6.

## A5 Other dynamical models

The differences between LU and SALT momentum equations – Eqs. (A6) and (A7) respectively – resemble the distinctions between other classical and geophysical fluid dynamic models. Table A2 recalls some of these classical models in SALT and LU stochastic formulations. The LU Euler equations can be found in Mémin (2014) without noise and in Resseguier (2017) including noise. The vorticity equation in each case is easily derived by taking the curl of the Euler equation. The LU QG equations under moderate influence of turbulence has been derived by Resseguier et al. (2017b) with Itō notation. In Table A2, we have written the same equation in Stratonovich notation. The SALT Euler, vorticity and quasi-geostrophic equations are derived by Holm (2015).

Both SALT and LU assumes the stochastic transport of tracers (e.g. temperature and salinity). For QG models, it implies the transport of buoyancy, $b$, at the surface ($z = 0$). Therefore, up to the drift correction $w^* - w$ in the streamfunction definition, SALT and potential vorticity (PV) coincide. In the SQG framework, the PV is assumed to be zero

in the ocean interior ($z < 0$). This leads to the same relationship between buoyancy and streamfunction, $\psi$, in LU and SALT dynamics. Therefore, in the Sect. 3 of this paper, we will choose this model to compare two parameterization (i.e. choice of $\sigma$) for SALT and LU methods in a common ground.

## A6 Invariants

To interpret these results which apply more generally, mutatis mutandis the KE is conserved by the LU scheme and its variants (Resseguier et al., 2017a), while the enstrophy, its generalizations (e.g. PV), the circulation, and the helicity are conserved by the SALT scheme and its variants (Holm, 2015). In the specific case of homogeneous small-scale velocity, the LU dynamics also conserve circulation and helicity in average (see Appendix A7 for circulation mean conservation; the proof of helicity mean conservation is similar). However, in the homogeneous case, this stochastic dynamics increases enstrophy mean and its generalizations (Resseguier et al., 2017b), while SALT models increase KE mean (Resseguier, 2017).

These distinctions between conservation of vorticity in the SALT approach and conservation of energy in the LU approach would be of critical importance when choosing a model for 2D or QG turbulence, as Kraichnan (1967) and Charney (1971) show that the conservation of energy and enstrophy lead to turbulence typified by an inverse energy cascade at large scales and a direct (potential) enstrophy cascade at small scales (Fox-Kemper and Menemenlis, 2008). In SQG, Blumen (1982) shows that a dual cascade results from conservation of depth-integrated energy and available potential energy on level boundaries. The former is not singled out for special treatment in the SALT or LU SQG formulation, but the latter is exactly conserved in both formulations.

## A7 Kelvin theorem under location uncertainty

The Kelvin theorem describes the variation of the circulation, which is defined as follows for LU dynamics as follows:

$$\Gamma = \oint_{C(t)} w(x, t) \cdot dl = \oint_{C(0)} w(X_t(l_0), t)^T J_M dl_0, \quad (A10)$$

where $l = X_t(l_0)$ is the Lagrangian path labelled by the initial position $l_0$, $C(t) = X_t(C(0))$ is a material loop at time
<table><caption>**Table A2.** Some LU and SALT dynamics models with Stratonovich notations.</caption><thead><tr><th>Equations</th><th>LU</th><th>SALT</th></tr></thead><tbody><tr><td>3D Euler equation in a rotating frame (with 3D ∇)</td><td>ρ(D/Dt + f ×)w + f × σdot = (ρg - ∇p), D<sub>t</sub>ρ = 0.</td><td>ρ(D/Dt + ∇dot<sub>t</sub><sup>T</sup> + f ×)w* + f × σdot = (ρg - ∇p), D<sub>t</sub>ρ = 0.</td></tr><tr><td>Incompressible 3D vorticity equation (with 3D ∇)</td><td>ω = ∇ × w, Dω/Dt = (ω · ∇)dot - ∑<sub>q=1</sub><sup>d</sup> ∇(dot<sub>t</sub>)<sub>q</sub> × ∇w<sub>q</sub>, ∇ · w* = 0, ∇ · σ = 0.</td><td>ω* = ∇ × w*, Dω*/Dt = (ω* · ∇)dot, ∇ · w* = 0, ∇ · σ = 0.</td></tr><tr><td>Incompressible 2D vorticity equation (with 2D ∇)</td><td>ζ = ∇⊥ · w, Dζ/Dt = -tr(∇⊥dot<sub>t</sub><sup>T</sup> ∇w<sup>T</sup>), ∇ · w* = 0, ∇ · σ = 0.</td><td>ζ* = ∇⊥ · w*, Dζ*/Dt = 0, ∇ · w* = 0, ∇ · σ = 0.</td></tr><tr><td>Quasi-geostrophic equations (with 2D ∇ and Δ)</td><td>w = ∇⊥ψ, ∇ · w* = 0, ∇ · σ = 0, PV = Δψ + f + (f<sub>0</sub>/N)<sup>2</sup>∂<sub>z</sub><sup>2</sup>ψ, f<sub>0</sub>∂<sub>z</sub>ψ = b, DPV/Dt = (w* - w) · ∇f dt - tr(∇⊥dot<sub>t</sub><sup>T</sup> ∇w<sup>T</sup>), for z < 0 Db/Dt = 0 for z = 0.</td><td>w* = ∇⊥ψ, ∇ · σ = 0, PV = Δψ + f + (f<sub>0</sub>/N)<sup>2</sup>∂<sub>z</sub><sup>2</sup>ψ, f<sub>0</sub>∂<sub>z</sub>ψ = b, DPV/Dt = 0, for z < 0 Db/Dt = 0 for z = 0.</td></tr><tr><td>Surface quasi-geostrophic equation (with 2D ∇ and Δ)</td><td>w = ∇⊥ψ, ∇ · w* = 0, ∇ · σ = 0, (-Δ)<sup>1/2</sup>ψ = 1/Nb, Db/Dt = 0.</td><td>w* = ∇⊥ψ, ∇ · σ = 0, (-Δ)<sup>1/2</sup>ψ = 1/Nb, Db/Dt = 0.</td></tr></tbody></table>

t and $J_M = J_M(l_0, t) = [\nabla X_t^T]^T(l_0)$ is the Jacobian matrix of the flow. The time differentiation of the circulation involves the time variation of the Jacobian matrix ($dJ_M = [\nabla \dot{x}_t^T]^T J_M$) as follows:

$$ \begin{aligned} \frac{d\Gamma}{dt} &= \oint_{C(0)} \left( \frac{Dw}{Dt}^T J_M + w^T [\nabla \dot{x}_t^T]^T J_M \right) dl_0 \\ &= \oint_{C(t)} \left( -\frac{1}{\rho} \nabla p + [\nabla \dot{x}_t^T] w \right) \cdot dl \\ &= \oint_{C(t)} (\nabla \dot{x}_t^T w) \cdot dl. \end{aligned} \quad (A11) $$

This equation is the equivalent of the Reynolds transport theorem for vorticity Fox-Kemper and Pedlosky (2004) but in the stochastic framework. This noise term is a priori not certified, since it is a Stratonovich noise. However, in the homogeneous case, its ensemble mean cancels. Indeed, using Itō notations, we simply need to compute a quadratic cross variation. This calculus is possible in Lagrangian coordinates (i.e. when functions are composed by $x \mapsto X_t$) by noticing that $Dw/Dt$ is a term in $dt$, $D\sigma/Dt = (\dot{x}_t \cdot \nabla)\sigma$,

and $DJ_M = dJ_M = [\nabla \dot{x}_t^T]^T J_M$ as follows:

$$ \begin{aligned} \frac{d}{dt} \mathbb{E}\{\Gamma\} &= \frac{1}{dt} \mathbb{E} \oint_{C(0)} \left( w^T [\nabla(\sigma \circ dB_t)^T]^T J_M \right) dl_0, \\ &= \frac{1}{2} \mathbb{E} \oint_{C(0)} \sum_k \frac{D}{Dt} \langle w^T [\nabla \sigma_{\bullet k}^T]^T J_M dl_0, (B_t)_k \rangle, \\ &= \frac{1}{2} \mathbb{E} \oint_{C(0)} \sum_k w^T \left( (\sigma_{\bullet k} \cdot \nabla) + [\nabla \sigma_{\bullet k}^T]^T \right) \\ & \quad [\nabla \sigma_{\bullet k}^T]^T J_M dl_0, \\ &= \frac{1}{2} \mathbb{E} \oint_{C(t)} \sum_{pqk} w_q \nabla \cdot \partial_p \underbrace{(\sigma_{qk} \sigma_{\bullet k})}_{\text{cst.}} dl_p, \\ &= 0. \end{aligned} \quad (A12) $$

Therefore, the mean LU circulation is conserved in the homogeneous case.
## Appendix B: Effective resolution and inertial range

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

$$F(\|k\|) = \exp(-\nu t \|k\|^{2p}). \qquad (B1)$$

If the hyperviscosity $\nu$ is well chosen, we may expect that at the Shannon resolution $\pi/\Delta x = \kappa_M$, only 10% of the energy is left by the filter, i.e.

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

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 (B3)$$

The ratio $\kappa_M/\kappa_M$ can then be derived from Eqs. (B1), (B2), and (B3).
*Data availability.* No experimental data were collected nor used in this study. All data used were numerically generated from model implementations. The complete code can be found at the GitHub repository https://github.com/vressegu/LU_SALT_SelfSim (last access: 27 March 2020, GitHub, 2020).

*Author contributions.* VR and WP decided to perform this SALT versus LU numerical study. BFK obtained the funding and supervised the project. VR and BFK designed the ADSD methodology and its heterogeneous variants. VR developed the code and performed the simulations. WP helped to reproduce the data-driven parameterization of Cotter et al. (2019a). All authors prepared the paper.

*Competing interests.* One of the author is currently employed by a private company named SCALIAN.

*Acknowledgements.* We would like to warmly thank Darryl D. Holm, Dan Crisan, Colin Cotter, Igor Shevchenko, Bertrand Chapron, and Etienne Mémin for helpful discussions and for having enabled the collaborations which have lead to this paper.

*Financial support.* This research has been supported by the National Science Foundation Division of Ocean Sciences (grant no. 1350795), the Office of Naval Research Office of Naval Research Global (grant no. N00014-17-1-2963), and the Engineering and Physical Sciences Research Council (grant no. EP/N023781/1).

*Review statement.* This paper was edited by Wansuo Duan and reviewed by two anonymous referees.

## References

Bachman, S. D., Fox-Kemper, B., and Pearson, B.: A scale-aware subgrid model for quasi-geostrophic turbulence, J. Geophys. Res.-Oceans, 122, 1529–1554, 2017.

Berner, J., Achatz, U., Batte, L., Camara, A. D. L., 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. N., Penland, C., Storch, J.-S. V., Sakradzija, M., Weniger, M., Weisheimer, A., Williams, P. D., and Yano, J.-I.: Stochastic Parameterization: Towards a new view of Weather and Climate Models, B. Am. Meteorol. Soc., 98, 565–588, 2017.

Blumen, W.: Uniform potential vorticity flow: part I. theory of wave interactions and two-dimensional turbulence, J. Atmos. Sci., 35, 774–783, 1978.

Blumen, W.: Wave-Interactions in Quasi-Geostrophic Uniform Potential Vorticity Flow, J. Atmos. Sci., 39, 2388–2396, 1982.

Cai, S., Mémin, E., Dérian, P., and Xu, C.: Motion estimation under location uncertainty for turbulent fluid flows, Exp. Fluids, 59, 8, https://doi.org/10.1007/s00348-017-2458-z, 2018.

Callies, J., Flierl, G., Ferrari, R., and Fox-Kemper, B.: The role of mixed-layer instabilities in submesoscale turbulence, J. Fluid Mech., 788, 5–41, 2016.

Chapron, B., Dérian, P., Mémin, E., and Resseguier, V.: Large-scale flows under location uncertainty: a consistent stochastic framework, Q. J. Roy. Meteor. Soc., 144, 251–260, 2018.

Charney, J. G.: Geostrophic Turbulence, J. Atmos. Sci., 28, 1087–1095, 1971.

Chekroun, M. D., Kondrashov, D., and Ghil, M.: Predicting stochastic systems by noise sampling, and application to the El Niño-Southern Oscillation, P. Natl. Acad. Sci. USA, 108, 11766–11771, 2011.

Constantin, P., Majda, A., and Tabak, E.: Formation of strong fronts in the 2-D quasigeostrophic thermal active scalar, Nonlinearity, 7, 1495–1533, 1994.

Constantin, P., Nie, Q., and Schörghofer, N.: Front formation in an active scalar equation, Phys. Rev. E, 60, 2858, https://doi.org/10.1103/PhysRevE.60.2858, 1999.

Constantin, P., Lai, M., Sharma, R., Tseng, Y., and Wu, J.: New numerical results for the surface quasi-geostrophic equation, J. Sci. Comput., 50, 1–28, 2012.

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

Cotter, C., Crisan, D., Holm, D., Pan, W., and Shevchenko, I.: Numerically Modelling Stochastic Lie Transport in Fluid Dynamics, SIAM Multiscale Modeling & Simulation, 17, 192–232, https://doi.org/10.1137/18M1167929, 2019a.

Cotter, C., Crisan, D., Holm, D. D., Pan, W., and Shevchenko, I.: A Particle Filter for Stochastic Advection by Lie Transport (SALT): A case study for the damped and forced incompressible 2D Euler equation, arXiv preprint arXiv:1907.11884, 2019b.

Cotter, C. J., Gottwald, G. A., and Holm, D. D.: Stochastic partial differential fluid equations as a diffusive limit of deterministic Lagrangian multi-time dynamics, P. Roy. Soc. A-Math. Phy., 473, 20170388, https://doi.org/10.1098/rspa.2017.0388, 2017.

Crisan, D. and Lang, O.: Well-posedness for 2D Euler equations with stochastic Lie transport noise, arXiv 1907.00451, 2018.

Crisan, D., Flandoli, F., and Holm, D. D.: Solution properties of a 3D stochastic Euler fluid equation, J. Nonlinear Sci., 29, 813–870, 2019.

Da Prato, G. and Zabczyk, J.: Stochastic equations in infinite dimensions, Cambridge University Press, https://doi.org/10.1017/CBO9781107295513, 1992.

Falkovich, G., Gawedzki, K., and Vergassola, M.: Particles and fields in fluid turbulence, Rev. Mod. Phys., 73, 913, https://doi.org/10.1103/RevModPhys.73.913, 2001.

Ferrari, R. and Nikurashin, M.: Suppression of Eddy Diffusivity across Jets in the Southern Ocean, J. Phys. Oceanogr., 40, 1501–1519, 2010.

Flandoli, F.: The interaction between noise and transport mechanisms in PDEs, Milan J. Math., 79, 543–560, 2011.

Fox-Kemper, B.: New Frontiers in Operational Oceanography, chap.: Notions for the Motions of the Oceans, GODAE Ocean-View, 27–73, https://doi.org/10.17125/gov2018, 2018.

Fox-Kemper, B. and Menemenlis, D.: Can Large Eddy Simulation Techniques Improve Mesoscale-Rich Ocean Models?, in: Ocean Modeling in an Eddying Regime, edited by: Hecht, M. and Ha-
sumi, H., AGU Geophysical Monograph Series, 177, 319–338, 2008.

Fox-Kemper, B. and Pedlosky, J.: Wind-driven barotropic gyre I: Circulation control by eddy vorticity fluxes to an enhanced removal region, J. Mar. Res., 62, 169–193, 2004.

Frank, J. and Gottwald, G.: Stochastic homogenization for an energy conserving multi-scale toy model of the atmosphere, Physica D, 254, 46–56, 2013.

Frisch, U.: Turbulence: the legacy of AN Kolmogorov, Cambridge university press, https://doi.org/10.1017/CBO978113917066670666, 1995.

Gawedzky, K. and Kupiainen, A.: Anomalous scaling of the passive scalar, Phys. Rev. Lett., 75, 3834–3837, 1995.

Gay-Balmaz, F. and Holm, D. D.: Stochastic Geometric Models with Non-stationary Spatial Correlations in Lagrangian Fluid Flows, J. Nonlinear Sci., 28, 873–904, https://doi.org/10.1007/s00332-017-9431-0, 2018.

GitHub repository: https://github.com/vressegu/LU_SALT_SelfSim, last access: 27 March 2020.

Gordon, A. L.: Indian-Atlantic transfer of thermocline water at the Agulhas retroflection, Science, 227, 1030–1033, 1985.

Gottwald, G. and Melbourne, I.: Homogenization for deterministic maps and multiplicative noise, P. Roy. Soc. A-Math. Phy., 469, 20130201, https://doi.org/10.1098/rspa.2013.0201, 2013.

Hallberg, R. and Gnanadesikan, A.: The role of eddies in determining the structure and response of the wind-driven Southern Hemisphere overturning: Results from the Modeling Eddies in the Southern Ocean (MESO) project, J. Phys. Oceanogr., 36, 2232–2252, 2006.

Held, I., Pierrehumbert, R., Garner, S., and Swanson, K.: Surface quasi-geostrophic dynamics, J. Fluid Mech., 282, 1–20, 1995.

Holm, D. D.: Variational principles for stochastic fluid dynamics, P. Roy. Soc. A-Math. Phy., 471, 20140963, https://doi.org/10.1098/rspa.2014.0963, 2015.

Isern-Fontanet, J., Chapron, B., Lapeyre, G., and Klein, P.: Potential use of microwave sea surface temperatures for the estimation of ocean currents, Geophys. Res. Lett., 33, L24608, https://doi.org/10.1029/2006GL027801, 2006.

Isern-Fontanet, J., Lapeyre, G., Klein, P., Chapron, B., and Hecht, M. W.: Three-dimensional reconstruction of oceanic mesoscale currents from surface information, J. Geophys. Res., 113, C09005, https://doi.org/10.1029/2007JC004692, 2008.

Janjić, T., McLaughlin, D., Cohn, S. E., and Verlaan, M.: Conservation of mass and preservation of positivity with ensemble-type Kalman filter algorithms, Mon. Weather Rev., 142, 755–773, 2014.

Keating, S., Smith, S., and Kramer, P.: Diagnosing lateral mixing in the upper ocean with virtual tracers: Spatial and temporal resolution dependence, J. Phys. Oceanogr., 41, 1512–1534, 2011.

Klein, P., Isern-Fontanet, J., Lapeyre, G., Roullet, G., Danioux, E., Chapron, B., Le Gentil, S., and Sasaki, H.: Diagnosis of vertical velocities in the upper ocean from high resolution sea surface height, Geophys. Res. Lett., 36, L12603, https://doi.org/10.1029/2009GL038359, 2009.

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

Kolmogorov, A. N.: The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Cr Acad. Sci. URSS, 30, 301–305, 1941.

Kraichnan, R.: Small-scale structure of a scalar field convected by turbulence, Phys. Fluids, 11, 945–963, 1968.

Kraichnan, R.: Anomalous scaling of a randomly advected passive scalar, Phys. Rev. Lett., 72, 1016, https://doi.org/10.1103/PhysRevLett.72.1016, 1994.

Kraichnan, R. H.: Inertial Ranges in Two-Dimensional Turbulence, Phys. Fluids, 16, 1417–1423, 1967.

LaCasce, J. and Mahadevan, A.: Estimating subsurface horizontal and vertical velocities from sea-surface temperature, J. Mar. Res., 64, 695–721, 2006.

Lapeyre, G.: Surface quasi-geostrophy, Fluids, 2, 7, https://doi.org/10.3390/fluids2010007, 2017.

Lapeyre, G. and Klein, P.: Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory, J. Phys. Oceanogr., 36, 165–176, 2006a.

Lapeyre, G. and Klein, P.: Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory, J. Phys. Oceanogr., 36, 165–176, 2006b.

Leith, C.: Atmospheric predictability and two-dimensional turbulence, J. Atmos. Sci., 28, 145–161, 1971.

Lilly, J. M., Sykulski, A. M., Early, J. J., and Olhede, S. C.: Fractional Brownian motion, the Matérn process, and stochastic modeling of turbulent dispersion, Nonlin. Processes Geophys., 24, 481–514, https://doi.org/10.5194/npg-24-481-2017, 2017.

Lim, S. and Teo, L.: Generalized Whittle–Matérn random field as a model of correlated fluctuations, J. Phys. A-Math. Theor., 42, 105202, https://doi.org/10.1088/1751-8113/42/10/105202, 2009.

Lu, F., Lin, K. K., and Chorin, A. J.: Data-based stochastic model reduction for the Kuramoto–Sivashinsky equation, Physica D, 340, 46–57, 2017.

Majda, A. and Kramer, P.: Simplified models for turbulent diffusion: Theory, numerical modelling, and physical phenomena, Phys. Rep., 314, 237–574, 1999.

Majda, A., Timofeyev, I., and Vanden-Eijnden, E.: Models for stochastic climate prediction, P. Natl. Acad. Sci. USA, 96, 14687–14691, 1999.

Majda, A. J.: Statistical energy conservation principle for inhomogeneous turbulent dynamical systems, P. Natl. Acad. Sci. USA, 112, 8937–8941, 2015.

Margolin, L. G., Rider, W. J., and Grinstein, F. F.: Modeling turbulent flow with implicit LES, J. Turbul., 7, N15, https://doi.org/10.1080/14685240500331595, 2006.

Mémin, E.: Fluid flow dynamics under location uncertainty, Geophys. Astro. Fluid, 108, 119–146, https://doi.org/10.1080/03091929.2013.836190, 2014.

Mikulevicius, R. and Rozovskii, B.: Stochastic Navier–Stokes equations for turbulent flows, SIAM J. Math. Anal., 35, 1250–1310, 2004a.

Mikulevicius, R. and Rozovskii, B.: Stochastic Navier–Stokes equations for turbulent flows, SIAM J. Math. Anal., 35, 1250–1310, 2004b.

Mitchell, L. and Gottwald, G.: Data assimilation in slow fast systems using homogenized climate models, J. Atmos. Sci., 69, 1359–1377, 2012.

Orszag, S.: Analytical theories of turbulence, J. Fluid Mech., 41, 363–386, 1970.
Pearson, B. and Fox-Kemper, B.: Log-Normal Turbulence Dissipation in Global Ocean Models, Phys. Rev. Lett., 120, 094501, https://doi.org/10.1103/PhysRevLett.120.094501, 2018.

Pearson, B., Fox-Kemper, B., Bachman, S. D., and Bryan, F. O.: Evaluation of scale-aware subgrid mesoscale eddy models in a global eddy-rich model, Ocean Model., 115, 42–58, 2017.

Penland, C.: A Stochastic Approach to Nonlinear Dynamics: A Review (extended version of the article – “Noise Out of Chaos and Why it Won’t Go Away”), B. Am. Meteorol. Soc., 84, 925–925, 2003.

Prévôt, C. and Röckner, M.: A concise course on stochastic partial differential equations, in: Lecture Notes in Mathematics, 1905, 148 pp., Springer-Verlag, Berlin Heidelberg, 2007.

Resseguier, V.: Mixing and fluid dynamics under location uncertainty, PhD thesis, Université Rennes 1, 2017.

Resseguier, V., Mémin, E., and Chapron, B.: Geophysical flows under location uncertainty, Part I: Random transport and general models, Geophys. Astro. Fluid, 111, 149–176, https://doi.org/10.1080/03091929.2017.1310210, 2017a.

Resseguier, V., Mémin, E., and Chapron, B.: Geophysical flows under location uncertainty, Part II: Quasi-geostrophic models and efficient ensemble spreading, Geophys. Astro. Fluid, 111, 177–208, https://doi.org/10.1080/03091929.2017.1312101, 2017b.

Resseguier, V., Mémin, E., and Chapron, B.: Geophysical flows under location uncertainty, Part III: SQG and frontal dynamics under strong turbulence, Geophys. Astro. Fluid, 11, 209–227, https://doi.org/10.1080/03091929.2017.1312102, 2017c.

San, O., Staples, A. E., and Iliescu, T.: Approximate deconvolution large eddy simulation of a stratified two-layer quasigeostrophic ocean model, Ocean Model., 63, 1–20, 2013.

Sapsis, T. and Majda, A.: A statistically accurate modified quasilinear Gaussian closure for uncertainty quantification in turbulent dynamical systems, Physica D, 252, 34–45, 2013.

Sardeshmukh, P. D. and Sura, P.: Reconciling non-Gaussian climate statistics with linear dynamics, J. Climate, 22, 1193–1207, 2009.

Sardeshmukh, P. D., Compo, G. P., and Penland, C.: Need for caution in interpreting extreme weather statistics, J. Climate, 28, 9166–9187, 2015.

Scott, R. B. and Straub, D. N.: Small Viscosity Behavior of a Homogeneous, Quasigeostrophic, Ocean Circulation Model, J. Mar. Res., 56, 1225–1258, 1998.

Smagorinsky, J.: General circulation experiments with the primitive equation: I. The basic experiment, Mon. Weather Rev., 91, 99–165, 1963.

Tandeo, P., Autret, E., Chapron, B., Fablet, R., and Garello, R.: SST spatial anisotropic covariances from METOP-AVHRR data, Remote Sens. Environ., 141, 144–148, 2014.

Vallis, G.: Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation, Cambridge University Press, https://doi.org/10.1017/9781107588417, 2006.

Voosen, P.: The Earth Machine, Science, 361, 344–347, 2018.

Wang, J., Flierl, G. R., LaCasce, J. H., McClean, J. L., and Mahadevan, A.: Reconstructing the ocean’s interior from surface data, J. Phys. Oceanogr., 43, 1611–1626, 2013.

Williams, C. K. and Rasmussen, C. E.: Gaussian processes for machine learning, in: Adaptive Computation and Machine Learning series, The MIT Press, 272 pp., 2006.