A Kinematic Kinetic Energy Backscatter Parametrization: From Implementation to Global Ocean Simulations

Ocean models at eddy‐permitting resolution are generally overdissipative, damping the intensity of the mesoscale eddy field. To reduce overdissipation, we propose a simplified, kinematic energy backscatter parametrization built into the viscosity operator in conjunction with a new flow‐dependent coefficient of viscosity based on nearest neighbor velocity differences. The new scheme mitigates excessive dissipation of energy and improves global ocean simulations at eddy‐permitting resolution. We find that kinematic backscatter substantially raises simulated eddy kinetic energy, similar to an alternative, previously proposed dynamic backscatter parametrization. While dynamic backscatter is scale‐aware and energetically more consistent, its implementation is more complex. Furthermore, it turns out to be computationally more expensive, as it applies, among other things, an additional prognostic subgrid energy equation. The kinematic backscatter proposed here, by contrast, comes at no additional computational cost, following the principle of simplicity. Our primary focus is the discretization on triangular unstructured meshes with cell placement of velocities (an analog of B‐grids), as employed by the Finite‐volumE Sea ice‐Ocean Model (FESOM2). The kinematic backscatter scheme with the new viscosity coefficient is implemented in FESOM2 and tested in the simplified geometry of a zonally reentrant channel as well as in a global ocean simulation on a 1/4° mesh. This first version of the new kinematic backscatter needs to be tuned to the specific resolution regime of the simulation. However, the tuning relies on a single parameter, emphasizing the overall practicality of the approach.


Introduction
Horizontal viscosity used in numerical ocean models is thought to roughly represent the effect of unresolved scales on the resolved flows. In most cases, harmonic or biharmonic viscosity operators are used with viscosity coefficients either scaled with an appropriate power of mesh cell size or relying on the Smagorinsky (1963) or Leith (1996) parametrizations (see, e.g., Bachman et al., 2017;Fox-Kemper & Menemenlis, 2008). The latter are flow-aware parametrizations in which the amplitude of the viscosity coefficient depends on the strain rate or on vorticity gradients, usually derived under the assumption of fully developed turbulence with energy or enstrophy cascades to small scales (Fox-Kemper & Menemenlis, 2008). However, the premises of 3-D homogeneous isotropic turbulence (Smagorinsky) or homogeneous isotropic quasi-geostrophic turbulence obeying the −3 spectral law (Leith) are not observed in a spatially and temporally local sense in realizations of eddying ocean flows. Therefore, strictly speaking, the reasoning based on the turbulent scaling laws only provides a heuristic argument. Originally, the Smagorinsky and Leith coefficients of viscosity were derived for harmonic operators but can be generalized to biharmonic operators as, for example, proposed by Griffies and Hallberg (2000) for the biharmonic Smagorinsky coefficient of viscosity. The Leith coefficient of viscosity can be generalized to include the dependence on the gradients of horizontal divergence (Fox-Kemper & Menemenlis, 2008) or gradients of quasi-geostrophic potential vorticity (QG Leith; see Bachman et al., 2017). Bachman et al. (2017) and Pearson et al. (2017) show that use of the harmonic Leith coefficient of viscosity may lead to a reduced overall dissipation in high-resolution simulations compared to a constant-coefficient biharmonic viscosity.
In turbulent flows, energy is transferred between large and small scales. Numerical modeling of turbulent flows relies on discretizations of the equations of motion and hence introduces a finite cutoff scale beyond which motion is not resolved. The term physical energy backscatter accounts for a missing energy transfer from the unresolved scales to the resolved flow (e.g., Frederiksen & Davies, 1997). Various forms of diagnostics have been applied to develop both deterministic and stochastic parametrizations to simulate the backscatter from unresolved (mesoscale) eddies to the resolved flow in idealized ocean models (e.g., Berloff, 2018;Grooms et al., 2015;Kitsios et al., 2013;Mana & Zanna, 2014). Besides physical backscatter, the idea of compensating excessive energy dissipation of viscous closures through a numerically motivated backscatter parametrization has first appeared in atmospheric modeling (e.g., Berner et al., 2009;Shutts, 2005). It reduces overdissipation through an additional, antidissipative operator while maintaining numerical model stability.
In this study, we address the numerical kind of backscatter, henceforth referred to as an "(energy) backscatter parametrization". It is particularly relevant for ocean models at eddy-permitting resolutions because the scales at which viscous dissipation is noticeable may occupy the entire range of resolved scales, as illustrated, for example, in Soufflet et al. (2016). This affects the release of available eddy potential energy, further weakening the eddy kinetic energy (EKE). For this reason, the concept of energy backscatter parametrizations was extended to ocean models as a means to reduce the accompanying EKE dissipation (Bachman, 2019;Jansen & Held, 2014;Jansen et al., 2015Jansen et al., , 2019Juricke et al., 2019Juricke et al., , 2020Klöwer et al., 2018). According to this concept, the dissipation of EKE is diagnosed, and partly compensated for at larger scales, either deterministically through negative viscosity or stochastically through random forcing of the momentum equation. Since the dissipated energy must be returned at scales essentially larger than the scales of dissipation, a common approach is to combine biharmonic viscosity and negative harmonic viscosity. Juricke et al. (2019) proposed to use a spatially smoothed negative viscosity term to increase the scales at which energy is scattered back to the resolved flow. Here we develop and implement a simplified version of this concept, henceforth referred to as kinematic backscatter. Different from the parametrization of Juricke et al. (2019)-henceforth called dynamic backscatter-kinematic backscatter does not involve additional prognostic variables in its simplest form. It only attempts to reduce dissipation and ignores any feedback that may control the amplitude of the reduction in dissipation. It also acts instantaneously and therefore does not include the memory effect of the dynamic backscatter. Whether this memory, and if so, in what form and with what time constants, might give even better results-particularly for inhomogeneous flows-is not easy to assess in detail and will be subject to future work.
Energy backscatter parametrizations aim at returning the excessively dissipated energy back to the flow. However, backscatter parametrizations are paired with a viscosity closure and the amount of excessive dissipation depends on the chosen closure. A more spatially selective or localized estimate of viscosity coefficients limits dissipation to very specific areas and may lead to a reduction of overall dissipation. The Smagorinsky and Leith parametrizations are examples of such spatially localized, flow-aware parametrizations. We propose a simple estimate of the viscosity coefficient that relies on the velocity difference across the face (i.e., edge if viewed from the surface) of a velocity control volume and demonstrate that, in certain cases, it can further reduce dissipation on its own when compared to the above-mentioned viscosity schemes.
Part of our interest in numerical viscosity operators-and energy backscatter parametrizations in particular -comes from the fact that some unstructured-mesh discretizations used presently for large-scale ocean modeling maintain spurious branches of wave dispersion related to the geometry of placing discrete degrees of freedom on triangular meshes (e.g., Danilov et al., 2017;Korn, 2017). The spurious branches may lead to enhanced grid-scale variability. Controlling it may require a somewhat stronger viscous dissipation than commonly used in codes designed for regular quadrilateral meshes, which exacerbates the problem of overdissipation if mesh resolution is only eddy-permitting. We sought an implementation of a viscosity operator as part of a new energy backscatter parametrization that would allow us to control spurious grid-scale oscillations without excessively damping the resolved dynamics.
In this paper, we explore the simplest version of kinematic backscatter in conjunction with the new viscosity coefficient. We mention extensions to this concept but leave the analysis of other options for the future. Furthermore, we specifically target the Finite-volumE Sea ice-Ocean Model (FESOM2) which utilizes unstructured triangular meshes and a cell placement of velocities. Although our detailed considerations are specific to this particular geometry, the concept of the parametrization is general enough to be applicable elsewhere after necessary adjustments.
The structure of the paper is as follows. We begin with explaining the general idea of kinematic backscatter on a 1-D example in section 2. We introduce the new viscosity coefficient in section 3. The implementation of kinematic backscatter on planar triangular grids-together with the implementation of the appropriate viscosity operator and the new viscosity coefficient-is presented in section 4. In section 5, we compare the performance of kinematic backscatter to some classical viscosity closures and to the dynamic backscatter of Juricke et al. (2019) in the channel setting of Soufflet et al. (2016). We consider an eddy-permitting mesh at 20-km resolution (section 5.1) as well as a barely eddy-resolving mesh at 10-km resolution (section 5.2). Results for a global, eddy-permitting ocean simulation are presented in section 6. Section 7 closes with a discussion and conclusions.

