A Breakdown in Potential Vorticity Estimation Delineates the Submesoscale‐to‐Turbulence Boundary in Large Eddy Simulations

Most submesoscale motions—fronts, eddies, filaments, and even large internal waves—are sufficiently rapidly rotating and stratified as to be strongly influenced by potential vorticity dynamics. Below the Ozmidov scale where turbulence overturns and isotropizes, potential vorticity is not commonly considered. Here, it is shown that in Large Eddy Simulations, the velocity gradients, buoyancy gradients, and potential vorticity are strongly influenced by grid‐scale processes. Grid‐scale processes in Large Eddy Simulations, as opposed to those in Direct Numerical Simulations, imply that spuriously noisy potential vorticity variance will become increasingly dominant as resolution increases—analogous to ultraviolet catastrophe. A solution, the prefiltered potential vorticity, is shown to be effective in linking the potential vorticity dynamics of the submesoscale to the nearly‐isotropic turbulent fluxes beyond the Ozmidov scale, and a derivation is provided for a set of closed conservation equations for use in interpreting potential vorticity dynamics in Large Eddy Simulations. This diagnostic approach is exceptional in that Large Eddy Simulation analysis and hydrostatic ocean modeling with parameterized turbulence analysis are harmonized.


Introduction
Potential vorticity (PV), the fundamental correlation between vorticity and stratification, is a key active variable of quasigeostrophic (QG) theory and one of the most widely used properties in large-scale oceanic and atmospheric dynamics (Pedlosky, 2013). PV is conserved via material advection for adiabatic and inviscid flows (or nearly so for atmospheric and oceanic flows), is invertible in QG flows to govern the subinertial flow, is the critical parameter governing stability in symmetric instability (Hoskins, 1974), and even when not materially conserved, the evolution of PV is tractable via frictional and diabatic fluxes (Haynes & McIntyre, 1987;Marshall & Nurser, 1992;L. N. Thomas, 2005). The central question to be addressed here is: Does 3-D turbulence, when partly resolved in a Large Eddy Simulation, count as material advection of PV or a frictional and diabatic process?
Different flow scales exhibit different properties. Large and slow flows tend to be governed by strong stratification and planetary rotation; however, as smaller, faster scales are approached-from mesoscale to submesoscale to 3-D turbulence-inertial forces come to dominate as the flow passes through isotropic turbulent properties before finally reaching the Kolmogorov (1941Kolmogorov ( , 1991 and Batchelor (1959) scales, where molecular motions are sufficiently dominant as to restrain the production of variance through (continuum) turbulence. As model capacities are increasing, multiscale simulations now have the ability to partially resolve 3-D turbulence together with larger scale flows. Since PV is strongly dependent on stratification and rotation, the question arises as to its robustness in an isotropic turbulent regime-where all mean vector quantities tend to zero because no particular directions or vectors are allowed in isotropic statistics-and mean buoyancy gradients (i.e., stratification) and mean vorticity (i.e., rotation) are zero or at best ephemeral and only partially resolved in Large Eddy Simulations (LES).
In this study, we demonstrate that PV holds different scale-dependent statistics as a 3-D turbulent regime is approached than those expected from larger mesoscale and submesoscale dynamics, and this statistical transition may effectively delineate the small-scale boundary of the submesoscale: the smallest scale of classic geophysical fluid dynamics where influences of rotation and stratification are strong. Hamlington et al. (2014) define the small-scale end of the submesoscale as the scale where nonhydrostatic effects become dominant, which is a definition based on the asymptotic behavior of reduced equations consistent with the change in PV statistics, but the latter delineation is easily detected as a statistical and dynamical marker of this scale.
Section 2 introduces the theory and behavior of PV in a 3-D turbulence regime using spectral analysis, and the implications of these results are presented in section 3 using multiscale simulations which span the submesoscale and 3-D turbulence regime. The consequences of Reynolds averaging are discussed in section 4, and this analysis is implemented in section 5 to obtain an accurate PV-along with its governing equations consistent with the needs of large-scale PV analysis-by prefiltering the velocity and buoyancy fields. Discussion and conclusions are presented in section 6.

PV and Spectral Analysis
Ertel PV, q, in nonhydrostatic, Boussinesq approximation fluids, is usefully defined by the product of the buoyancy gradient and total vorticity (Pedlosky, 2013): q ¼ ðf þ ωÞ · ∇b: (1) The relative velocity vector is u, ω = ∇ ×u is the relative vorticity, f is the planetary vorticity, and b is the buoyancy anomaly. The buoyancy is related to the potential temperature field θ by b = −gαθ under a constant salinity linear equation of state in oceanic flows or an anelastic approximation in atmospheric conditions, g being gravity acceleration and α the thermal expansion coefficient. This simplified form is that of the model studied here, although a more complete discussion of PV-including methods for treating nonlinear equations of state, background stratification gradients, and nonaligned rotation and gravity vectors-can be found in Fox-Kemper (2018).
The key utility of PV comes from two aspects: its material conservation and the invertibility of PV and buoyancy to find the balanced velocity fields (McWilliams, 2006). At the mesoscale, subinertial unbalanced motions are weak, so geostrophic turbulence (Rhines, 1979) is governed entirely by the advection of PV and buoyancy by the turbulent flows, which in turn then provide the turbulent flows from the PV and buoyancy distributions through invertibility. Large-scale PV dynamics are key to understanding the general circulation and stratification (Haynes & McIntyre, 1987;Marshall & Nurser, 1992). At the submeoscales, the Rossby number is larger, so unbalanced motions such as internal waves are strong enough to have unbalanced velocities that rival the balanced velocity field and may nonlinearly interact (Bakas & Farrell, 2009), so submesoscale turbulence is a complicated combination of wavelike and geophysical turbulence (McWilliams, 2016). Yet, many submesoscale features such as fronts, filaments, and eddies have strongly constrained motions where PV dynamics remain relevant (L. N. Thomas, 2005Thomas, , 2008; L. Thomas & Ferrari, 2008;L. N. Thomas & Lee, 2005;Wenegrat & McPhaden, 2016;Siegelman, 2020). Obviously, even as models resolve smaller and smaller phenomena in larger and larger domains, retaining a dynamical equivalent of submesoscale, mesoscale, and large-scale PV is critical to linking these theoretical analyses to interpretation of the latest model.
A second important theoretical result regarding PV is that its material conservation in inviscid, adiabatic flow extends to the full Navier-Stokes equations if diffusion and viscous transfers are correctly interpreted (Bachman et al., 2017;Haynes & McIntyre, 1987;Pedlosky, 2013). Of course, in realistic geophysical contexts, stratification and rotation effects are leading order on the length and time scale of tens of kilometers and days and are negligibly small at scales approaching the Kolmogorov (1941Kolmogorov ( , 1991 and Batchelor (1959) time and length scales. Thus, these scales are typified by order 1 Reynolds and Péclét numbers but huge Rossby and Froude numbers, thus most of the motions at these dissipative scales are unbalanced and not predictable by PV inversion.
The variance of PV can be studied by examining the power spectral density of its components across scales (following, among others, Frisch, 1995). Namely, key quantities to examine are the isotropized (i.e., integrated over all angles spanning a horizontal plane and including the polar coordinate metric factor) power spectral density of kinetic energy E, potential energy B, enstrophy Z and potential enstrophy Q, which upon integration over wavenumber magnitude are proportional to the horizontal area-integrated-and-squared velocity, buoyancy, buoyancy gradient, vorticity, PV, and tracer concentration, respectively. It should be noted that here horizontal power spectra and scales are emphasized, as most of the flows under consideration are transversely isotropic-that is, strongly anisotropic in the vertical versus horizontal directions but approximately isotropic in the horizontal over all scales, due to the vertical alignment of gravity and planetary rotation (for simplicity, the simulations are on a nontilted f-plane). Only at the smallest resolved scales does a 3-D isotropic spectrum become useful. The spectral density of a passive tracer C with concentration c will also be discussed later. For brevity, the power spectral density will just be referred to as the spectrum below, and their relation to area-integrated variance can be written out explicitly as Note that a Cartesian x-y domain and Fourier transform are assumed, again anticipating the model discretization to be analyzed. A proportionality is assumed between the integral of each power spectrum and the physical space area integral, because there are a variety of conventions as to how to choose the units of the power spectrum and thus the coefficient of proportionality may vary by convention. This proportionality holds under all conventions, however. In this manuscript, following Hamlington et al. (2014), after consistent normalization by the right-hand side of Equations 2a-2f, all comparable power spectra are rescaled by a convenient factor (often the value of the spectrum evaluated at the lowest wavenumber), which eliminates any ambiguity in definitions of the coefficients of the power spectra. This normalizing factor will be made explicit in each figure.
At the smallest scales between the Nyquist wavenumber (1/2 the sampling rate or π/Δx in angular units) and 1= ffiffi ffi 2 p times it, the grid scale acts as a filter introducing horizontal anisotropy affecting the diagnosed isotropized spectrum. However, it will be seen in the LES that the spectral effects of the subgrid scheme used in the model (Sullivan et al., 1994) extend to even larger scales. Likewise, at the largest scales, the square domain size affects the potential isotropy of the lowest wavenumbers and the wavenumber resolution of spectra in different directions.
By differentiation, the enstrophy (vorticity) spectrum is related to the kinetic energy (velocity) spectrum by Z(κ) = κ 2 E(κ) (Kraichnan, 1967), and the buoyancy gradient spectrum is related to the buoyancy spectrum by κ 2 B(κ) (Fox-Kemper et al., 2011). Note that the terminology here is somewhat inconsistent as enstrophy and energy are squared quantities, while buoyancy gradient and buoyancy are not. For consistency, one could refer to the set (2a)-(2f) as the velocity, buoyancy, buoyancy gradient, vorticity, PV, and tracer spectra, or the kinetic energy, enstrophy, buoyancy variance, and buoyancy gradient variance, potential enstrophy, and tracer concentration variance spectra. Here we will use preferentially use the former set, which is more concise on most occasions here, but in all cases, the list (2a)-(2f) indicate the correct formulation.
In the mesoscales, the interior or Charney (1971) quasigeostrophic regime where PV dynamics are best understood and are expected to govern spectra, the PV, buoyancy, and velocity spectra all scale as κ −3 , which means that the buoyancy gradient and vorticity spectra scale as κ −1 . The QG turbulence velocity spectra in this regime closely resemble those of 2-D turbulence (Kraichnan, 1967), but of course 2-D turbulence does not have stratification or PV (e.g., Bachman et al., 2017;Fox-Kemper & Menemenlis, 2008). Surface QG has a more complex spectral dependence and interpretation (Callies & Ferrari, 2013;Bodner et al., 2020).
However, in isotropic 3-D turbulence with weak temperature gradients or passive scalar gradients, the Kolmogorov (1941)-Corrsin (1951 scalings for buoyancy and velocity spectra all scale as κ −5/3 . The Ozmidov (1965) scale is generally taken as the scale at which energy density is high enough that buoyancy stratification can be considered weak in this sense. Many oceanographers prefer the Thorpe (1977) scale in identification of near-surface overturning and mixing (Dillon, 1982), but here the focus on cascade theories makes the Ozmidov scale much more natural. Thus, the buoyancy gradient and vorticity spectra are taken to scale as κ 1/3 on scales finer than the Ozmidov scale. Since PV is the covariance of the vorticity and buoyancy gradient, these scaling properties also extend to the PV field, where we expect a slope ranging between κ 1/3 (effectively assuming that buoyancy behaves as a passive scalar stirred by dominantly unbalanced motions) and κ 2/3 (when the buoyancy gradient and vorticity are perfectly correlated-e.g., motion dominated by PV inversion). These observations are not new and have been discussed in several other studies (see, e.g., Herring et al., 1994;Kunze, 2019 and references therein). However, here the emphasis is the significance of a positive slope that is terminated at the grid-scale of a LES, which implies that the smallest scales contain the largest variance per unit wavenumber and are likely to dominate the integrated vorticity, buoyancy gradient, and especially PV variance. This somewhat familiar property of PV has important implications to the study of turbulent fluxes in LES, as will be explained below. Hamlington et al. (2014) and  describe simulations that are interpreted as either extremely high-resolution submesoscale simulations or extremely large-domain nonhydrostatic LES of boundary layer turbulence, where rotation and stratification become dominant at the largest scales (see also Callies & Ferrari, 2018;Matheou et al., 2017;Pham & Sarkar, 2018;Sullivan & McWilliams, 2018;Skyllingstad & Samelson, 2012;Taylor & Ferrari, 2011;Whitt & Taylor, 2017). Hamlington et al. (2014) and  describe LES of the spindown of a temperature front in the presence of submesoscale eddies, winds, and waves, and these simulations are analyzed here. The surface wave-averaged Boussinesq equations with and without Stokes drift wave forcing are solved to simulate the interactions between submesoscale processes and boundary layer turbulence, including Langmuir turbulence. Figure 1 illustrates a snapshot of the temperature field anomaly θ − θ near the surface (z = −3 m) after the development of submesoscale instabilities, where θ ¼ 290:16 K is the background temperature field.

Modeled PV in Multiscale Simulations
Importantly for this study, these simulations span over three decades of horizontal resolution, allowing roughly a decade of classic submesoscale turbulence-exhibiting a quasi-balanced regime scaling of E(κ), B(κ) ∝ κ −2 -and a near-decade of boundary layer turbulence-exhibiting a Kolmogorov-Corrsin like scaling of E(κ), B(κ) ∝ κ −5/3 (circularly integrated spectra shown in Figure 2). In these simulations, the mixed layer is found to be 50 m in some regions and shallower in others. Here we choose to examine the spectral properties at three different depths: near surface (z = −3 m), within the mixed layer (z = −25 m), and below the mixed layer (z = −60 m).
The Ozmidov (1965) wavenumber k Oz = (N 3 /ε) 1/2 and ε are shown in Figure 3 for the no-Stokes drift simulations, where ε is the dissipation rate and N is the Brunt-Väisälä frequency. These simulations use the Sullivan et al. (1994) subgrid scheme with an explicit dependence on the estimated dissipation rate ε, thus ε can be calculated from the prognostic subgrid turbulent kinetic energy which is then used to calculate k Oz . The black circles in Figure 3 mark the depths of the spectra displayed in Figure 2. As depth increases, the Ozmidov length scale decreases, indicating that the transition to 3-D turbulence is found at smaller scales, as expected below the mixed layer, declining fast below z = −50 m. The Ozmidov wavenumber, thus estimated, proves to be a good approximation of when the u, v, w components of the velocity spectra attain a similar magnitude and behavior, a necessary criterion for isotropy. Note that at z = −60 m, the Ozmidov length scale is at roughly 8 m, closely approaching the grid scale of 4.5 m, and the interpretation of the 3-D turbulence regime, especially with regard to gradients, is likely to be inaccurate. Also note that at the near-surface z = −3 m location, near-surface suppression of the vertical velocity magnitude is apparent even below the Ozmidov scale, which is perhaps a deficiency of the LES (Ducros et al., 1998) or a sign of the similarity scaling enforced by the Sullivan et al. (1994) scheme. In any LES, isotropization is unlikely to be fully realized, and the resolution of this particular model is marginal as an LES to accomodate the large domain (Van Roekel et al., 2012); nonetheless, it is anticipated that the PV dynamics that occur in this model will be mostly generic due to the theoretical considerations from the previous section.
A critical modeling viewpoint is that LES, or more generally modeling of flows without a scale separation between parameterized unresolved phenomena and the grid scale, is meaningful only when the resolved flow variables are dominated by the larger eddies or larger, resolved scales of the power spectrum, or are at least uniformly distributed across spectral space-that is, a white power spectrum. If the smaller-scale end of the variable in question dominates, then an "ultraviolet catastrophe" occurs wherein the most important variability in that variable is contributed from the flow that occurs on scales smaller than the resolution, and any deficiencies of the subgrid scheme contaminate the basin-integrated variance of this variable at leading order (2a)-(2f).
The components of the vorticity and buoyancy gradient spectra are shown in Figure 4. These components are not yet isotropic, nor are they all oriented at a positive slope. The aforementioned limited resolution beyond the Ozmidov scale is one reason for this, as is the fact that a deliberate effort in subgrid modeling is required to have good spectral properties in these quantities near the grid scale (Johnson & Meneveau, 2018). Harcourt et al. (2002) suggest some remedies involving improved LES filtering that may also improve the spectral behavior of velocity gradients near the grid scale that were not utilized in these simulations. However, even with these deficiencies, the PV spectrum still reveals the issue at the center of this study. The spectrum of PV, that is, the product of spectral quantities in Figure 4, is shown in Figure 5. At all depths, at least one of the PV components in Figure 4 displays a positive slope nearing the smallest scales, as expected to occur when sufficiently below the Ozmidov scale. At the larger scales, the near-surface spectra displays a relatively flat slope, as expected for the submesoscale range (Capet, Klein et al., 2008). These spectral properties (flat at large scales, positive near the highest wavenumber) extend to the total PV spectra shown in Figure 5. Note that the PV is a convolution between the vorticity and buoyancy gradient when computed in wavenumber space, which highlights the correlation between these fields and may well vary across scales, so this result is not trivial.
The limited extent of the positive slope in this marginal resolution LES is theoretically predicted to extend up to the Kolmogorov microscale at least for the vorticity and buoyancy gradients independently, if not for the PV (Corrsin, 1951;Johnson & Meneveau, 2018), which is κ K ¼ ν 3 ϵ À Á 1=4 ≈ 4 mm for this simulation. At that scale, turbulent production of vorticity variance is matched by molecular viscous dissipation. By comparing this scale to the approximately 4 m horizontal resolution of this model, if the spectrum of PV continues to     Figure 4 for PV in the Stokes drift forced case (light magenta) and no-Stokes drift case (dark magenta). Thin lines denote the sharp cutoff prefiltered PV, and dotted lines denote the Gaussian prefiltered PV for both cases, as discussed in section 5. Shaded region is the same as in Figure 2. All spectra are normalized by E q (k min ) at z = −3 m of the no-Stokes case.
scale near k 1/3 , then the spectrum of PV will be 20 times higher at κ K than it is at the grid scale here. If the correlation between vorticity and buoyancy gradients persists (leading to a spectrum proportional to k 2/3 ), then the spectrum will be 460 times larger at κ K than at this model's grid scale. Figure 5 includes also the spectra of the Stokes drift forced runs where small-scale vertical transport and mixing, due to the induced Langmuir turbulence, is 4 times greater (Hamlington et al., 2014). For these runs, much of the variance is located at smaller scales and within the mixed layer and even surpasses that of the larger scale. We would expect similar results if finer resolution Stokes-free runs had been available, resolving further into the 3-D turbulence regime.
As discussed above, in interior QG simulations, the buoyancy gradient and vorticity spectra scale as κ −1 , and for near-surface submesoscale, the spectra are white, scaling as κ 0 . However, in boundary layer LES simulations, the Kolmogorov-Corrsin scalings of the buoyancy gradient and vorticity spectra scale as κ 1/3 : an ultraviolet catastrophe in both fields. LES is inconsistent if unresolved turbulence dominates over turbulence at the grid-scale, that is, when unresolved statistics dominate those of resolved large eddies that dominate many statistics of the turbulence. In this case, the PV over all resolved scales cannot be meaningfully examined, because the PV variance is dominated by unresolved scales, and even the resolved PV variance is dominated by scales that are contaminated by spurious subgrid scheme effects. Johnson and Meneveau (2018) argue that this effect can be ameliorated for the velocity gradients, for example, near grid-scale vorticity statistics can be improved through reworking LES subgrid schemes, but so far no such scheme has been developed to affect both buoyancy and momentum equations in such a way as to produce the correct near grid-scale behavior of PV.
However, LES are valuable because they perform well at predicting large-scale momentum and buoyancy variability-that is, that of the "large eddies." As PV is a derived quantity-with a dynamical equation formed by differentiation and combination of the buoyancy and momentum equations-if the buoyancy and momentum are pre-lowpass-filtered before calculating the PV field and its budget, then that system can be meaningfully related to the submesoscale and mesoscale PV budgets. Furthermore, by locating the filter scale as to average and smooth over all variability with Kolmogorov-Corrsin scalings, evaluating the prefiltered velocity and prefiltered buoyancy will exclude all of the variability in this regime and the problematic PV variability. Thus, prefiltering of the velocity and buoyancy up to the largest scale of the boundary layer turbulence is required for PV dynamics in LES to be meaningful.

Prefiltering and Reynolds Averaging
In discrete spectral space, a Reynolds average or smoothing filter can be constructed from using summation over resolved scales from the largest wavenumber based on the basin dimension to a cutoff wavenumber K (Frisch, 1995;Buzzicotti et al., 2018). In complex topography or heterogeneous flows, spectral transforms are not convenient, so a variety of other filtering approaches can be taken (with important caveats, Aluie et al., 2018;Bachman et al., 2017;Hudgins et al., 1993). Defining an average based on this cutoff, (3b) where · < denotes the discrete (horizontal) Fourier transform, and κ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k 2 þ l 2 p is the wavenumber magnitude. The essential aspect of this work is noting that because PV is the product of the vorticity and buoyancy gradient, and thus inherently nonlinear, the calculation of PV does not commute with averaging, and thus the spectra of PV are nontrivially related to the spectra of vorticity variance and buoyancy gradient variance.
By differentiation, the enstrophy spectrum is related to the energy spectrum, and the buoyancy gradient spectrum is related to the buoyancy spectrum.
This result will not be sensitive to the details of the filtering scheme in the mesoscale geophysical regime, where even the enstrophy and buoyancy gradient spectra are expected to be dominated by large scales; in the submesoscale and 3-D turbulence regimes, the choice of spectral filter will affect the result (Buzzicotti et al., 2018).
Note that the buoyancy gradient, vorticity, and velocity do not have a mean value in isotropic 3-D turbulence, as by definition an isotropic flow implies that there is no statistically meaningful directionality of any mean field: Thus, mean values of these fields will vanish, allowing LES to predict their larger scales without undue contamination by isotropic small-scale turbulence. However, due to the nonlinearilty of PV, their contribution to the PV does not vanish and may dominate the small-and even large-scale PV variance.
It is important to note that for the small-scale contributions to dominate the variance of a field, its small-scale spectrum must be both blue and deep-that is, spanning a sufficient range of scales-so as to have sufficient variance per unit wavenumber at small scales to dominate the overall contribution to integrals such as (2c)-(2e). The present theory of rotating, stratified turbulence in seawater anticipates no meaningful deviations from the cascade behavior at scales slightly smaller than the Ozmidov scale where isotropization occurs (in flows where N > f) and those near the Kolmogorov scale where viscosity becomes dominant (for fluids where viscosity is larger than diffusivity), so if the spectrum is blue just beyond the Ozmidov scale, then it is expected to remain so until the molecular viscosity takes over near the Kolmogorov scale (Kunze, 2019).
Note that to be consistent with Reynolds averaging, here the simulations are processed using a sharp cutoff with no dealiasing (Orszag, 1971) rather than a smoother roll-off, to preserve (i) the idempotent property that defines a Reynolds average, · ¼ · and · ′ ¼ 0, where overbar indicates an average and primes denote the perturbation from the average, and (ii) the property that temporal and spatial derivatives commute with the cutoff filter. An example of a smooth Gaussian filter will be presented in the following section for comparison.
However, to retain the generality of the equations, we will avoid assuming that the averaging or filtering operation is a Reynolds average or that the spatial derivatives commute with the filter. This involves using the total turbulent stress divergence form ∇ · uu − ∇ · u u rather than the Reynolds stress divergence ∇ · R ≡ ∇ · u ′ u′. If the filter commutes with differentiation but is not idempotent, then the total turbulent stress tensor that results includes the Leonard stress (L ≡ u u − u u) of resolved-but-turbulent motions, the cross stress (C ≡ u′ u þ uu′), and the Reynolds stress contributions (Leonard, 1975). In practice, only the Leonard stress will be explicit, as much of the cross and Reynolds stresses will be parameterized in LES. The same approach is taken with the buoyancy fluxes.
The filter will be assumed to commute with time derivatives, but the commutation of the filter with spatial derivatives will not be assumed, which is difficult to achieve in practice-a sharp spectral cutoff does so, as does an ensemble average. Instead, only one additional constraint on the filter will be required, which is where P is the pressure. It will also be assumed that which makes the prefiltered PV materially conserved as well as conserved in flux form. Note that (6) and (7) are much less restrictive than assuming that all derivatives commute with the filter.
With these definitions in mind, the consideration of the prefiltered PV equation begins with the buoyancy and momentum equations.

Prefiltered PV
The Boussinesq governing equations for the momentum and buoyancy are F and D are the frictional and diabatic terms, respectively, which in an LES-like model are the parameterized unresolved turbulent fluxes (e.g., Deardorff, 1980;Moeng, 1984) and in DNS are molecular viscosity and diffusivity. f is the Coriolis parameter vector, which is not required to be vertical or constant in space, although it will be both in the simulations analyzed and assumed to be slowly varying enough that it is unaffected by averaging. Consistent with the Boussinesq approximation, it is also assumed that ∇ · u = 0.
Forming the vorticity equation and the gradient of the buoyancy equation from (8a)-(8b), the PV equation is formed first by multiplying these times the appropriate factors.
(9b) and summing to form the product rule, which cancels out the vortex stretching term in square brackets, leaving behind the forcing terms in curly brackets, It will also be assumed that ∇ × k = 0, which is true in flat domains where k (the direction of gravity) is spatially constant and on a sphere where k is radially oriented, to make ∇b · ∇ × ðbkÞ ¼ ∇b · b∇ × k ð Þ¼0.
To summarize, the resulting Boussinesq PV equation is which can also be written as where J ¼ uq − F × ∇b − ðf þ ωÞD are the "J-fluxes" of PV. In both forms of the PV equation, the PV evolves according to both the frictional and diabatic effects, F; D, and the advective fluxes, u q, but only in (12) is it apparent that all of these effects can be written in flux form (Haynes & McIntyre, 1987;Marshall & Nurser, 1992).
However, in the 3-D turbulence regime, the curl of the momentum equation and the gradient of the buoyancy equation are dominated by the gridscale noise, so the resulting PV equation is contaminated. By contrast, the prefiltered PV is constructed from the prefiltered velocity and buoyancy fields, to avoid the ultraviolet catastrophes.q Note thatq ≠ q, as the velocity and buoyancy are filtered before formingq, whereas q is filtered after the PV is formed from the unfiltered fields. As PV is a nonlinear combination of buoyancy and momentum, these operations do not commute. Also, as PV is nonlinear, care to avoid aliasing unresolved buoyancy gradient or vorticity wavenumbers into the resolved range when calculating the prefiltered PV is needed, dealiasing of quadratic nonlinearities is guaranteed if the filter eliminates one third or more of the resolved wavenumbers (Orszag, 1971). However, the prefiltering advocated from a physical perspective goes far beyond this numerical constraint: In the simulations analyzed, the optimal filter wavenumber is near the 15th out of 2,049. As with kinetic energy of filtered fields , an exact budget for the prefiltered PV can be constructed and analyzed. Figure 6 shows a plane view of the near-surface PV before and after prefiltering of the velocity and buoyancy fields using a sharp cutoff filter and a standard Gaussian filter with the same cutoff scale at 400 m. Highlighted by the shaded region in Figure 5, the 400 m cutoff scale is chosen here to correspond with the loss of hydrostatic balance found in Hamlington et al. (2014), which also corresponds to the scale of change in slope of gradients ( Figure 4) and of PV near the surface and within the mixed layer ( Figure 5). The unfiltered PV is contaminated by small-scale variance that distracts from the larger scale features which appear after the prefiltering is applied (thick vs. thin and dotted lines in Figure 5). The advantages of choosing a sharp cutoff are discussed in the previous section, yet one disadvantage is the ringing effect as evident in Figure 6, which is otherwise smoothed by the Gaussian filter. In the Stokes drift forced case, where small-scale mixing is very intense near the surface, the PV is contaminated heavily at all scales. The prefiltered PV in Figure 5 exhibits a flat slope throughout most of the submesoscale. This particular property of Stokes forced mixing will be examined more carefully in a future study. Furthermore, the transports of the unfiltered PV cannot be obtained cleanly from the model, as will be shown next.
The filtered versions of the momentum and buoyancy Equations 8a and 8b are where the superscript plus symbol indicates the frictional and diabatic effects including the total turbulent Figure 6. Plane view of a snapshot of the near-surface (z = −3 m) PV field without prefiltering of the velocity and buoyancy fields (left) and after using a sharp cutoff filter (center) and Gaussian filter (right) on the buoyancy and velocity fields separately before calculating the PV as in (13). The cutoff scale here is 400 m, which corresponds to the shaded regions in Figures 2, 4, and 5.
stress and total turbulent buoyancy flux and any potential non-commutation of filter and spatial derivatives. It is assumed that temporal and spatial derivatives commute with the filtering operation, which also guarantees that ∇ · u ¼ 0. Following the same procedure as that to form (11), but using (14a)-(14b) in place of (8a)-(8b), the prefiltered PV equivalent of (11) is where D Dt includes advection by ū only . Equation 15 can also be written in flux form, Figure 7. Same as in Figure 5 for the J-flux spectra in the no-Stokes drift case. Shown is the total J-flux (top row) and its components: qu (second row), F × ∇b (third row), and ðf þ ωÞD (bottom row). Thick lines denote the unfiltered spectra, thin lines the sharp cutoff prefiltering, and dotted lines the Gaussian prefiltering as described in section 5. All spectra are normalized by E x · J (k min ) at z = −3 m.
Note that no commutation of spatial derivatives and averages is assumed in this derivation (only ∇ · u ¼ 0), and the prefiltered PV Equation 15 has the expected form of a PV conservation equation, with the "frictional" and "diffusive" operators (that actually also contain turbulent transports): As (15) and (16) have the expected form, the prefiltered PV behaves as expected for all geophysical applications.
These are the operators that need to be assessed in LES to understand PV injection at the surface, not just the molecular or subgrid friction and diffusion. It is worthy to note that if the filtering operator is standard Reynolds averaging, Equations 17a and 17b arrive at a similar expression as in Keyser and Rotunno (1990), following the original work by Staley (1960) and Shapiro (1976Shapiro ( , 1978. However, our emphasis here is on the particular application in the Kolmogorov regime using LES, which is not addressed by the previous authors. The J-fluxes as written in (16) remain meaningful if the dissipative operators in (17a) and (17b) are used, and even though F, D, ∇ · (uu), and ∇ · (u b) may be noisy, they are of exactly the form solved for by the modeled momentum and buoyancy equations at full resolution, including any oddities of discretization, dealiasing, etc. Under the prefiltered PV approach, these noisy terms are never converted into J-fluxes before filtering, which leads to large aliasing contaminating even the large-scale J-fluxes (Figure 7). To remove aliasing that will arise in calculating the advective J-fluxes (including the implied cubic nonlinearity in PV and energy budgets), the filter wavenumber should be half of the cutoff wavenumber (which is π over the gridscale) or below, equivalent to a smoothing operator over 4Δx or greater. This rule is stricter than filtering at the wavenumber that is two thirds of the cutoff wavenumber allowed under quadratic nonlinearity, equivalent to smoothing over 3Δx or greater for quadratic nonlinearity (Orszag, 1971). By contrast, if one tries to calculate the advective J-fluxes that result from simply averaging (11), that is, ∇ · uq and ∇ · u ½∇ × u þ f · ∇b ð Þ , the calculation involves differentiating a triple correlation that the model never predicts directly, is noisy, and depends on scales outside of those resolved in an LES-that is, aliasing. To a degree this cannot managed by normal dealiasing filters because it is a cubic instead of a quadratic nonlinearity, requiring filtering of not just one third of resolved wavenumbers, but one half (Orszag, 1971). Figure 7 shows that even the large-scale J-fluxes are contaminated by aliasing between the small-scale and unresolved modes. Likewise, the frictional J-fluxes involve products of buoyancy gradients, vorticity, with noisy frictional terms which are not consistent with the discrete model prognostic equations.

Discussion and Conclusions
Thus, to meaningfully study PV in LES, via material conservation or via J-fluxes, a prefiltered PV as in (13) with a physically selected filter scale ≥ 4Δx is the most meaningful quantity, as it avoids spurious effects, aliasing, and noise. Momentum and buoyancy fluxes on scales below the submesoscale, including those of 3-D turbulence, should be treated alongside parameterizations of momentum and buoyancy fluxes and bundled into the frictional and diabatic terms: F þ ; D þ as written in (17a) and (17b). By this method, correct interpretations of the PV evolution equations, (11)-(12), Kelvin's circulation theorem, the impermeability theorem (Haynes & McIntyre, 1987;Marshall & Nurser, 1992), and PV invertibility on the scales where it is the dominant constraint will be recovered.
In interior QG simulations, the PV, buoyancy, and velocity spectra all scale as κ −3 , which means that the buoyancy gradient and vorticity spectra still scale as κ −1 , and so are well behaved even without filtering. Note that if the 3-D turbulence scale is not resolved (e.g., in a hydrostatic general circulation model), the turbulent parameterizations F and D already contain the correct terms without additional prefiltering. Thus, the distinction between the full PV and prefiltered PV is not important in interior QG.
In near-surface submesoscale simulations, buoyancy and (horizontal) velocity spectra all scale as κ −2 . Thus, the buoyancy gradient and velocity gradient spectra are white, scaling as κ 0 . However, the vertical vorticity spectrum is likely to resemble the interior QG scaling even as the surface velocity divergence spectrum does not (Callies & Ferrari, 2013;Chereskin et al., 2019;Rocha et al., 2016), which means that PV may be slightly red and thus meaningfully calculated without prefiltering, although prefiltering may still be desired to reduce aliasing in calculating advective terms. The parameterizations of boundary layer turbulence in such models are included in the dissipative terms.
In Direct Numerical Simulations, the molecular viscous and Kolmogorov and Batchelor scales are thoroughly resolved (with at least four gridpoints). In this case, prefiltering across the dissipative scales will not affect the result (e.g., Herring et al., 1994), and cubic nonlinearity de-aliasing will suffice.
In boundary layer LES simulations, the Kolmogorov-Corrsin scalings for buoyancy and velocity spectra all scale as κ −5/3 , and the turbulence is isotropic. Thus, the buoyancy gradient and vorticity spectra scale as κ 1/3 : an ultraviolet catastrophe in both fields. Furthermore, LES is inconsistent with the assumption that no important turbulence effects occur at scales near to or finer than the grid scale. In this case, the full PV over all resolved scales cannot be meaningfully examined, because its variance is dominated by unresolved scales and its advective J-fluxes are dominated by triple correlations that depend on unresolved scales (Clark et al., 1979;Leonard, 1975). Even though backscatter parameterizations may improve velocity gradient statistics (Johnson & Meneveau, 2018), they have not been designed to correctly inject PV variance (if this is indeed possible). However, averaging or filtering the PV affects interpretation (Keyser & Rotunno, 1990). Thus, for the case of the Hamlington et al. (2014) and  and related LES-includingsubmesoscale simulations, a course for how to manage the interpretation of PV analysis covering the transition from the submesoscale to the boundary layer scalings emerges. Prefiltering of velocity and buoyancy up to the O(H) scale, that is, the largest scale of the boundary layer turbulence, is required before calculation of the PV, and here a meaningful construction of the resulting PV dynamics is given in (15)-(16). Although alternative expressions of PV may be useful for diagnostic purposes (e.g., Morel et al., 2019), prefiltering is still necessary to avoid contamination from the smallest scales in either diagnostic method. In other circumstances, it may be more convenient to estimate this scale by the Ozmidov scale, by where the unfiltered PV spectrum inflection point lies, by the scale where nonhydrostatic effects (i.e., acceleration terms in the vertical momentum equation) become pronouced, by where the aspect ratio of motions in 3-D is order 1, or by where the u,v,w spectra all converge and then proceed on smaller scales with a k −5/3 slope. Particular classes of less well-understood turbulence, such as that emanating alongside coherent structures (Taylor & Ferrari, 2009) or forced toward anisotropy by convection (Rincon, 2006), stratification (Kunze, 2019), or Stokes forces as in Langmuir turbulence (Hamlington et al., 2014;, merit a variety of these checks be made during the selection of a filter scale or procedure. Indeed, in Hamlington et al. (2014), the agreement among estimates of the isotropization scale differs between the simulation including Stokes forces and one without it (e.g., Figure 5).
While the Boussinesq equations with buoyancy have been used here in the tradition of oceanic boundary layer LES (e.g., McWilliams et al., 1997), there is no reason why atmospheric boundary layer LES (e.g., Matheou et al., 2017) or cloud-resolving LES (e.g., Yamaguchi et al., 2017) would not also suffer the same issues with identifying the relevant PV dynamics or energy budgets. Other issues have been identified with PV intended to be used for atmospheric turbulence simulations (Herring et al., 1994), but here the focus is the method for producing the prefiltered PV, why it is the right choice, and its evolution equations.