A note on ‘Toward a stochastic parameterization of ocean mesoscale eddies’

Porta Mana & Zanna (2014) recently proposed a subgrid-scale parameterization for eddy-permitting quasigeostrophic models. In this model the large-scale fluid is represented as a non-Newtonian viscoelastic medium, with a subgridstress closure that involves the Lagrangian derivative of large-scale quantities. This note derives this parameterization, including the nondimensional proportionality coefficient, using only two statistical assumptions: that the subgrid-scale term is locally homogeneous and decorrelates rapidly in space. The parameterization is then verified by comparing against eddy-resolving quasigeostrophic simulations, independently reproducing the results of Porta Mana and Zanna in a simpler model.


Introduction
Continuing improvement in the spatial resolution of operational global ocean models has led to recent interest in subgrid-scale parameterizations appropriate to models that partially resolve mesoscale eddy dynamics. Fox-Kemper and Menemenlis (2008) advocate a nonlinear viscosity based on Leith's (Leith, 1996) adaptation to quasi-2D dynamics of Smagorinsky's (Smagorinsky, 1963) successful large eddy simulation (LES) approach. Whereas Leith's nonlinear viscosity is based on the idea of an inertial range with a constant downscale flux of enstrophy, Jansen and Held (2014) and Jansen et al. (2015) rely on the idea of an inertial range with zero flux of energy to develop a nonlinear negative-viscosity approach similar to Sukoriansky et al. (1996). The deterministic approaches of Fox-Kemper and Menemenlis (2008); Jansen and Held (2014); Jansen et al. (2015) are complemented by stochastic approaches that model the energy transfer between resolved and unresolved scales as a random process; such stochastic models have been largely based on empirical knowledge of sub-grid eddy statistics (e.g. Berloff, 2005b;Grooms and Majda, 2013;Kitsios et al., 2013;Jansen and Held, 2014;Grooms et al., 2015).
Porta Mana and Zanna (2014) proposed a novel eddypermitting parameterization not based on LES ideas like those above. They performed thorough multiscale statistical analysis of an eddy-resolving quasigeostrophic (QG) gyre simulation, similar to the 'dynamically consistent' diagnostic framework from Berloff (2005a), studying in particular the component of the time tendency of the large-scale potential vorticity (PV) that is induced by subgrid-scale terms. Finding that existing parameterizations did not match their data well, they proposed and empirically verified an accurate parameterization of the form Eq. (5) below. This parameterization relates the subgrid-scale term to the Lagrangian time derivative of the large-scale potential vorticity. They showed that an analogy between the parameterization, which includes a time-tendency of largescale quantities, and the theory of non-Newtonian fluids of second grade or 'Rivlin-Ericksen fluids' (Rivlin and Ericksen, 1997;Dunn and Fosdick, 1974;Truesdell and Rajagopal, 2010) can be drawn. This parameterization has been successfully implemented in quasigeostrophic models, showing improvement in the mean flow, its variability and energy transfer between scales (Zanna et al., 2017). The parameterization is currently being developed for primitive-equation ocean models (Anstey and Zanna, 2017;Zanna et al., 2017). Anstey and Zanna (2017) show that some properties of QG turbulence are adequately captured by the parametrization. The aim of this paper is to obtain the parameterization of Porta Mana and Zanna (2014), including the nondimensional constant, using only assumptions of local homo-geneity and rapid spatial decorrelation of the subgridscale term.