The 1-D Consideration
To illustrate the idea of kinematic backscatter, we will explain the concept for a uniform one-dimensional mesh. In this setting, the ideas can be formulated in spectral language. The key ideas carry over to the two-dimensional setting, but we caution the reader that, particularly with nonconstant coefficients of viscosity and nonquadrilateral or nonuniform grids, spectral language may be neither directly applicable nor an adequate description.
On a uniform one-dimensional mesh, velocity points u n are located at cell centers x n = nh, where n is the index of the velocity point and h is the distance between these points. We begin with the case of a constant viscosity coefficient ν. Then the simplest viscous operator is the finite difference Laplacian The operator in (1) that is, the power of viscous forces is negative if summed over all mesh cells. This is more readily seen in the Fourier representation, which will be used in this section. Assuming u n ¼ u a e ikxn , where u a is the amplitude and k the wavenumber, (1) takes the form where is the spectral symbol of V. The dissipative character of V corresponds to the nonpositivity of the spectral symbol V k . At the largest resolvable wavenumber jkj ¼ π=h, we have V k ¼ −4ν=h 2 ; when jkhj is small, As a second key operator, we introduce the filter 10.1029/2020MS002175

Journal of Advances in Modeling Earth Systems
It is obtained by first defining values at the cell boundaries x n + 1/2 by averaging the velocity of cells n and n + 1 then averaging back from cell boundaries to cell centers. Its spectral symbol is The key idea of kinematic backscatter is to modify the viscous operator such that it is antidissipative on some range of scales. The most simple such modification is given by where α is a numerical factor. The spectral symbol of V ð1Þ is V Since the spectral symbol of the smoothing operator F k ¼ 0 for jkj ¼ π=h, this selection will not affect dissipation at the grid scale. It will, however, reduce dissipation at all resolved wavenumbers because F k ≥ 0 and may even lead to an overall antidissipation (positivity of V ð1Þ k ) in some range of k if α > 1. Kinematic backscatter, therefore, consists-in an area-averaged sense-of a dissipative component, V, which is the viscosity operator, and an antidissipative backscatter component, −αFV, whose strength depends on the coefficient α.
For F defined by (5), F − 1 is, up to the factor 4/h 2 , the discrete Laplacian and Thus, in the particular case when α ¼ 1, (7) becomes V ð1Þ ¼ −ðh 2 =4νÞ V 2 , that is, the two components of (7) reduce to a biharmonic operator with the biharmonic coefficient of viscosity νh 2 /4. When α > 1, we enter the actual backscatter regime because the combined operator (7) is antidissipative for some range of k. The stability of the resulting scheme depends on the shape of the kinetic energy spectrum, hence, on the efficiency of nonlinear interactions to distribute the backscattered energy across scales. Stability should be warranted if α is sufficiently close to one, as discussed in Appendix A1. The precise range of stability is difficult to ascertain theoretically; finding it requires experimentation in practice.
There is a connection between the constant coefficient α and the variable coefficient R dis used in the dynamic backscatter of Juricke et al. (2019). The coefficient R dis controls how much energy enters the subgrid energy budget, which in turn defines how much energy is available for reinjection. A spatially local backscatter coefficient is then scaled by this subgrid energy budget. The kinematic backscatter coefficient α also scales the backscatter component of (7), but in an instantaneous sense, where it affects how much of the dissipated energy at each location is scattered back straightaway.
One can do further steps and introduce We remark that in the special case 1 − α þ β ¼ 0, the right-hand side factorizes as follows: In other words, due to (8), the choice (10) implements the simplest possible kinematic backscatter on top of a biharmonic viscosity. In particular, when β ¼ 1 in (10), the operator reduces to a triharmonic viscosity so that the backscatter regime is characterized by β > 1.
Since the coefficient of viscosity will be chosen as a function of the velocity field in the next section, (11) is a nonlinear operator; in particular, the simple relation (8) ceases to hold and the resulting operators V ð1Þ or V ð2Þ with, respectively, α ¼ 1 or β ¼ 1 will differ from their biharmonic or triharmonic counterparts. The advantage of using V ð1Þ or V ð2Þ as substitutes for biharmonic or triharmonic operators is the absence of higher-order derivatives in computations. Its drawback is the lack of (spatially) local momentum conservation for it ceases to be a divergence of a flux. This is not expected to create a problem: The total, domain-averaged momentum is conserved, and we only redistribute viscous sinks or sources without artificially creating them.
Most of the 1-D discussion carries over to discretizations on general meshes. In particular, α ¼ 1 still marks the boundary between purely dissipative behavior and antidissipative behavior on some range of wavenumbers. This can be shown as follows. Any consistent averaging operator F must preserve constants and dissipate any spatially varying function. Thus, ð1 − α FÞ is a positive operator so long as α ≤ 1; for any α > 1, ð1 − αFÞ acting on a constant field is negative and thus will remain negative if the field varies on sufficiently large scales. As α is further increased, a successively larger range of scales will be affected.
Further, in the limit of small wavenumbers, we can expand the spectral symbol of F about wavenumber zero. It can be shown that in this limit the spectral symbol of ð1 − FÞ corresponds to the Laplacian if the grid and F are sufficiently isotropic. Thus, for α ¼ 1 and a constant coefficient of viscosity, V (1) still corresponds to a biharmonic operator. The details of behavior at large wavenumbers will be different but do not change the main idea. We will therefore use the concepts and terminology developed here as motivation for the implementation of kinematic backscatter on triangular grids in the context of horizontal discretizations of general ocean circulation models.

