Research article 16 Apr 2020
Research article  16 Apr 2020
Datadriven versus selfsimilar parameterizations for stochastic advection by Lie transport and location uncertainty
 ^{1}Lab, SCALIAN DS, Espace Nobel, 2 Allée de Becquerel, 35700 Rennes, France
 ^{2}Department of Mathematics, Imperial College London, London SW7 2AZ, UK
 ^{3}Department of Earth, Environmental and Planetary Sciences (DEEPS), Brown University, Providence, RI 02912, USA
 ^{4}Institute at Brown for Environment and Society (IBES), Brown University, Providence, RI 02912, USA
 ^{1}Lab, SCALIAN DS, Espace Nobel, 2 Allée de Becquerel, 35700 Rennes, France
 ^{2}Department of Mathematics, Imperial College London, London SW7 2AZ, UK
 ^{3}Department of Earth, Environmental and Planetary Sciences (DEEPS), Brown University, Providence, RI 02912, USA
 ^{4}Institute at Brown for Environment and Society (IBES), Brown University, Providence, RI 02912, USA
Correspondence: Valentin Resseguier (valentin.resseguier@scalian.com)
Hide author detailsCorrespondence: Valentin Resseguier (valentin.resseguier@scalian.com)
Stochastic subgrid parameterizations enable ensemble forecasts of fluid dynamic systems and ultimately accurate data assimilation (DA). Stochastic advection by Lie transport (SALT) and models under location uncertainty (LU) are recent and similar physically based stochastic schemes. SALT dynamics conserve helicity, whereas LU models conserve kinetic energy (KE). After highlighting general similarities between LU and SALT frameworks, this paper focuses on their common challenge: the parameterization choice. We compare uncertainty quantification skills of a stationary heterogeneous datadriven parameterization and a nonstationary homogeneous selfsimilar parameterization. For stationary, homogeneous surface quasigeostrophic (SQG; QG) turbulence, both parameterizations lead to highquality ensemble forecasts. This paper also discusses a heterogeneous adaptation of the homogeneous parameterization targeted at a better simulation of strong straight buoyancy fronts.
Geophysical fluid dynamics and other turbulent flows cover a wide range of spatial and temporal scales (see, e.g. FoxKemper, 2018). However, numerical simulations and observations are limited in the range of scales that can be simulated 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 regularizations (e.g. nugget in optimal interpolation, Tandeo et al., 2014) are incorporated into numerical simulations and measurement processes. The limited resolution of these numerical simulations introduces unknown errors which are continuously 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, observations (data) are integrated into the model state of the numerical simulation limited by imperfect physics and imperfect initial conditions. DA proceeds by alternating between forecast and analysis cycles. In each analysis cycle, observations are combined with the results from a prediction model (the forecast) to produce a socalled analysis. The analysis is expected to be more accurate than either predictions based on physical models without incorporating observations or predictions from physicsfree 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 precision and sampling; and wellknown uncertainty quantification (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 quantify because often a “perfect model”, exact solution, or direct 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), stochastic 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 FoxKemper (2018) highlight the relevance of skewsymmetric 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 wellposedness 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 timehomogeneous subgrid parameterization and extended to the inhomogeneous case in GayBalmaz 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 wellposedness result for the 3D SALT Euler equations is proved in Crisan et al. (2019), and a global wellposedness 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 twolayer 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 largescale component, w, and a timeuncorrelated component, $\mathit{\sigma}\dot{\mathit{B}}=\mathit{\sigma}\mathrm{d}{\mathit{B}}_{t}/\mathrm{d}t$. The latter is Gaussian (conditionally on some largescale 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 infinitedimensional linear operator, σ, to a ddimensional space–time white noise, $\dot{\mathit{B}}$. Holm and coauthors rely on a different but equivalent notation. It directly explicates the spectral decomposition of the timeuncorrelated velocity spatial covariance. Specifically, the smallscale velocity reads $\mathit{\sigma}\left(\mathit{x}\right)\dot{\mathit{B}}={\sum}_{p}{\mathit{\xi}}_{p}\left(\mathit{x}\right)\mathrm{d}{W}_{p}\left(t\right)/\mathrm{d}t$, where (W_{p})_{p} is a set of independent 1D Brownian motions and (ξ_{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 σ or (ξ_{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 the 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 nondatadriven choices. Specifically, we propose a nonstationary and tuningfree improvement of the work of Resseguier et al. (2017b) based on selfsimilar assumptions and possible heterogeneous modulation. Links with the Smagorinsky subgrid parameterization and nonlinear energy fluxes are also discussed.

Section 2.2 discusses a datadriven 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 quasigeostrophic (SQG) models, their numerical setup, and simulation parameters. These are followed by detailed discussions of uncertainty quantification skills in Sect. 3.4.
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 quasigeostrophic 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; Callies 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; IsernFontanet 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.1 Parametric and selfsimilar model for the unresolved velocity
In this section, we propose parametric nondatadriven 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 parameterfree 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 Smagorinskylike subgrid parameterizations.
2.1.1 Parameterfree and nonstationary homogeneous model
We first define the statistics of the smallscale 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, $\mathrm{tr}\left(\mathit{a}\right)=\mathbb{E}{\u2225\mathit{\sigma}\mathrm{d}{\mathit{B}}_{t}\u2225}^{\mathrm{2}}/\mathrm{d}t$, 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, Kubotype formulas may be seen as what physicists call absolute diffusivities. More generally, since the variance of a timecontinuous 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 timeuncorrelated velocity. Thus, keeping a spectral approach, we define – for any spatiotemporal field – an absolute diffusivity spectral density 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:
is the eddy turnover time at the scale 1∕κ and v_{κ} is at characteristic velocity at this scale. Accordingly, we have
If in addition we assume a KE selfsimilar distribution,
and we obtain
where $r=(s+\mathrm{3})/\mathrm{2}$.
We aim at defining the unresolved velocity ADSD from the largescale velocity. For this purpose, we will assume the selfsimilar model (Eq. 4) is valid at all spatial scales. At each time step, we compute the ADSD of the largescale 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 timedependent because they depend on w which is. More precisely, we estimate the coefficients C_{w} and r_{w} in a wavenumber interval which approximately represents a inertial range of fully resolved scales (i.e. before the spectrum rolloff). 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
Therefore, the unresolved ADSD is chosen to compensate the resolved ADSD rolloff – introduced by the deterministic subgrid parameterization – at small scales. Specifically, the unresolved ADSD is set to
As previously, ${f}_{\mathrm{BP}}^{\mathrm{2}}$ is a bandpass filter between κ_{m} and κ_{M}. In practice, we set κ_{M} to the theoretical resolution,
and κ_{m} to the effective resolution, which is estimated as follows:
where p is the order of the Laplacian used as a deterministic subgrid tensor (i.e. $\frac{Db}{Dt}=\mathit{\nu}(\mathrm{\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 κ_{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 bandpass 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 κ_{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 ${\mathbf{A}}_{\mathit{\sigma}\dot{\mathit{B}}}$ through Eq. (6) at large wavenumbers – larger than κ_{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}_{\mathit{\sigma}\dot{\mathit{B}}}$ from the largescale velocity statistics. From that ADSD, it is possible to sample the unresolved velocity, as explained in the following.
A 2D homogeneous divergencefree smallscale velocity can be constructed by filtering a 1D white noise $\dot{B}$ with a 2D divergencefree filter $\stackrel{\mathrm{\u02d8}}{\mathit{\sigma}}={\mathbf{\nabla}}^{\perp}{\stackrel{\mathrm{\u02d8}}{\mathit{\psi}}}_{\mathit{\sigma}}$ (Resseguier et al., 2017b).
where ⋆ denotes a spatial convolution. This velocity field can be easily sampled in Fourier space.
where ${\widehat{\mathrm{d}B}}_{t}$ is the spatial Fourier transform of dB_{t} and $\frac{\mathrm{d}{B}_{t}}{\sqrt{\mathrm{\Delta}t}}$ is a discrete scalar whitenoise process of unit variance in space and time. Thus, the smallscale velocity is defined by the Fourier transform of the streamfunction kernel, $\widehat{{\stackrel{\mathrm{\u02d8}}{\mathit{\psi}}}_{\mathit{\sigma}}}$.
In order to link the unresolved ADSD to the kernel $\stackrel{\mathrm{\u02d8}}{\mathit{\sigma}}$ which defines the unresolved velocity, we note that
where μ(Ω) is the surface of the spatial domain Ω and θ_{k} is the angle of the wave vector k. From Eqs. (6)–(11), we can finally express the unresolved velocity as follows:
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 (timecorrelated) unresolved velocity ADSD, whereas the KE spectra of the simulated and true unresolved velocities differ. Indeed, the true unresolved velocity correlation time spectral distribution τ(κ) is not restricted to the time step dt.
The tuningfree 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 nonstationary 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 nonstationary and heterogeneous, and so regionbyregion and seasonbyseason tuning is impossible to do with quality checks in place.
We believe that nonstationarity may yield a fair, yet idealized, comparison. So, we here compare the two parameterizations in a nonstationary 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 selfsimilar 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 nonlinear processes. In a similar way, Figs. 3 and 4 show that our stochastic dynamics enable a morerealistic eddy distribution than deterministic simulations. However, a homogeneous smallscale velocity may also perturb the tracer isolines which should remain still (i.e. which remain still in highresolution 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 realworld 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 smallscale velocity is required. Note that the heterogeneity discussed here needs to be nonstationary 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, κ, cannot depend on the position, x, the constant, C_{w}, and/or the spectrum slope, r, should do. A spatially varying spectrum slope is probably difficult to estimate. Hence, we restrict ourselves to a spatially varying constant, C_{w}, and a spatially homogeneous spectrum slope. The constant may also varies with the time and the wavenumber. According to the Kolmogorov theory (e.g. Frisch, 1995) and Eq. (3):
where ϵ_{F} is the energy flux through the spatial scales and $p=\mathrm{1}/\mathrm{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∕κ toward smaller scales and can be computed as follows:
where q is the transported (up to possible source terms) quantity, ${g}_{\mathit{\kappa}}^{<}$ is the lowpass filtered version of g (setting to zero the Fourier modes of g which have frequencies larger than κ), and $\stackrel{\mathrm{\u203e}}{\u2022}$ 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 thirdorder moment. It is very important because it describes the cascade of the flow by nonlinear 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 velocity heterogeneities. We simply modulate the unresolved ADSD (Eq. 6) by the heterogeneous ratio ${\mathit{\u03f5}}_{\mathrm{F}}^{p}/\stackrel{\mathrm{\u203e}}{{\mathit{\u03f5}}_{\mathrm{F}}^{p}}$, averaged over the resolved inertial range wavenumbers.
where $\mathcal{P}={\mathbf{I}}_{\mathrm{d}}{\mathrm{\Delta}}^{\mathrm{1}}\mathbf{\nabla}{\mathbf{\nabla}}^{T}$ is the projector onto the space of freedivergence functions. This parameterization is physically meaningful, since, locally in space, a stronger direct cascade at large scales (larger ϵ_{F} and thus larger C_{w}) suggests that the unresolved velocity (large) should maintain this cascade by folding smallerscale tracer structures. Furthermore, considering that the energy flux is a thirdorder 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 divergencefree velocity and the ensuing properties (e.g. energy conservation), the modulated velocity is projected onto the space of freedivergence functions, using the operator 𝒫. Because of that we do not consider the advection correction ${\mathit{w}}^{*}\mathit{w}=\frac{\mathrm{1}}{\mathrm{2}}(\mathbf{\nabla}\mathbf{\cdot}\mathbf{a}{)}^{T}$ of the LU formalism. Indeed, here the variance tensor has the simple form $\mathbf{a}=\frac{\mathrm{1}}{d}\text{tr}\left(\mathbf{a}\right){\mathbf{I}}_{\mathrm{d}}$. As such, the advection correction is a gradient field and is hence removed by the projection onto the space of freedivergence 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 secondorder statistics (e.g. KE spectra) cannot distinguish between eddies and filaments. Bispectra may overcome this drawback, since they express threepoint 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. FoxKemper and Menemenlis, 2008; Bachman et al., 2017) provide a path to developing deterministic and dissipative scaleaware subgrid models. Typically, these models result in a Laplacian dissipation which involves a heterogeneous eddy diffusivity or viscosity coefficient, ν_{Sm}. Aside from their heuristic theoretical justification, these Smagorinskytype subgrid terms are formally equivalent to the turbulent dissipation of our stochastic model: $\mathbf{\nabla}\mathbf{\cdot}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathbf{a}\mathbf{\nabla}q\right)$ (Resseguier et al., 2017a), where a=2ν_{Sm}I_{d} and q is a transported quantity. This similarity suggests that a Smagorinskytype 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,
where ϵ_{D} is the dissipation, a secondorder moment related to the molecular or turbulent diffusion. To develop a Smagorinskytype model, the resolved flux of the energylike conserved invariant is equated to the dissipation of the energy invariant such that
From there, one can obtain an eddy diffusivity or viscosity coefficient proportional to $\Vert \mathbf{\nabla}q{\Vert}^{h}$. For an SQG flow, the exponent is $h=\mathrm{1}/\mathrm{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 secondorder 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 higherorder statistics. Figure 4 (middle right) illustrates this statement. The following unresolved velocity parameterization is used there:
Along sharp straight fronts, the dissipation will be larger. Accordingly, the associated modulation $\Vert \mathbf{\nabla}b{\Vert}^{\mathrm{1}/\mathrm{4}}$ enhances the stochastic folding where one would need it to weaken.
2.2 Datadriven 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, (ξ_{k})_{k}, involved in the unresolved velocity definition,
If we denote λ_{k} as the L^{2} norm of those EOFs, their normalized versions $({\mathit{\xi}}_{k}={\mathit{\xi}}_{k}/{\mathit{\lambda}}_{k}{)}_{k}$ are the eigenvectors of the selfadjoint operator,
defined by the smallscale velocity covariance [σ(x)σ^{T}(y)∕dt]. (λ_{k}∕dt)_{k} is the set of eigenvalues of this operator.
The datadriven 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 infinitedimensional Brownian motion, σ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 highresolution velocity, v, on a fine meshgrid (512^{2}),

a lowresolution velocity, $\stackrel{\mathrm{\u203e}}{\mathit{v}}$, on a coarse meshgrid (64^{2}). This velocity field is a spatially lowpassfiltered version of f.
Then, two types of Lagrangian path are defined:

a “highresolution flow”, X_{ij}(t_{0},t), defined by the highresolution velocity, u:
$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\mathrm{d}{\mathit{X}}_{ij}}{\mathrm{d}t}}\left(t\right)=\mathit{u}\left(t,{\mathit{X}}_{ij}({t}_{\mathrm{0}},t)\right)\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{\mathit{X}}_{ij}({t}_{\mathrm{0}},{t}_{\mathrm{0}})={\mathit{x}}_{ij},and\end{array}$$with ${\mathit{x}}_{ij}=\left(i\mathrm{\Delta}x,j\mathrm{\Delta}y\right)$.

a “lowresolution flow”, ${\stackrel{\mathrm{\u203e}}{\mathit{X}}}_{ij}({t}_{\mathrm{0}},t)$, defined by the lowresolution velocity $\stackrel{\mathrm{\u203e}}{\mathit{u}}$:
$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\mathrm{d}{\stackrel{\mathrm{\u203e}}{\mathit{X}}}_{ij}}{\mathrm{d}t}}({t}_{\mathrm{0}},t)=\stackrel{\mathrm{\u203e}}{\mathit{u}}\left(t,{\stackrel{\mathrm{\u203e}}{\mathit{X}}}_{ij}({t}_{\mathrm{0}},t)\right)\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{\stackrel{\mathrm{\u203e}}{\mathit{X}}}_{ij}({t}_{\mathrm{0}},{t}_{\mathrm{0}})={\mathit{x}}_{ij}.\end{array}$$
2.2.2 Candidate for the increment realization
In order to estimate the EOFs, (ξ_{i})_{i}, involved in Eq. (A5), we assume that we can observe increments
We will interpret the following residual flow increments as a realization of the above:
2.2.3 Preprocessing
The increments are supposed to be centred and divergence free (as $\mathbf{\nabla}\mathbf{\cdot}{\mathit{\xi}}_{i}=\mathrm{0}$, ∀i). Therefore, after computing the residual flow increments, ${\left({\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}\mathit{X}}}_{ij}^{m}\right)}_{m}$, they are centred,
with the estimator
and projected onto the space of divergencefree functions,
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).
In order to properly define the EOFs, we must add the following orthogonal constraint:
Finally, after estimating the (ξ_{k})_{k} offline, the ensemble forecast can generated online with Eq. (22).
3.1 Surface quasigeostrophic 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,
where $\frac{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, $\widehat{\mathit{\psi}}$, is related to the Fourier transform of the buoyancy, $\widehat{b}$, by the following SQG relationship:
where N is the stratification and 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
in the SALT context, and $\mathit{w}={\mathbf{\nabla}}^{\perp}\mathit{\psi}={\mathbf{\nabla}}^{\perp}(\mathrm{\Delta}{)}^{\mathrm{1}/\mathrm{2}}(b/N)$, in the LU one. Besides, in the LU SQG, the advecting drift w^{⋆} has to be divergencefree. We enforce this constraint by projecting the drift correction $({\mathit{w}}^{*}\mathit{w})$ onto the space of freedivergence functions. So, we have
where $\mathcal{P}={\mathbf{I}}_{\mathrm{d}}{\mathrm{\Delta}}^{\mathrm{1}}\mathbf{\nabla}{\mathbf{\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 $\mathrm{0}={\text{PV}}_{\mathrm{slt}}={\mathbf{\nabla}}^{\perp}\mathbf{\cdot}{\mathit{w}}^{*}+f+{\left({f}_{\mathrm{0}}/N\right)}^{\mathrm{2}}{\partial}_{z}^{\mathrm{2}}\mathit{\psi}$, ${\mathit{w}}^{*}={\mathbf{\nabla}}^{\perp}\mathit{\psi}$ and $b={f}_{\mathrm{0}}{\partial}_{z}\mathit{\psi}$. For the LU dynamics, $\mathrm{0}={\text{PV}}_{\mathrm{lu}}={\mathbf{\nabla}}^{\perp}\mathbf{\cdot}\mathit{w}+f+{\left({f}_{\mathrm{0}}/N\right)}^{\mathrm{2}}{\partial}_{z}^{\mathrm{2}}\mathit{\psi}$, $\mathit{w}={\mathbf{\nabla}}^{\perp}\mathit{\psi}$ and $b={f}_{\mathrm{0}}{\partial}_{z}\mathit{\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 $({\mathit{w}}^{*}\mathit{w})$. So, we simulate the stochastic transport (Eq. 34) coupled with the SQG relation (Eq. 36). The unresolved velocity statistics encoded in σ are specified either by the selfsimilar method of Sect. 2.1.1 or by the datadriven method of Sect. 2.2.
3.1.3 Flow simulation
We perform a highresolution simulation (512^{2} spatial grid) of the deterministic SQG model with the following initial condition:
The source term, S, involves a hyperviscosity, a linear drag, and an additive stationary forcing such that
The parameters of the simulations are summed up in Table 1. The influence of the initial condition remains for about a 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 datadriven method are learned. From t=100 d to t=130 d, the deterministic highresolution simulation is used as a reference. In this time interval, several stochastic lowresolution (64^{2}) simulations are performed and compared. These simulations are initialized at t=100 d with the reference highresolution simulation projected at low resolution (i.e. keeping only the Fourier modes associated with a coarseresolution 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, ΔT. We first describe the convergence with the number of snapshots and then – as a remark – the convergence with the Lagrangian advection time, Δ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 nonstationary 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, ΔT, is a critical parameter for the definition of the EOF. Indeed, for small advection time, ΔT, the length of an increment of a Lagrangian flow path is proportional to the velocity and to the advection time. This is the socalled ballistic regime (Falkovich et al., 2001; Vallis, 2006). In particular, residual flow increments squared norm would be proportional to ΔT^{2}; the variance tensor estimator would be proportional to ΔT; and the estimated EOFs would be proportional to $\sqrt{\mathrm{\Delta}T}$. With a larger advection time, Δ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{\mathrm{\Delta}T}$; i.e. the Lagrangian velocity acts as a white noise in time. This is the socalled diffusive regime. Figure 6 illustrates this convergence with the average tensor ${a}_{\mathrm{0}}=\frac{\mathrm{1}}{d}\stackrel{\mathrm{\u203e}}{\mathit{tr}\left(\mathbf{a}\right)}$ with ΔT. In this paper, the Lagrangian advection time ΔT is computed from a CFL (Courant–Friedrichs–Lewy condition) at the coarse resolution of 64^{2} – 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 spatially aliased, since a large part of those increments lives at small spatial scales but are spatially sampled on a coarse spatial grid, x_{ij}. Figure 7 confirms this idea. If all the $\mathrm{2}\times {\mathrm{64}}^{\mathrm{2}}\mathrm{1}=\mathrm{8191}$ EOFs are considered, the ADSD of the unresolved velocity reveals a strong spatial aliasing.
Both the datadriven and the nondata 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 datadriven method involves some large spatial scales. One could think that these largescale components would disappear if the flow X and $\stackrel{\mathrm{\u203e}}{\mathit{X}}$ are integrated in a Eulerian way rather than in a Lagrangian way. However, with the advection time, Δ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 datadriven model (see Sect. 2.2) or the parametric and selfsimilar model (see Sect. 2.1). For the datadriven 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 lowresolution (64^{2}) projection of the reference deterministic highresolution (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 lowresolution deterministic SQG simulations and the lowresolution (64^{2}) projection of the reference deterministic highresolution simulation. As already pointed out by Resseguier et al. (2017b) for freedecaying turbulence, the realizations of the LU dynamics are no worse than a lowresolution deterministic simulation. For the shortterm simulations, our stochastic subgrid parameterizations have often weak improvements on the lowresolution 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×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×128. In Fig. 8, the resolution is coarser (64×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 location across simulations of larger vortices in Fig. 10, while the filaments between the vortices have lost all coherence. In Fig. 11, the larger, morecoherent features persist in the ensemble mean of the coarseresolution 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^{2}). The first ensemble is generated with the datadriven model for the unresolved velocity (see Sect. 2.2), while the second ensemble is generated with the parametric and selfsimilar model for the unresolved velocity (see Sect. 2.1). For the datadriven 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 datadriven and the selfsimilar parameterizations achieve this goal everywhere and for every time. The pointwise biases of the datadriven method, of the selfsimilar 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 datadriven method. In Fig. 12, the represented biases and error estimations (mean ±1.96 times the pointwise standard deviation^{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 freeturbulence flow (Resseguier et al., 2017b). This close analysis 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 nonlinear 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×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 ensemblebased confidence intervals at 95 % and at 50 % and compare them with the reference value. Figure 13 display 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 200ensemblemember 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.
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 largescale velocity and an unresolved timeuncorrelated 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 enstrophy 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 which apply to both SALT and LU frameworks. Parameterization for SALT or LU means choosing a spatial covariance for the unresolved smallscale velocity (i.e. choosing σ). The stationary homogeneous selfsimilar parameterization of Resseguier et al. (2017b) has been improved to make it nonstationary and tuningfree. Spectral properties are learned “on the fly” from the resolved largescale velocity as is common in the practice of largeeddy simulation. A heterogeneous modulation of that parameterization – based on the energy flux through scales – is also proposed. This modulation naturally comes into play when considering Kolmogorovlike energy cascades and spectra. Because the energy flux quantifies the energy cascade, which is assumed to be constant across scales, the modulated unresolved velocity acts only where there is a largescale energy cascade. This parameterization improvement enables a better simulation of straight strong fronts in SQG dynamics. However, the usual convenient approximation produced in the Smagorinsky diagnosis of the energy flux – the dissipation – cannot be used accurately here. We also recall the stationary heterogeneous datadriven parameterization method of Cotter et al. (2019a, 2018). Here, small Lagrangian displacement increments are computed at distinct resolutions from highresolution simulation outputs. The spectral decomposition of the Lagrangian displacements covariance leads to a light and convenient representation of the unresolved velocity statistics: the Lagrangian displacement EOFs.
Finally, two tuningfree parameterizations have been numerically compared: the new nonstationary homogeneous selfsimilar parameterization for LU and the stationary heterogeneous datadriven 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 equivalent. Single realizations of the LU test case are found to be – at least – as good as the result of the corresponding deterministic simulation at the same resolution. For both parameterizations, the ensemble forecasts are found to predict the right amplitudes and positions of numerical errors. The uncertainty skills of the ensemble forecasts remain impressively good even with only 20 realizations.
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 simplification has enabled a clearer parameterization comparison. An interesting dual study would be the comparison of distinct SALT and LU models but with a fixed parameterization (i.e. fixed σ model).
For this purpose, the simplest dynamics to consider would be 2D incompressible Navier–Stokes equations. The 3D or quasigeostrophic (QG) dynamics would be appropriate, although 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 inertial cascades of energy, enstrophy, and potential enstrophy in these systems and how they are affected by the SALT versus LU assumptions. Nonetheless, even an understanding of these cascades is probably not sufficient to objectively conclude which framework is more appropriate in general, as the answer will likely depend on the application, model resolution, 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 secondorder 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 smallscale features are required to exhibit systematic property transport. One interesting example is the oceanic winddriven 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 largescale streamlines (i.e. affecting the interpretation of the Kelvin circulation theorem) is needed (FoxKemper and Pedlosky, 2004). A crossstreamline 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 smallscale transport (Gordon, 1985) as are the “eddy cannons” 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 “largescale velocity” in practice. Again, this probably depends on the situation. Even though 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 physicaldomainbased implementation, as is an important step in evaluating schemes for operational use, e.g. Pearson et al. (2017). Indeed, Sect. 2.1 details this tuningfree parameterization in the Fourier space only. A convenient dual of selfsimilar spectra (i.e. spectra scaling as $\Vert \mathit{k}{\Vert}^{\mathit{\alpha}}$) 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 $\stackrel{\mathrm{\u02d8}}{\mathit{\sigma}}$ – appearing in $\mathit{\sigma}\mathrm{d}{\mathit{B}}_{t}={\stackrel{\mathrm{\u02d8}}{\mathit{\sigma}}}^{*}\mathrm{d}{B}_{t}$ – of a possible future 2D physicaldomainbased 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 thirdorder 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 ongoing 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 learningfree heterogeneous and nonstationary 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 LUbased 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^{6}) degrees of freedom were successfully assimilated into the parameterized system of O(10^{4}) degrees of freedom. We expect that many other data assimilation studies will follow.
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 largescale underresolved 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
where $\mathit{\sigma}\dot{\mathit{B}}$ is timeuncorrelated and w^{*} is the Stratonovich largescale drift term. Formally, for a spatial domain Ω⊂ℝ^{d}, the process (t↦B_{t}) is a cylindrical ${{\mathbf{I}}_{\mathrm{d}}}_{{\left({\mathcal{L}}^{\mathrm{2}}\left(\mathrm{\Omega}\right)\right)}^{d}}$Wiener process (see Da Prato and Zabczyk, 1992, and Prévôt and Röckner, 2007, for more information on infinitedimensional Wiener processes and cylindrical 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 smallscale velocity goes to zero.
A2 Notation correspondences
The Lagrangian equation can also be written with Itō notations as follows:
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 – scalartracer denoted q, the SALT prescribes the same type of evolution equation as the LU models.
where the material derivative, $\frac{Dq}{\mathrm{d}t}$, (in Stratonovich notations) is simply
with
In Itō form, the transport Eq. (A3) makes explicit the turbulent diffusion and the centred antisymmetric 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 largescale linear momentum, ρw^{*}, whereas the LU Euler derivation assumes the transport of the Itō largescale linear momentum, ρ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
whereas the SALT version is
Unlike in the SALT model, the largescale transported velocity of the LU Euler, w, differs from the largescale 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 ${\mathit{w}}^{*}\mathit{w}=\frac{\mathrm{1}}{\mathrm{2}}{\left(\mathbf{\nabla}\mathbf{\cdot}\mathit{a}\right)}^{T}$ as a modification of the largescale advection. Note that this modification cancels for a homogeneous smallscale velocity. Otherwise, whether the Itō drift,
or the Stratonovich drift,
should be (randomly) transported is still an open question. A related question is how to interpret the largescale 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 largescale advection difference between SALT and LU modelling, the SALT additional term, $\mathbf{\nabla}{\dot{\mathit{x}}}_{t}^{T}{\mathit{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 quasigeostrophic 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 ${\mathit{w}}^{*}\mathit{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, ψ, 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 σ) 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 smallscale 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 (FoxKemper and Menemenlis, 2008). In SQG, Blumen (1982) shows that a dual cascade results from conservation of depthintegrated 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:
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 t and ${\mathbf{J}}_{M}={\mathbf{J}}_{M}({\mathit{l}}_{\mathrm{0}},t)={\left[\mathbf{\nabla}{\mathit{X}}_{t}^{T}\right]}^{T}\left({\mathit{l}}_{\mathrm{0}}\right)$ is the Jacobian matrix of the flow. The time differentiation of the circulation involves the time variation of the Jacobian matrix ($\mathrm{d}{\mathbf{J}}_{M}={\left[\mathbf{\nabla}{\dot{\mathit{x}}}_{t}^{T}\right]}^{T}{\mathit{J}}_{M}$) as follows:
This equation is the equivalent of the Reynolds transport theorem for vorticity FoxKemper and Pedlosky (2004) but in the stochastic framework. This noise term is a priori not centred, 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↦X_{t}) by noticing that $\frac{D\mathit{w}}{Dt}$ is a term in dt, $\frac{D\mathit{\sigma}}{Dt}=({\dot{\mathit{x}}}_{t}\mathbf{\cdot}\mathbf{\nabla})\mathit{\sigma}$, and $\frac{D{\mathit{J}}_{M}}{Dt}=\frac{\mathrm{d}{\mathit{J}}_{M}}{\mathrm{d}t}={\left[\mathbf{\nabla}{\dot{\mathit{x}}}_{t}^{T}\right]}^{T}{\mathit{J}}_{M}$ as follows:
Therefore, the mean LU circulation is conserved in the homogeneous case.
Let us assume the simulated evolution law is ${D}_{t}b=\mathit{\nu}(\mathrm{\Delta}{)}^{p}b\mathrm{d}t$. The deterministic subgrid model $\mathit{\nu}(\mathrm{\Delta}{)}^{p}b$ acts, in a finite time t, as a lowpass filter. In Fourier space, this filter is
If the hyperviscosity ν is well chosen, we may expect that at the Shannon resolution $\mathit{\pi}/\mathrm{\Delta}x={\mathit{\kappa}}_{M}$, only 10 % of the energy is left by the filter, i.e.
A ratio smaller than 10 % may lead to an overdamped simulation. Moreover, the precise value of this ratio does not influence much our final estimate.
We may define the effective resolution as the scale κ=κ_{m}, where the deterministic subgrid model influence is negligible. There, we may expect the filter to be equal to 95 %, i.e.
The ratio κ_{m}∕κ_{M} can then be derived from Eqs. (B1), (B2), and (B3).
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).
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 datadriven parameterization of Cotter et al. (2019a). All authors prepared the paper.
One of the author is currently employed by a private company named SCALIAN.
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.
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. N000141712963), and the Engineering and Physical Sciences Research Council (grant no. EP/N023781/1).
This paper was edited by Wansuo Duan and reviewed by two anonymous referees.
Bachman, S. D., FoxKemper, B., and Pearson, B.: A scaleaware subgrid model for quasigeostrophic turbulence, J. Geophys. Res.Oceans, 122, 1529–1554, 2017. a, b, c
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. a
Blumen, W.: Uniform potential vorticity flow: part I. theory of wave interactions and twodimensional turbulence, J. Atmos. Sci., 35, 774–783, 1978. a
Blumen, W.: WaveInteractions in QuasiGeostrophic Uniform Potential Vorticity Flow, J. Atmos. Sci., 39, 2388–2396, 1982. a
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/s003480172458z, 2018. a
Callies, J., Flierl, G., Ferrari, R., and FoxKemper, B.: The role of mixedlayer instabilities in submesoscale turbulence, J. Fluid Mech., 788, 5–41, 2016. a
Chapron, B., Dérian, P., Mémin, E., and Resseguier, V.: Largescale flows under location uncertainty: a consistent stochastic framework, Q. J. Roy. Meteor. Soc., 144, 251–260, 2018. a
Charney, J. G.: Geostrophic Turbulence, J. Atmos. Sci., 28, 1087–1095, 1971. a
Chekroun, M. D., Kondrashov, D., and Ghil, M.: Predicting stochastic systems by noise sampling, and application to the El NiñoSouthern Oscillation, P. Natl. Acad. Sci. USA, 108, 11766–11771, 2011. a
Constantin, P., Majda, A., and Tabak, E.: Formation of strong fronts in the 2D quasigeostrophic thermal active scalar, Nonlinearity, 7, 1495–1533, 1994. a, b, c, d, e
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. a, b, c, d
Constantin, P., Lai, M., Sharma, R., Tseng, Y., and Wu, J.: New numerical results for the surface quasigeostrophic equation, J. Sci. Comput., 50, 1–28, 2012. a, b, c, d
Cotter, C., Crisan, D., Holm, D. D., Pan, W., and Shevchenko, I.: Modelling uncertainty using circulationpreserving stochastic transport noise in a 2layer quasigeostrophic model, arXiv preprint arXiv:1802.05711, 2018. a, b, c, d, e, f
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. a, b, c, d, e, f, g, h, i, j, k
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. a, b
Cotter, C. J., Gottwald, G. A., and Holm, D. D.: Stochastic partial differential fluid equations as a diffusive limit of deterministic Lagrangian multitime dynamics, P. Roy. Soc. AMath. Phy., 473, 20170388, https://doi.org/10.1098/rspa.2017.0388, 2017. a, b, c, d
Crisan, D. and Lang, O.: Wellposedness for 2D Euler equations with stochastic Lie transport noise, arXiv 1907.00451, 2018. a
Crisan, D., Flandoli, F., and Holm, D. D.: Solution properties of a 3D stochastic Euler fluid equation, J. Nonlinear Sci., 29, 813–870, 2019. a
Da Prato, G. and Zabczyk, J.: Stochastic equations in infinite dimensions, Cambridge University Press, https://doi.org/10.1017/CBO9781107295513, 1992. a
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. a, b, c
Ferrari, R. and Nikurashin, M.: Suppression of Eddy Diffusivity across Jets in the Southern Ocean, J. Phys. Oceanogr., 40, 1501–1519, 2010. a
Flandoli, F.: The interaction between noise and transport mechanisms in PDEs, Milan J. Math., 79, 543–560, 2011. a
FoxKemper, B.: New Frontiers in Operational Oceanography, chap.: Notions for the Motions of the Oceans, GODAE OceanView, 27–73, https://doi.org/10.17125/gov2018, 2018. a
FoxKemper, B. and Menemenlis, D.: Can Large Eddy Simulation Techniques Improve MesoscaleRich Ocean Models?, in: Ocean Modeling in an Eddying Regime, edited by: Hecht, M. and Hasumi, H., AGU Geophysical Monograph Series, 177, 319–338, 2008. a, b
FoxKemper, B. and Pedlosky, J.: Winddriven barotropic gyre I: Circulation control by eddy vorticity fluxes to an enhanced removal region, J. Mar. Res., 62, 169–193, 2004. a, b
Frank, J. and Gottwald, G.: Stochastic homogenization for an energy conserving multiscale toy model of the atmosphere, Physica D, 254, 46–56, 2013. a
Frisch, U.: Turbulence: the legacy of AN Kolmogorov, Cambridge university press, https://doi.org/10.1017/CBO978113917066670666, 1995. a, b, c
Gawedzky, K. and Kupiainen, A.: Anomalous scaling of the passive scalar, Phys. Rev. Lett., 75, 3834–3837, 1995. a
GayBalmaz, F. and Holm, D. D.: Stochastic Geometric Models with Nonstationary Spatial Correlations in Lagrangian Fluid Flows, J. Nonlinear Sci., 28, 873–904, https://doi.org/10.1007/s0033201794310, 2018. a
GitHub repository: https://github.com/vressegu/LU_SALT_SelfSim, last access: 27 March 2020. a
Gordon, A. L.: IndianAtlantic transfer of thermocline water at the Agulhas retroflection, Science, 227, 1030–1033, 1985. a
Gottwald, G. and Melbourne, I.: Homogenization for deterministic maps and multiplicative noise, P. Roy. Soc. AMath. Phy., 469, 20130201, https://doi.org/10.1098/rspa.2013.0201, 2013. a
Hallberg, R. and Gnanadesikan, A.: The role of eddies in determining the structure and response of the winddriven Southern Hemisphere overturning: Results from the Modeling Eddies in the Southern Ocean (MESO) project, J. Phys. Oceanogr., 36, 2232–2252, 2006. a
Held, I., Pierrehumbert, R., Garner, S., and Swanson, K.: Surface quasigeostrophic dynamics, J. Fluid Mech., 282, 1–20, 1995. a
Holm, D. D.: Variational principles for stochastic fluid dynamics, P. Roy. Soc. AMath. Phy., 471, 20140963, https://doi.org/10.1098/rspa.2014.0963, 2015. a, b, c, d
IsernFontanet, 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. a
IsernFontanet, J., Lapeyre, G., Klein, P., Chapron, B., and Hecht, M. W.: Threedimensional reconstruction of oceanic mesoscale currents from surface information, J. Geophys. Res., 113, C09005, https://doi.org/10.1029/2007JC004692, 2008. a
Janjić, T., McLaughlin, D., Cohn, S. E., and Verlaan, M.: Conservation of mass and preservation of positivity with ensembletype Kalman filter algorithms, Mon. Weather Rev., 142, 755–773, 2014. a
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. a, b
Klein, P., IsernFontanet, 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. a
Klyatskin, V.: Stochastic equations through the eye of the physicist: Basic concepts, exact results and asymptotic approximations, Elsevier, 2005. a, b
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. a
Kraichnan, R.: Smallscale structure of a scalar field convected by turbulence, Phys. Fluids, 11, 945–963, 1968. a
Kraichnan, R.: Anomalous scaling of a randomly advected passive scalar, Phys. Rev. Lett., 72, 1016, https://doi.org/10.1103/PhysRevLett.72.1016, 1994. a
Kraichnan, R. H.: Inertial Ranges in TwoDimensional Turbulence, Phys. Fluids, 16, 1417–1423, 1967. a
LaCasce, J. and Mahadevan, A.: Estimating subsurface horizontal and vertical velocities from seasurface temperature, J. Mar. Res., 64, 695–721, 2006. a
Lapeyre, G.: Surface quasigeostrophy, Fluids, 2, 7, https://doi.org/10.3390/fluids2010007, 2017. a
Lapeyre, G. and Klein, P.: Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory, J. Phys. Oceanogr., 36, 165–176, 2006a. a
Lapeyre, G. and Klein, P.: Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory, J. Phys. Oceanogr., 36, 165–176, 2006b. a
Leith, C.: Atmospheric predictability and twodimensional turbulence, J. Atmos. Sci, 28, 145–161, 1971. a
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/npg244812017, 2017. a
Lim, S. and Teo, L.: Generalized Whittle–Matérn random field as a model of correlated fluctuations, J. Phys. AMath. Theor., 42, 105202, https://doi.org/10.1088/17518113/42/10/105202, 2009. a
Lu, F., Lin, K. K., and Chorin, A. J.: Databased stochastic model reduction for the Kuramoto–Sivashinsky equation, Physica D, 340, 46–57, 2017. a
Majda, A. and Kramer, P.: Simplified models for turbulent diffusion: Theory, numerical modelling, and physical phenomena, Phys. Rep., 314, 237–574, 1999. a
Majda, A., Timofeyev, I., and VandenEijnden, E.: Models for stochastic climate prediction, P. Natl. Acad. Sci. USA, 96, 14687–14691, 1999. a
Majda, A. J.: Statistical energy conservation principle for inhomogeneous turbulent dynamical systems, P. Natl. Acad. Sci. USA, 112, 8937–8941, 2015. a
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. a
Mémin, E.: Fluid flow dynamics under location uncertainty, Geophys. Astro. Fluid, 108, 119–146, https://doi.org/10.1080/03091929.2013.836190, 2014. a, b
Mikulevicius, R. and Rozovskii, B.: Stochastic Navier–Stokes equations for turbulent flows, SIAM J. Math. Anal., 35, 1250–1310, 2004a. a
Mikulevicius, R. and Rozovskii, B.: Stochastic Navier–Stokes equations for turbulent flows, SIAM J. Math. Anal., 35, 1250–1310, 2004b. a
Mitchell, L. and Gottwald, G.: Data assimilation in slow fast systems using homogenized climate models, J. Atmos. Sci., 69, 1359–1377, 2012. a
Orszag, S.: Analytical theories of turbulence, J. Fluid Mech., 41, 363–386, 1970. a
Pearson, B. and FoxKemper, B.: LogNormal Turbulence Dissipation in Global Ocean Models, Phys. Rev. Lett., 120, 094501, https://doi.org/10.1103/PhysRevLett.120.094501, 2018. a
Pearson, B., FoxKemper, B., Bachman, S. D., and Bryan, F. O.: Evaluation of scaleaware subgrid mesoscale eddy models in a global eddyrich model, Ocean Model., 115, 42–58, 2017. a
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. a, b
Prévôt, C. and Röckner, M.: A concise course on stochastic partial differential equations, in: Lecture Notes in Mathematics, 1905, 148 pp., SpringerVerlag, Berlin Heidelberg, 2007. a
Resseguier, V.: Mixing and fluid dynamics under location uncertainty, PhD thesis, Université Rennes 1, 2017. a, b
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. a, b, c, d
Resseguier, V., Mémin, E., and Chapron, B.: Geophysical flows under location uncertainty, Part II: Quasigeostrophic models and efficient ensemble spreading, Geophys. Astro. Fluid, 111, 177–208, https://doi.org/10.1080/03091929.2017.1312101, 2017b. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q
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. a
San, O., Staples, A. E., and Iliescu, T.: Approximate deconvolution large eddy simulation of a stratified twolayer quasigeostrophic ocean model, Ocean Model., 63, 1–20, 2013. a
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. a
Sardeshmukh, P. D. and Sura, P.: Reconciling nonGaussian climate statistics with linear dynamics, J. Climate, 22, 1193–1207, 2009. a
Sardeshmukh, P. D., Compo, G. P., and Penland, C.: Need for caution in interpreting extreme weather statistics, J. Climate, 28, 9166–9187, 2015. a
Scott, R. B. and Straub, D. N.: Small Viscosity Behavior of a Homogeneous, Quasigeostrophic, Ocean Circulation Model, J. Mar. Res., 56, 1225–1258, 1998. a
Smagorinsky, J.: General circulation experiments with the primitive equation: I. The basic experiment, Mon. Weather Rev., 91, 99–165, 1963. a, b
Tandeo, P., Autret, E., Chapron, B., Fablet, R., and Garello, R.: SST spatial anisotropic covariances from METOPAVHRR data, Remote Sens. Environ., 141, 144–148, 2014. a
Vallis, G.: Atmospheric and oceanic fluid dynamics: fundamentals and largescale circulation, Cambridge University Press, https://doi.org/10.1017/9781107588417, 2006. a, b, c
Voosen, P.: The Earth Machine, Science, 361, 344–347, 2018. a
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. a
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. a
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.
 Abstract
 Introduction
 Parameterization of SALT and LU models: the statistics of the unresolved velocity
 Numerical simulations and uncertainty quantification
 Conclusion
 Discussion
 Appendix A: Theoretical motivations
 Appendix B: Effective resolution and inertial range
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Parameterization of SALT and LU models: the statistics of the unresolved velocity
 Numerical simulations and uncertainty quantification
 Conclusion
 Discussion
 Appendix A: Theoretical motivations
 Appendix B: Effective resolution and inertial range
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References