Theoretical Development
In a quasigeostrophic (QG) model, the equation for potential vorticity (PV) evolution in a single layer is where the QG PV in the i th layer is q n , F n denotes forcing (e.g. via wind stress), and D n denotes a dissipation operator acting to remove enstrophy at small scales. The QG PV in the n th layer includes planetary, relative and stretching vorticity terms, such that it is linearly related to the streamfunction ψ n in all layers, and to the relative vorticity ∇ 2 ψ n . The advective derivative takes the form where u n = (−∂ y ψ n , ∂ x ψ n ). We will assume that a system of equations of the form (1) governs the dynamics at all scales, and approximate solutions can be computed using numerical simulations with sufficiently high resolution. The computational cost of these simulations can be prohibitive. Eddypermitting models use spatial resolution high enough to permit, but not to completely resolve typical mesoscale eddies.
The goal is to design a numerical method that accurately simulates the resolvable scales in a model with eddy-permitting resolution, with a grid-box roughly equal to the Rossby radius of deformation. To that end we begin by constructing a set of partial differential equations governing the resolvable part of the true solution. We therefore apply a time-independent spatial low-pass filter denoted · to equation (1) The dissipation term D * n in the large-scale evolution is often not equal to the low-pass filtered original dissipation term D n ; for example, the viscosity is often increased to help keep solutions smooth. The eddy source term has the form Porta Mana and Zanna (2014) ran well-resolved simulations of QG dynamics with different forcing and dissipation, applied a low-pass spatial filter to the results, and diagnosed the eddy source terms S * n directly. They compared it to a variety of parameterizations, and discovered that the data is in excellent agreement with a new parameterization of the form where ∆x is the grid size of the eddy-permitting numerical model, which is related to the length scale of the low-pass spatial filter ·. The parametrization relies on analogy between the truncated turbulent stresses and non-Newtonian stresses, arguing the need for some spatial coherence and infinitesimal memory (Rivlin and Ericksen, 1997). The goal of this section is to present a derivation of the above parameterization using the following two fundamental assumptions: • The eddy source term is rapidly-decorrelating in space.
• The eddy source term is locally-spatiallyhomogeneous.
The precise meaning of these assumptions will become clear in the course of the derivation. Without loss of generality, we consider for the remainder of this section only the top layer n = 1 and drop subscripts. First, apply the Laplacian to (3) Our assumption that S * is rapidly-decorrelating in space implies that it is dominated by small scales rather than by large-scale patterns. This assumption is supported by Fig. 5a in Porta Mana and Zanna (2014). The Laplacian of a field dominated by small scales is large, and we expect the Laplacian of the forcing term to be negligible by comparison. Dissipation specifically occurs at small scales, so it is not clear a priori that the Laplacian of the dissipation term should be smaller than the Laplacian of S * . Nevertheless, in our experiments described below this is the case, and in Porta Mana and Zanna (2014, Fig. 5d) the dissipation term is found to be much smaller than S * . Grooms et al. (2015) provide heuristic arguments and supporting evidence that S * has a Fourier spectrum growing as the 5 th power of the wavenumber, which is very strongly dominated by the small scales, evidently more so even than the dissipation term. Thus, we make the following approximation, which ignores the Laplacian of the forcing and dissipation terms in equation (6) ∇ 2 Dq Dt ≈ ∇ 2 S * .
We now consider the second-order centered finitedifference approximation to the Laplacian of S * , as used by Porta Mana and Zanna (2014) in their diagnostics. Let S * i, j be the value of S * at the location (i∆x, j∆x) on an equispaced grid in the top layer. The second-order centered finite-difference approximation to S * at the location (i∆x, j∆x) is We now seek a linear relationship between (∆x) 2 L i, j and S * i, j that will allow us to write (c∆x) 2 L i, j ≈ S * i, j , i.e. (c∆x) 2 ∇ 2 S * ≈ S * . Such a relationship immediately implies S * ≈ (c∆x) 2 ∇ 2 Dq/Dt.
The assumption of local homogeneity implies that S * i, j is a random variable with distribution approximately the same as its neighbors. The assumption of rapid spatial decorrelation implies that S * i, j is approximately uncorrelated with its neighbors. We can interpret this result by arguing that by taking the Laplacian of S * (equivalently of Dq/Dt, according to equation (7)), we are introducing information from neighbouring grid cells, that are approximately uncorrelated with that grid cell. The introduction of such information implies some random-process model for the eddy closure. We do not expect S * i, j to be completely uncorrelated with its neighbors, since that would imply a complete scale separation between resolved and unresolved scales, which is unrealistic in the eddy-permitting regime where resolved and unresolved scales are both part of an inertial range of scales; nevertheless, to simplify analysis we assume that the correlations are small enough to be negligible.
The eddy source term S * i, j and the scaled finitedifference Laplacian ∆x 2 L i, j are jointly-distributed random variables, and under these assumptions we can derive their covariance matrix, which has the form where σ 2 is the variance of S * i, j and its neighbors. The diagonal entries are the variances of S * i, j and ∆x 2 L i, j , and the off-diagonal entries are the cross-covariance.
This covariance matrix is associated with an ellipse whose axes are aligned with the eigenvectors of the covariance matrix. The eccentricity 0 ≤ ε ≤ 1 of an ellipse quantifies how flat it is: if ε = 0 the ellipse is a circle, and if ε = 1 then the ellipse is simply a line segment. Zero eccentricity would imply that there is no linear relationship between S * i, j and ∆x 2 L i, j , whereas unit eccentricity would imply that the two variables are perfectly correlated. The eccentricity is related to the ratio of the two eigenvalues of the covariance matrix, which are σ 2 (21 ± 5 √ 17)/2, and which give an eccentricity ε = 17 1/4 10 21 Thus, under the above assumptions S * i, j is very closely correlated with ∆x 2 L i, j .
The linear relationship between S * i, j and ∆x 2 L i, j is given by the eigenvector of the covariance matrix that is associated with the major axis of the ellipse (the larger eigenvalue). This eigenvector is which implies the linear relationship Making use of the approximation in equation (7) leads to The above derivation could be repeated with a finitedifference approximation to any linear differential or integral operator. But the assumption of local homogeneity requires the finite-difference stencil to be local. For example, if the Laplacian were replaced by a biharmonic operator in the above derivation the stencil would widen and the local homogeneity assumption would become less accurate when applied over a wider stencil. Similarly, using a higher-order approximation to the Laplacian would result in a wider stencil and a different linear relationship (i.e. a different constant coefficient) between S * and ∆x 2 ∇ 2 Dq/Dt, but the assumption of local homogeneity would again be less accurate over the wider stencil.