A New Simple Form for the Coefficient of Viscosity
Kinematic backscatter is always tied to a suitable form of the coefficient of viscosity. In the following, we describe a new, simple choice for a flow-dependent coefficient that is cheap to compute, has good spatial selectivity, and helps the suppression of grid modes on triangular meshes. While designed in the context of the triangular quasi-B-grid, it is not confined to this type of mesh.
To begin with, any of the traditional flow-aware and scale-selective parametrizations for the viscosity coefficient can be chosen to replace a constant viscosity coefficient ν. However, these generally require computations which are defined on a stencil much wider than just nearest neighbors. To reduce dissipation, it is sensible to use flow-dependent viscosity coefficients that rely on the smallest possible stencil so that the overall area where viscosity is high will be kept at minimum. Furthermore, for FESOM2 the coefficient needs to be able to efficiently control velocity differences at grid scales. For this reason, we propose the simple form (given here for a two-dimensional mesh) where c and c′ are the grid cells on which the discretized velocities are placed, ℓ c′c is the length of the edge between the cells, γ 0 is a small background velocity amplitude of about 0.001 m/s, needed in regions of low flow speed or weak shear, and γ 1 is a small dimensionless coefficient. Both γ 0 and γ 1 are subject to tuning. Such tuning will be discussed for γ 1 in section 5. We find that γ 1 ≈ 0.03-0.1 works well. As it is based on velocity differences, (12) is reminiscent of the Smagorinsky viscosity and the numerical coefficient γ 1 can be compared to C/π 2 of the Smagorinsky parametrization (e.g., Fox-Kemper & Menemenlis, 2008). As will become clear in section 5, the advantage of this viscosity coefficient is that its calculation is based on the same variables as the viscosity operator, so that the computation of the coefficient as part of an operator Vu does not cause additional computational costs. Although (12) is not invariant to solid body rotations, this is the drawback of the whole operator on triangular meshes as we will see in section 4.
Up to the small background part, (12) is a flow-aware parametrization that strongly weights the differences in across-edge velocities. This design is motivated by the numerical consideration that the viscous closure should be able to counteract grid modes-see the discussion in Juricke et al. (2019).

Triangular Meshes
We now extend the idea of kinematic backscatter to the triangular quasi-B-grid employed by FESOM2. Some details of the viscosity operator are already described, for example, in Danilov et al. (2017) and Juricke et al. (2019). Here, we use this operator together with the new coefficient of viscosity described in section 3, so that we present a complete, self-contained description of the scheme.

Discrete Harmonic Operators on Triangular Meshes
Physically, a viscous force V is the divergence of the tensor of viscous Here, we only discuss viscosity in the horizontal, so that ∇ denotes the horizontal gradient operator. In the case of classic har- If ν is variable, (13) still defines a dissipative operator, but it is not invariant under solid-body rotations. For numerical reasons, it is nonetheless convenient to use a discrete form of (13) to define the viscous operator on unstructured meshes. For example, on triangular or hexagonal unstructured C-grids, the full stress tensor σ cannot be computed in terms of quantities naturally defined on the C-grid (Gassmann, 2018), while which is equivalent to (13) for constant ν, involves only the computation of curl and divergence, which are naturally defined quantities.
For FESOM2, which uses a quasi-B-grid on horizontal triangles, it would be possible to compute σ. However, its stencil would extend beyond the nearest neighbor cells and fail to include the contributions from the nearest velocities with sufficient weights to control grid modes. In contrast, the operator (13) is easily computed in a form emphasizing the nearest neighbor contributions . Higher-order viscous operators can be constructed by iterating discrete forms of (13) or (14).
Based on these considerations, noting that in FESOM2, the discrete velocities are placed at cell centers, we define an approximate harmonic viscous operator V via where NðcÞ is the set of all neighbor cells of cell c (i.e., the cells sharing a common edge), S c is the area of triangle c, ℓ c′c is now the length of the edge between triangular cells c′ and c, and r c′c is the distance between the centroid of cell c and the centroid of cell c′( Figure 1).
When the mesh cells are equilateral triangles, ðu c′ − u c Þ=r c′c is the velocity gradient projection on the outward normal (from c) at the edge between c and c′, so that the right-hand side of (15) is the sum of viscous fluxes leaving cell c, thus giving a consistent discretization of the harmonic viscosity operator (13). Although this interpretation fails on general unstructured meshes, (15) conserves momentum (if viscous stresses are zero on boundaries) and dissipates kinetic energy. It is therefore a numerical filter of grid-scale motion, acting directly upon differences between the nearest neighbor velocities. While this interpretation may seem undesirable to a physicist, it causes, in practice, no visible numerical instabilities such as grid-scale patterns. Moreover, the use of harmonic or biharmonic viscous operators to maintain a turbulent energy cascade is not based on first-principle physics in the first place and neither is the selection of the viscous coefficients (see, e.g., Fox-Kemper & Menemenlis, 2008). The approach advocated here is, therefore, at the same level of modeling as any other viscous closure, motivated by the necessity to effectively control the smoothness of cell velocities and subject to empirical evaluation.

Journal of Advances in Modeling Earth Systems
Since the geometric quantities in (15) and the mean viscosity are defined on edges, they can be incorporated into a generalized edge viscosity ν c′c , which is symmetric between c and c′, so that (15) takes the form Summing over all cells c (with appropriate boundary conditions), ∑ c S c ðVuÞ c ¼ 0 because the difference between u c′ and u c appears with opposite signs in the expressions for cell c′ and c. This means that momentum is globally conserved. Furthermore, taking the total viscous dissipation ∑ c S c u c · ðVuÞ c , we observe that the contribution from c′ and c appears twice, one time as ν cc′ u c · ðu c′ − u c Þ and the other time as ν c′c u c′ · ðu c − u c′ Þ, summing to −ν c′c ju c′ − u c j 2 . This means that kinetic energy is dissipated. The free-slip boundary conditions in this approach imply that boundary edges are not included in the computations. For no-slip boundary conditions, ghost velocities u c′ across boundary edges are used, directed opposite to the interior velocities u c .
The viscosity coefficient ν cc′ of (12) introduced in section 3 is readily applied with the viscous operator (16). The advantage of (12) becomes apparent now, since it is based on the same data as the right-hand side of (16) so that the coefficient (12) can be computed as part of the operator Vu at no extra cost.
In our experiments, see section 5 below, we find that the viscous operator (16) with coefficient (12) is less dissipative than other common viscosity closures. Yet it is still overdissipative. Thus, we now turn to the implementation of kinematic backscatter on the triangular mesh to further reduce total dissipation.

Averaging and Kinematic Backscatter
In our implementation, we use the same filter as in Juricke et al. (2019). It is realized as the composition of two averaging operations, denoted X and C. First, area-weighted cell values are averaged to common vertices. Second, the vertex values are averaged back to cell centroids. If a c is a quantity defined at cells, its vertex values are obtained as where CðνÞ denotes the set of cells containing vertex v. The area S v ¼ ∑ c ∈ CðvÞ ðS c =3Þ in the denominator is that of the median-dual control volume around vertex v. Each cell contributes with weight S c /3 to each of its vertices. We note that where the sums are over all cells and vertices, respectively. This shows that X preserves area integrals of a. Similarly, for a quantity b v defined on vertices, its averaged cell values are computed as where VðcÞ is the set of vertices of cell c. Vertex values are associated with one third of the cell area, so the factor 1/3 in this case agrees with area weighting and also implies that the area integral of b is preserved in this operation.
We then define the filter by and implement kinematic backscatter via the same formulas as in section 2, This corresponds-in a spatially averaged sense-to a classical dissipative component V and the new antidissipative part −α FV. The viscosity coefficient used for the operator V is given by (12).

Journal of Advances in Modeling Earth Systems
We were able to run numerically stable simulations using (21) with α ≤ 1.5 on eddy-permitting or barely eddy-resolving meshes, that is, with cell dimensions on the order of the first baroclinic Rossby radius. In these simulations, the new scheme led to a noticeable increase in the simulated EKE (see sections 5 and 6). As discussed in Appendix A1, the admissible value of α is expected to depend on the spectral slope of flow EKE, and finding it is a matter of experimenting. This is an obvious complication with using kinematic backscatter in the form (21), for it relies only on the general idea of returning energy to resolved scales, without controlling how much should be returned at a particular location. This means that careful tests of the functioning of the kinematic backscatter parametrization are necessary.