Experimental Configuration & Results
We test the above analysis in a two-layer, doublyperiodic QG model on an f -plane forced by an imposed mean shear. The nondimensional governing equations are The linear components of q and ψ are associated with an imposed uniform zonal baroclinic shear. The term multiplied by c d in equation (15) is a standard quadratic drag where the imposed zonal velocity in the lower layer −x has been subtracted from u 2 (which includes the mean flow) before computing the drag. The equations have been nondimensionalized using the deformation radius as a length scale, and the imposed zonal velocity as a velocity scale. We set c d = 0.1 and ν 4 = 0.08192. The domain is square and has nondimensional width 32 π. Approximate solutions are computed using 256×256 nonzero Fourier modes and a fourth-order semi-implicit Runge-Kutta method as described by Grooms and Majda (2014). The time step is 0.01. The grid size is 0.39, so there are just more than two grid points per deformation radius, a heuristic rule-of-thumb for eddy-resolving computations. Figure 1a shows a snapshot of the upperlayer PV q 1 from the eddy-resolving simulation. Once the simulation has reached a statistical equilibrium, 500 snapshots of the simulation with temporal spacing 0.2 are saved for diagnostic analysis.
To investigate the eddy source term S * in the upper layer, we compute The eddy-permitting biharmonic viscosity is set to ν 4 = 4ν 4 . The terms q 1 , ψ 1 , u i · ∇q i and ν 4 ∇ 4 ψ 1 are computed by applying an equal-weight average over 4 × 4 sets of grid points from the eddy-resolving simulation. The term u 1 · ∇q i is computed using the second-order energy-and enstrophy-conserving Arakawa (1966) jacobian, and the term ν 4 ∇ 4 ψ 1 is computed by iterating the standard second-order centered finite difference approximation.
The variance of ν 4 ∇ 6 ψ is 2.76 and the variance of ∇ 2 S * is 487, showing that it is not unreasonable to approximate equation (6) by equation (7). Figure 1b shows the empirical joint probability density of S * and ∇ 2 Dq 1 /Dt (which is computed using the same discrete Laplacian from section 2). The line S * = − (0.4494∆x) 2 ∇ 2 Dq 1 /Dt is plotted as a dotted line, and the linear regression S * = − (0.473∆x) 2 ∇ 2 Dq 1 /Dt is plotted as a dashed line. The best-fit regression matches the theoretical prediction reasonably well, and the elliptical contours of probability density are quite 'thin,' though not so thin as suggested by the theoreticallypredicted eccentricity of 0.995. The mismatches between theory and experiment are likely due to the neglect of spatial correlations in the eddy source term S * .

Conclusions
We have provided an a priori derivation of the parameterization proposed by Porta Mana and Zanna (2014) using only two statistical assumptions on the eddy source term S * : that its variance changes minimally from one coarse grid point to the next (local homogeneity), and that its value at one coarse grid point is approximately uncorrelated with its value at neighboring grid points (rapid spatial decorrelation). We have further verified the predictions of the theory using an eddyresolving simulation, thereby independently reproducing some of the results of Porta Mana and Zanna (2014). The derivation, being based purely on statistical arguments, provides no direct connection to the theory of non-Newtonian fluids. Other studies have however attempted to explain at a physical and mathematical level why the large scales behave in a manner analogous to non-Newtonian fluids (Foiaş et al., 2001;Porta Mana and Zanna, 2014;Anstey and Zanna, 2017).