Differences Between Kinematic and Dynamic Backscatter
The new kinematic backscatter is compared to classical viscosity closures as well as the computationally more demanding dynamic backscatter formulation of Juricke et al. (2019), which differs from kinematic backscatter in the following respects: 1. The amplitude of the coefficient α is related to a so-called subgrid energy equation which tracks the amount of dissipated and backscattered energy at each grid point and scales the strength of backscatter according to the available subgrid energy; this adds both a flow-dependent scaling of the backscatter coefficient and a temporal memory to the scheme 2. The viscosity operator is biharmonic, and its coefficient is also dynamically adjusted, based on the Upwind setup described in section 5 below 3. For dynamic backscatter, several iterations of the smoothing filter are applied to the backscatter component (i.e., the second term) in (21) to enhance the scale separation between energy backscatter and viscous dissipation.
The performance of dynamic backscatter was already analyzed and discussed by Juricke et al. (2019Juricke et al. ( , 2020. Dynamic backscatter is computationally more demanding as it (i) applies an additional subgrid energy equation, (ii) uses several iterations of the smoothing filter, and (iii) necessitates a time step reduction for global setups on eddy-permitting meshes. Kinematic backscatter can run with the same time step as the classical viscosity closures which we will use as control in the following sections.

Sensitivity Studies With a Baroclinic Instability Test Case
In this section, we explore the performance of kinematic backscatter in the idealized geometry of a zonally reentrant channel. Our intention is to test and then choose the optimal value of the parameter α, controlling the strength of the backscatter component in (21) and of the scaling coefficient γ 1 in the computation of the viscosity coefficient (12).
We follow the channel setup of Soufflet et al. (2016): The channel extends 500 km in the zonal direction, 2,000 km in the meridional direction and has a uniform depth of 4,000 m. The β-plane approximation with Coriolis frequency 10 −4 s −1 in the center and β ¼ 1:6 · 10 −11 m −1 s −1 is used. Temperature is the only active tracer, and the equation of state for density is linearized with respect to temperature. The channel is initialized with a zonally invariant temperature distribution and a meridional temperature gradient which induces an eastward jet in the center of the channel. The initial zonal velocity field is in geostrophic balance. Jet destabilization is triggered by a small temperature perturbation. To achieve and maintain a quasi-stationary turbulent regime, the zonally averaged temperature and velocity fields are relaxed to the basic state that was used as the unperturbed initial state. The restoring timescale is set to 50 days, as in Soufflet et al. (2016). The Rossby radius of deformation is about 30 km in the center and ranges from about 35 km in the south to 25 km in the north. For more details on the setup, see Soufflet et al. (2016).
The reason for following this setup is (i) that it is well documented and (ii) has a mean flow fixed through relaxation. Because of the controlled mean flow, an increase in resolution and therefore increased eddy-activity does not alter the structure of the mean flow. While eddy-mean flow interactions might be important, our sensitivity studies focus solely on the eddy component of the flow.
This test case exhibits a near-surface submesoscale instability which is starting to be resolved at around 2-to 5-km resolution, changing the dynamics in the channel. Since the focus of this study is on mesoscale instabilities and their improved representation by kinematic backscatter, we will only run sensitivity studies on resolutions of 20 and 10 km. The 20-km resolution setup falls into the eddy-permitting range, which is the range where ocean kinetic energy backscatter has the highest potential for improvements (see, e.g., Juricke et al., 2019). The 10-km resolution is used as an eddy-resolving control as well as a setup to investigate how the kinematic backscatter scales with resolution. As less kinetic energy is removed via dissipation at higher resolution and the mesoscale eddy field becomes largely resolved, the strength of the backscatter component of the kinematic backscatter is expected to decrease as resolution is increased.
For the numerical setup of this study, we use FESOM2. The two different surface meshes are made of equilateral triangles with sides of 20 and 10 km, respectively, with 40 unevenly spaced vertical layers. For the viscous closure, we compare the following cases.
• Two different, classical viscosity closures: harmonic Leith viscosity ("Leith") and a biharmonic viscosity setup which simulates the dissipation of a third-order upwind scheme, with small additional background harmonic viscosity ("Upwind"). These schemes are described in detail in Juricke et al. (2019). • First-order kinematic backscatter as described by (21) with the new viscosity coefficient described by (12) for different values of α between 0 and 1.5. Reducing α reduces the amount of energy that is scattered back on larger scales. For α ¼ 0 the second term in (21) drops out and what is left is a harmonic viscosity operator with the viscosity coefficient (12). We also pick two values for the scaling coefficient, γ 1 , of (12): (i) a setup with a lower coefficient γ 1 ¼ 1 35 ≈ 0:03 ("KBackα") and (ii) a setup with a larger coefficient γ 1 ¼ 0:1 ("KBackα*") which turns out to be more stable at the higher 10-km resolution. Note that lower and larger here characterize only the scaling factor γ 1 and not the actual magnitude of the viscosity coefficient (12). Simulations are run for 10 years where the first year is then omitted as spin-up for all temporal averages. The main diagnostics are vertical profiles of EKE and fluctuations of vertical velocity w (see Figure 2). For each depth, layer-averaged, 9-year temporal means for EKE are calculated. Similarly, w fluctuations are computed as the root mean square of vertical velocity anomalies where brackets denote horizontal layer averaging, the bar denotes time averaging, and w′ ¼ w − w is the temporal anomaly of w. We also discuss the spatial structure of dissipation (see Figures 4 and 5 below) and total amount of energy dissipation. We remark that the vertical profiles of EKE and variance of vertical velocity w simulated with FESOM2 ( Figure 2) and those presented by Soufflet et al. (2016, their Figure 4) for NEMO and ROMS are very similar but not quite the same. These small differences are not relevant for this study.

Results for the 20-km Mesh
We first present our results for the eddy-permitting 20-km mesh. In our discussion, we focus on kinematic backscatter with coefficient α ¼ 1:5 and the smaller scaling coefficient γ 1 ≈ 0.03. A broader range of choices is discussed in sections 5.2 and 5.3 together with the results for the 10-km mesh.
At 20-km resolution, all schemes show a similar structure of the vertical EKE profile, with high levels of energy at the surface and rapidly decaying energy toward the bottom (Figure 2a). The backscatter runs DBack and KBack1.5 are substantially more energetic than the two reference simulations with classical viscosity closures, which have similar levels of EKE between them, with Upwind slightly more energetic than Leith. The relative increase of EKE between the two backscatter simulations and the reference simulations is strongest at depths between 1,500 m and the bottom. Both backscatter schemes are slightly more energetic than the higher-resolution 10-km control simulation with upwind biharmonic viscosity, Upwind10. While the latter is eddy-resolving, its EKE level is still not saturated as the submesoscales are not yet resolved (see the discussion in Soufflet et al., 2016).
Differences in the fluctuations of vertical velocity between the simulations follow a similar pattern as the differences in the EKE profiles (Figure 2b). The qualitative structure of the profile is the same between all simulations, with peaks of w RMS at around 1,000-m depth. The backscatter schemes enhance w RMS . Compared to the higher-resolution control, Upwind10, levels of w RMS lie slightly below for DBack and somewhat above for KBack1.5, especially at depths between 500 to 1,000 m.
We conclude that both backscatter schemes produce a level of EKE close to the control Upwind10.
Fluctuations of vertical velocity are enhanced, but the qualitative structure of the vertical profiles is unaltered.
The increased levels of kinetic energy in KBack are accompanied by reduced total energy dissipation tendencies ( Figure 3). Here, total energy dissipation is defined as the combination of energy injection due to the backscatter component, that is, −α FV in (21), and energy dissipation due to the viscosity component, that is, the operator V in (21). The power injected by the backscatter component is positive definite in an area-averaged sense, while the power due to the dissipation component is negative definite in an area-averaged sense. The energy dissipation tendencies in KBack are 2 to 10 times lower than those of the  (21), that is, u · ð−αFVuÞ, while tendencies of total dissipation correspond to the full operator (21), that is, u · V ð1Þ u. KBack1.5 has the least amount of dissipation and the largest amount of backscatter, followed by KBack1.2 and KBack1.0. Calculations rely on 9 years of data; tendencies of the backscatter component are only available for KBack.

Journal of Advances in Modeling Earth Systems
classical viscosity closures, with largest reductions at depth. As mentioned before, with decreasing α, less energy is scattered back and the total dissipation tendency is larger, that is, KBack1.0 is more dissipative than KBack1.5.
When α > 1, the viscous operator is antidissipative on some range of larger scales. To illustrate this, we can compare the tendencies of energy dissipation based on the original fields with the tendencies based on spatially filtered fields (Figure 4). For the unfiltered, original fields we compute the energy tendency by multiplying the dissipation tendencies with the velocity u. Thus, in the case of KBack, the energy tendency reads u · V ð1Þ u (Figures 4a and 4b). For the spatially filtered counterparts of the tendency of energy dissipation, we apply the filter X (17)-which averages area-weighted cell values to common vertices-to each factor in this product. Thus, in the case of KBack, the spatially filtered tendency of energy dissipation reads Xu · XV ð1Þ u (Figures 4c and 4d). This averaging is conservative for the individual fields but removes small-scale fluctuations. Computing energy tendencies based on filtered fields removes the impact of the grid scales on energy dissipation and highlights the behavior of the energy tendencies on all scales larger than these smallest scales.
The unfiltered tendencies of KBack1.5 are mostly negative when averaged in time, that is, the scheme dissipates energy. However, the corresponding energy tendencies based on filtered fields are predominantly Figure 4. Kinematic backscatter (KBack1.5) acts as mostly negative viscosity on larger scales compared to classical viscosity closures (Leith) which are dissipative also on larger scales: An annual mean (Year 9) of dissipation energy tendency for (a and c) Leith and (b and d) KBack1.5. The dissipation energy tendencies are calculated from original fields on triangles (a and b) and from filtered fields at vertices (b and d). Note that the filtering is performed before the calculation of energy tendencies and temporal averaging, that is, due to nonlinearity (c) cannot be obtained by simply filtering (a).

Journal of Advances in Modeling Earth Systems
positive, that is, the kinematic scheme energizes the flow on larger scales (Figure 4d), while it still maintains total energy dissipation based on the original fields (Figure 4b). The Leith simulation, on the other hand, dissipates energy even on larger scales (Figures 4a and 4c). As α is reduced, the spatially local tendencies based on filtered fields become increasingly mixed between positive (i.e., energy injection) and negative (i.e., energy dissipation) values (not shown). For α ¼ 1 , when kinematic backscatter would reduce to a biharmonic operator if the coefficient of viscosity were constant, most of the filtered fields are negative and therefore predominantly dissipative.
On the daily scale, without long-term temporal averaging, dissipation energy tendencies for filtered fields are mostly positive but also show negative values for KBack1.5 (Figure 5d); the opposite is true for the classical Leith viscosity scheme (Figure 5c). This illustrates that spatially local energy tendencies are not sign definite on daily time scales, with spatial patterns related to the frontal structure of the simulated flow. Even annual averaging does not make the energy tendency of the Leith scheme strictly negative at all locations ( Figure 4a); that is, it does not lead to a sign definite behavior of the dissipation operator everywhere. This sign indefinite behavior was also discussed by Juricke et al. (2019) in their formulation of the dynamic backscatter closure. In dynamic backscatter, energy dissipation tendencies act as a source term for the subgrid energy equation. Juricke et al. (2019) mention two options for the computation of the dissipation tendency: either using the actual energy dissipation tendency computed as the power of the full viscous force or using a sign definite form obtained by eliminating the flux part, which was used by Jansen et al. (2015) and Klöwer et al. (2018). Juricke et al. (2019) observed that using the actual energy tendency leads to higher and more realistic levels of EKE.
(a) (b) (c) (d) Figure 5. As in Figure 4 but for a daily instead of an annual mean.

Results for the 10-km Mesh
We now compare the different schemes for the same channel simulation on a 10-km mesh which is essentially eddy-resolving. In particular, we ask how the parameters of the scheme depend on resolution. As we have not built any explicit resolution scaling of the coefficient α into the derivation of the scheme, we would not expect kinematic backscatter to be scale-aware in a strict sense, but the extent of necessary retuning is of interest.
Indeed, a direct application of kinematic backscatter with parameters from the 20-km simulation shows an unphysically strong increase of w RMS in deep layers, around the depth of 3 km (see Figure 6b, and compare to Figure 4 of Soufflet et al., 2016). This increase in variance is due to waves propagating from the upper layers and is presumably due to the excitation of the Charney instability mode, which is not fully resolved on a 10-km mesh but is also not damped enough in some regions due to the minimal choice of the coefficient γ 1 in the computation of the viscosity coefficient (12). An increased coefficient, γ 1 ¼ 0:1, will filter grid-scale differences more strongly and largely remove the artifact in the profile of w RMS while having little impact on the profile of mean EKE (profiles labeled KBack1.5* in Figure 6).
Comparing kinematic backscatter to the classical viscosity parametrizations, we see that the relative increase in EKE provided by KBack1.5* (as well as KBack1.5) to Upwind is not as high as for KBack1.5 on the 20-km mesh. This is expected as we move from eddy-permitting to eddy-resolving resolution and the relative effect of the backscatter component decreases. While the EKE increase due to KBack1.5* is slightly lower than the increase provided by DBack, it is even slightly larger than the increase due to KBack1.5. At 10-km resolution we need to discard the latter due to the artifacts in its w RMS profile.

Choice of Parameters for Kinematic Backscatter
Prompted by the difference in behavior on the two meshes, we now carry out a limited study of the influence of the parameters α and γ 1 on the performance of kinematic backscatter. While the strength of the backscatter component is controlled by α, the new viscosity coefficient (12) can be scaled by the coefficient γ 1 and may on its own lead to reduced dissipation. We therefore consider our main diagnostics for α ¼ 0 (only the new viscosity coefficient is active), α ¼ 1 (approximately biharmonic behavior), and α ¼ 1:5 (well into the backscatter regime), for the two choices of γ 1 and for the two different meshes (Figure 7).
For γ 1 ≈ 0.03 at 20-km resolution, the largest contribution to energizing the flow comes from the new viscosity coefficient (12). Indeed, even with α ¼ 0, EKE is substantially increased compared to Upwind. A subsequent increase in the backscatter coefficient α leads to a further noticeable increase in EKE but weaker than the initial increase due to the new viscosity. 10.1029/2020MS002175

Journal of Advances in Modeling Earth Systems
For γ 1 ¼ 0:1 , in contrast, KBack0.0* has levels of EKE that are very close to the classical Upwind, and KBack1.5* reaches nearly the same level as KBack1.5. In this case, the increase of EKE is almost entirely due to the backscatter component in (21). The profiles of w RMS (Figures 7b and 7d) indicate that using a smaller scaling factor γ 1 ≈ 0.03 in the computation of the viscosity coefficient (12) in KBack leads to increased variance, especially at 10-km resolution. This can be seen as a sign that the smaller scaling coefficient γ 1 is close to the minimum that still provides realistic results at 20-km resolution and that it is insufficient at 10 km.
We highlight that kinematic backscatter is expected to be dissipative, on average, at all scales when α ≤ 1. But already the choice α ¼ 1 (corresponding to a biharmonic operator if the viscosity coefficient were constant) leads to a considerable increase of EKE for both values of γ 1 at 20-km resolution as compared to the simulation with α ¼ 0.
Furthermore, we remark that even though the scaling coefficient γ 1 and, consequently, the viscosity coefficient are larger for KBack1.5* compared to KBack1.5 at the same resolution, γ 1 also affects the scales at which the scheme scatters back. Thus, an ordering of the level of EKE by the value of γ 1 cannot be In summary, the relative effect of the backscatter component of (21) depends on the magnitude of the coefficient γ 1 , which is part of the viscosity coefficient (12), and the mesh resolution. The numerical stability even for very small values of γ 1 (e.g., γ 1 ≈ 0.03) is, at least partly, due to the spatially very localized viscosity coefficient (12) helping to efficiently remove noise at the grid scale. In addition to reducing total viscosity, the backscatter component further energizes the flow, even though the increase in EKE is not necessarily large.
However, a small value of γ 1 is not always a safe option, especially at higher resolution; for more general applications, a larger value of γ 1 is needed, as in KBack*. In that case, the influence of the backscatter component is much more noticeable as compared to the effect of the new viscosity coefficient. At 10-km resolution, EKE levels of the simulation KBack* with a larger coefficient γ 1 are similar to those of KBack ( Figure 6) but without the imprint of spurious waves in the profiles of w RMS . At eddy-permitting 20-km resolution, all options for kinematic backscatter produce more realistic EKE profiles and reduce energy dissipation while maintaining the general flow structure. Consequently, kinematic backscatter configurations with the strongest increase in EKE, KBack1.5 and KBack1.5*, are both safe options at eddy-permitting resolution. For

Journal of Advances in Modeling Earth Systems
higher resolutions that vary between eddy-permitting and eddy-resolving, we recommend to use the kinematic backscatter setup with the larger value of γ 1 , KBack1.5*.

Global Ocean Simulations at Eddy-Permitting Resolution
In a second set of experiments, we investigate the behavior of kinematic backscatter in realistic simulations: FESOM2 in a global configuration on a mesh with 1/4°nominal resolution. The mesh is derived from the NEMO ORCA25 mesh by splitting its scalar cells into two triangles. Its resolution varies from about 25 km around the equator to about 10 km in high latitudes. This resolution is called eddy-permitting: For a large portion of the globe, it does not allow eddies to be simulated explicitly, but the mesh triangles are smaller or close to the Rossby radius in many places so that eddies can be represented formally. The same mesh was used to test the dynamic backscatter parametrization DBack in Juricke et al. (2020).
We performed two experiments. The first uses the standard FESOM settings described in Scholz et al. (2019) with Leith viscosity ("Control" in the following). The GM parametrization is active in this experiment, but only a very small background thickness diffusion is present in elements smaller than 25 km (i.e., almost everywhere). The second experiment uses the new viscosity coefficient combined with kinematic backscatter with α ¼ 1:5 and γ 1 ≈ 0.03, and the GM eddy parametrization is switched off entirely. The parameters are as for KBack1.5 in the previous section, so that we retain this name here. Both experiments use COREII forcing (Large & Yeager, 2009).
For both experiments, we first performed a spin-up over the years 1979-2009. Each run was then continued over the complete COREII cycle , using the spin-up results as initial conditions. Only the final years 1979-2009 from the full COREII cycle are used for comparison with observational estimates.
The snapshots of the absolute velocity (Figures 8a and 8b) show that, in general, KBack1.5 has more energetic currents, with larger velocities even over "calm" areas in the centers of the ocean basins. Currents in the Southern Ocean show more spatial variability, and well-defined eddies form at the Agulhas retroflection in KBack1.5, while Control does not show Agulhas eddies in the southern Atlantic.
The temporal standard deviation of sea surface height ("SSH variability") is often used as a measure of eddy variability (e.g., Juricke et al., 2020;Sein et al., 2017). FESOM2, as well as many other ocean models, tends to considerably underestimate SSH variability in eddy-permitting configurations (e.g., Penduff et al., 2010;Sein et al., 2017). As mesoscale eddies cannot be fully resolved at eddy-permitting resolution, SSH variability is especially underestimated in regions of strong mesoscale eddy activity, that is, in the Southern Ocean and the western boundary currents (Sein et al., 2017). Thus, in general, an increase in variability will be an improvement when compared to observational estimates of AVISO (https://www.aviso.altimetry.fr Ducet et al., 2000;Le Traon et al., 1998). The use of kinematic backscatter combined with the new viscosity coefficient leads to a considerable increase of SSH variability over most energetically active areas of the oceans (Figures 8c-8e). This is especially noticeable in the Southern Ocean, the Gulf Stream, and the Kuroshio, with most pronounced changes for the Kuroshio and the Agulhas currents, where also the underestimation of Control is largest. The increase in these regions is over 0.1 m, which is substantial, but still not enough to reach the levels observed by satellite altimetry (cf. Figures 8c and 8d with 8f).
There are clear differences between the two experiments in terms of temperature biases with respect to the PHC (PHC 3.0, updated from Steele et al., 2001) climatology (Figure 9). At the surface, the infamous cold bias over the northwest corner of the Atlantic Ocean (see, e.g., Wang et al., 2017) is considerably reduced in KBack1.5, although it does not disappear completely. This is also true for the same bias at 100-m depth. Eddy-mean flow interactions play an important role in defining the path and the intensity of the Gulf Stream (e.g., Kang & Curchitser, 2015). The enhanced eddy activity in KBack1.5 leads to an improved simulation of the path of the Gulf Stream and its extension (cf. the discussion in Juricke et al., 2020). Furthermore, we observe a general decrease of the 100-m warm bias in the northern North Atlantic. There are also reductions of biases at 100-m depth along the Southern Ocean fronts, especially in the Atlantic sector, while the Gulf Stream gets warmer. At 500-m depth, the most striking difference is a reduction of the biases in the whole Southern Ocean. Temperatures in the North Atlantic also move closer to climatology, while differences in the Pacific Ocean are less pronounced. The cold bias at 1,000 m west of the Strait of Gibraltar largely disappears in KBack1.5, and the warm bias in the Southern Ocean is reduced.
Overall, kinematic backscatter with the new viscosity coefficient leads to improvements in the ocean dynamics that translate into reduced hydrographic biases. The bias reductions are similar to those obtained with the dynamic energy backscatter parametrization of Juricke et al. (2020) and are associated with enhanced eddy activity and associated changes in eddy heat transport. However, different from the dynamic backscatter of Juricke et al. (2020)-which required solving an additional subgrid energy equation as well as a time step reduction to 10 min for the global simulations-kinematic backscatter does not induce extra costs in terms of computation time. Here, we can keep the same, relatively large time step of 20 min for both Control and KBack1.5. Further, our implementation of kinematic backscatter comes at no additional cost for the computation of the operators. As KBack1.5 runs stably with the larger time step despite enhanced eddy activity, we conjecture that the time step reduction necessary for dynamic backscatter on global meshes in Juricke et al. (2020) is primarily due to an intermittently excessively large negative viscosity coefficient there, not due to the more vigorous eddy activity. How to transfer this positive result to dynamic backscatter requires further studies.

Journal of Advances in Modeling Earth Systems
Improvements with KBack1.5 compared to the previous FESOM2 default setup Control with Leith viscosity are substantial without additional cost and obvious drawbacks. While DBack may energize the flow even more, as discussed in section 5, this additional increase in EKE does not yet justify the increase in computational cost. However, further adjustments to DBack may change this assessment and ultimately challenge KBack1.5, especially since DBack automatically scales with resolution without additional tuning.

Discussions and Conclusions
We describe a simplified, "kinematic" implementation of an energy backscatter parametrization, together with a new flow-dependent parametrization for the viscosity coefficient. The new scheme is tested in an idealized channel and in a global mesh configuration. In both tests, it leads to considerable improvements in simulated ocean dynamics as well as hydrography in comparison with classical viscosity closures.
Viscosity affects a range of scales in the vicinity of the grid scales, commonly creating excessive dissipation on eddy-permitting meshes. Recognizing this, we propose to reduce the viscous force by subtracting its locally averaged value multiplied with some coefficient. In this case, viscous damping remains unchanged at grid scales but is reduced at larger scales or replaced by antidissipation.
We reiterate that the viscous operator with its new local coefficient of viscosity discussed here is not a physical viscosity in a strict sense. Rather, we regard it as a numerical filter which prevents accumulation of grid-scale noise. It preserves total momentum and dissipates kinetic energy but violates more subtle physical principles. However, this occurs at small scales where all operators deviate from their continuous counterparts. The utility of our approach is supported by the simulations shown in this paper.
We describe the implementation of kinematic backscatter for cell velocities on a triangular mesh and show that, in practice, subtraction of the mean viscosity multiplied with a coefficient larger than one (up to 1.5 in the case of FESOM2) is possible, which allows the simulated flows to be reenergized. In contrast to dynamic backscatter (see, e.g., Juricke et al., 2019), kinematic backscatter does not use a consistent subgrid energy budget but rather relies on the heuristics that (i) net dissipation must be reduced on eddy-permitting meshes, (ii) dissipation near the grid scale cannot be reduced for reasons of numerical stability, and (iii) nonlinear interactions will act to restore the cascade of energy across scales.
Kinematic backscatter uses a simple new viscosity coefficient based on velocity differences between the nearest neighbors. The coefficient can, therefore, be efficiently computed together with the necessary operators, resulting in a slight reduction in computational costs compared to the Leith scheme with its more complex viscosity coefficient. More importantly, kinematic backscatter comes at much lower computational cost compared to dynamic backscatter for two reasons. First, the current implementation of dynamic backscatter requires a time step smaller than is necessary with kinematic backscatter. Second, kinematic backscatter does not employ an additional prognostic equation for subgrid energy. While this is computationally less expensive, it has the drawback that there is currently no theory for the selection of the factor α which scales the backscatter component of the new scheme. The selection may depend on the shape of the energy spectrum of the simulated flow, thus, at least partially, on the energy scattered back to larger scales by the scheme itself for values of α > 1. We were able run eddy-permitting simulations stably with values for α up to 1.5, in the channel as well as for the global ocean. Adjustments might be required on different meshes.
The EKE increase provided by our kinematic backscatter parametrization depends on the strength of the chosen viscosity. The amplitude of the viscosity coefficient for kinematic backscatter can be controlled by a scaling coefficient γ 1 . The eddy-permitting configurations (the channel at 20 km and the realistic ocean at 1/4°mesh) run stably with a small γ 1 , that is, γ 1 ≈ 0.03. It creates much less dissipation than traditional viscosities such as the Leith scheme. This reduced dissipation contributes to the EKE increase for the 20-km resolution channel even more than the backscatter component. However, such a weak viscosity leads to a spurious behavior of the vertical velocity component w on a finer mesh (10-km resolution), and a higher γ 1 ¼ 0:1 is needed to eliminate it. For the viscosity coefficient with γ 1 ¼ 0:1, kinematic backscatter simulates nearly the same EKE levels as reached in kinematic backscatter runs with the smaller γ 1 . In these simulations with increased γ 1 , any increase in EKE is now almost entirely due to the backscatter component. Ultimately, EKE can be increased by tuning both the viscosity coefficient and the backscatter component of the kinematic backscatter scheme. The eventual increase in EKE is quite insensitive to the choice of the coefficient γ 1 when α is set to 1.5. We can conclude that with α ¼ 1:5, the increased dissipation due to a larger γ 1 is effectively compensated by similarly increased antidissipation, and a larger γ 1 only adjusts the minimum dissipation near the grid scale, which is necessary for model stability. In the future, more testing is required for all possible options for the viscosity and backscatter components (including biharmonic viscosities) to further optimize the scheme, by also taking into consideration additional aspects, such as spurious dianeutral mixing that might accompany an increase in EKE (Ilıcak et al., 2012). An analysis of the behavior of kinetic energy and energy dissipation close to the grid scales may provide further insights into the selection of γ 1 and α. The placement of velocities at triangle centers in FESOM2 creates aliasing in spectral analysis at grid scales, which is the reason why dissipation spectra were not presented here. The development of alternative methods for scale analysis is the subject of ongoing work.
Finally, we stress that energy backscatter parametrizations have a limited range of resolutions where they provide substantial benefits, namely, flows on eddy-permitting or barely eddy-resolving meshes. There is no strong argument in favor of using backscatter if the scales of EKE dissipation are well separated from the scales of EKE generation. However, kinematic backscatter with α ≤ 1 may still be useful as a means to combine harmonic and biharmonic viscosities in those cases.
Eddy-permitting models will remain an essential part of coupled climate models for years to come, and in many applications fully eddy-resolving ocean simulations will remain unaffordable. This new kinematic backscatter parametrization provides an efficient and easy-to-implement option to improve the simulation of eddy effects in global ocean simulations and bring eddy-permitting simulations closer to computationally much more expensive higher-resolution setups.