Ocean Kinetic Energy Backscatter Parametrization on Unstructured Grids: Impact on Global Eddy‐Permitting Simulations

In this study we demonstrate the potential of a kinetic energy backscatter scheme for use in global ocean simulations. Ocean models commonly employ (bi)harmonic eddy viscosities causing excessive dissipation of kinetic energy in eddy‐permitting simulations. Overdissipation not only affects the smallest resolved scales but also the generation of eddies through baroclinic instabilities, impacting the entire wave number spectrum. The backscatter scheme returns part of this overdissipated energy back into the resolved flow. We employ backscatter in the FESOM2 multiresolution ocean model with a quasi‐uniform 1/4° mesh. In multidecadal ocean simulations, backscatter increases eddy activity by a factor 2 or more, moving the simulation closer to observational estimates of sea surface height variability. Moreover, mean sea surface height, temperature, and salinity biases are reduced. This amounts to a globally averaged bias reduction of around 10% for each field, which is even larger in the Antarctic Circumpolar Current. However, in some regions such as the coastal Kuroshio, backscatter leads to a slight overenergizing of the flow and, in the Antarctic, to an unrealistic reduction of sea ice. Some of the bias increases can be reduced by a retuning of the model, and we suggest related adjustments to the backscatter scheme. The backscatter simulation is about 2.5 times as expensive as a simulation without backscatter. Most of the increased cost is due to a halving of the time step to accommodate higher simulated velocities.


Introduction
The accurate representation of the oceanic mesoscale eddy field is still a major challenge for current global ocean models, especially those that are used as ocean components of global climate models. Most of the kinetic energy in the ocean is contributed by the mesoscale eddy field (Ferrari & Wunsch, 2009). However, due to computational constraints, many current global ocean models intended for climate simulations still cannot afford resolutions at which mesoscale eddies and their life cycles are sufficiently well simulated. These models are usually either coarse with a resolution of about 1 • and lower (Taylor et al., 2012), in which case, all eddy effects are parametrized, or they are eddy permitting with a resolution between about 1/3 • and 1/10 • (Haarsma et al., 2016). To simulate eddies with some degree of realism, a resolution of 1/10 • or even higher-especially in the very high latitudes and in some areas on the continental shelves-is required: Hallberg (2013) suggests that a resolution with at least 2 grid points per deformation radius is needed, but Sein et al. (2017) demonstrated that even this is not necessarily sufficient. Thus, fully resolved simulations of the eddy field are computationally extremely costly, usually too costly to run them as part of climate models for decades to centuries. As a result, ocean eddy effects are generally underrepresented in climate For non-eddy-resolving models, the most common eddy parametrization-aside from isopycnal diffusivity-is the Gent-McWilliams (GM) parametrization (Gent & McWilliams, 1990) where an additional advection term in the tracer equation represents the mean eddy effect as a potential energy sink. For finer-scale eddy-permitting simulations, viscosity schemes in the momentum equation (e.g., Fox-Kemper et al., 2014) are generally parametrizing unresolved subgrid eddy effects, assuming-on average-a dissipative nature of these structures. Such closures are necessary to ensure model stability by preventing buildup of enstrophy at the grid scale (see, e.g., Danilov et al., 2019). They also dissipate an excessive amount of energy which-at eddy-permitting resolution-not only affects the small scales but can also reach up to the scales where eddies are generated. The consequence is a reduction in effective resolution and mean kinetic energy (MKE) as the eddies themselves and the baroclinic instabilities through which they are generated are damped (Soufflet et al., 2016).
A way to improve the representation of eddies without substantially increasing the model grid resolution is the use of unstructured-mesh models such as the Finite Element Sea ice-Ocean Model (FESOM) (Danilov et al., 2004;Wang et al., 2014). It was the first multiresolution global ocean model which was coupled to an atmosphere model, enabling climate simulations (Rackow et al., 2018;Sidorenko et al., 2015). By local refinement of resolution on triangular meshes, computational resources can be focused on regions where explicitly resolving important physical processes is expected to enhance the quality of the simulations, leaving the rest at coarser resolution (Sein et al., 2016;Sein et al., 2017). The effects of an explicit simulation of eddies on the representation of model biases has already been investigated with the variable meshes supported by FESOM, for example, in the Southern Ocean (Rackow et al., 2019) or the Fram Strait (Wekerle et al., 2017). With FESOM2, the dynamical core has been changed to a finite volume scheme Koldunov et al., 2019;Scholz et al., 2019), yielding much higher computational efficiency. However, even with locally increased resolution in simulations on a variable grid, impacts of a region with a coarser mesh are still observable in higher-resolution domains downstream (Danilov & Wang, 2015;Sein et al., 2016).
In this context, ocean kinetic energy backscatter gained a lot of attention in recent years, so far mostly explored in idealized settings, for example, for idealized two-layer quasi-geostrophic (Jansen & Held, 2014) or primitive equation configurations (Jansen et al., 2015) or for the shallow water equations (Klöwer et al., 2018). Backscatter reduces overdissipation by reinjecting a portion of the excessively dissipated kinetic energy into the resolved flow. Reinjecting this energy on larger scales feeds the inverse energy cascade which is partly resolved in eddy-permitting simulations. As a result, backscatter improves the representation of kinetic energy across the entire spectrum, hence, in combination with a classical viscosity closure, provides a more energetically consistent subgrid model which not only improves eddy variability but also the feedback of eddy activity on the mean flow. As the physical rate of enstrophy dissipation is maintained, backscatter proves to be numerically stable despite its antidiffusive behavior on intermediate scales. Juricke et al. (2019) successfully implemented backscatter in the full primitive equation FESOM2 model in an eddy-permitting channel setup. The present study is a natural extension of this work. We demonstrate that backscatter is viable for a global eddy-permitting ocean model. The backscatter scheme is identical to the one used by Juricke et al. (2019), with only very minor adjustments to two parameters, discussed below. Our analysis focuses on eddy variability, which backscatter is expected to improve, as well as on the global water mass and circulation structure, where one needs to ascertain that the simulation of water mass properties is not deteriorated through potentially overly strong mixing.
Our results not only show a clear improvement of eddy variability but also of mean SSH and general water mass properties. For the latter, the improvements are not as uniform, and we observe deterioration in some areas, such as the Antarctic coast line and continental shelves. However, reductions in climatological biases

Journal of Advances in Modeling Earth Systems
10.1029/2019MS001855 dominate. We also perform an initial sensitivity study to assess the impact of the choice of parameters and the potential for additional tuning.
The structure of this paper is as follows. Section 2 introduces the model and the experimental setup, consisting of two multidecadal simulations, a reference without backscatter, and a second simulation with backscatter, performed on a 1/4 • mesh. In section 3, we discuss the quality of the eddy-permitting simulation and the changes caused by the backscatter parametrization, focusing on biases in SSH variability, mean SSH, temperature, and salinity and also on general aspects of the simulated circulation and sea ice. In section 4, we estimate the sensitivity of the backscatter run to several tuning parameters. Section 5 summarizes our results and discusses their relevance and future perspectives for global ocean models.

Model Setup
For this study we use the current version of FESOM2 . FESOM2 is a finite volume ocean and sea ice model that uses a triangular surface grid. It utilizes a cell-vertex (quasi-B grid) discretization where horizontal velocities are computed on the centroids of triangles and the tracers, SSH, and vertical velocities are computed on vertices. Further details on model setup as well as temporal and spatial discretizations can be found in Danilov et al. (2017) and Scholz et al. (2019). The implementation of backscatter follows Juricke et al. (2019).

Experimental Setup
Two main simulations were carried out, one reference simulation without backscatter, denoted REF in the following, and one simulation with backscatter switched on, denoted BACK1. Both simulations use the same mesh derived from the NEMO model ORCA0.25 mesh (Bernard et al., 2006) by regularly splitting its quadrilateral cells into triangles, with a nominal resolution of approximately 1/4 • and 46 vertical layers. The simulations were started from a control simulation on 1 January 1979 and run for 31 years and use the same atmospheric forcing provided by COREII (Large & Yeager, 2009) which is available for 1948-2009. The control simulation is a spin up from 1948, where initial temperature and salinity data were provided by PHC (PHC 3.0, updated from Steele et al., 2001) and initial velocities were set to zero.
Simulation REF is an extension of the control and uses the same parameter settings: a time step of 20 min and the modified harmonic Leith viscosity with biharmonic background viscosity (see Fox-Kemper & Menemenlis, 2008, and the appendix of Juricke et al., 2019, for more details). The only difference to the control is that the GM eddy parametrization (Gent & McWilliams, 1990) for non-eddy-resolving grids was switched off entirely in REF, while the control used the resolution-scaled GM parametrization following Wang et al. (2014). Redi (1982) diffusion is used for the tracer equations. To parametrize vertical mixing, the KPP scheme was selected (Large et al., 1994).
Simulation BACK1 uses the default backscatter parametrization developed by Juricke et al. (2019). This parametrization maintains a scalar unresolved kinetic energy (UKE) budget for each grid cell. A large fraction of the energy that is locally dissipated by the viscosity operator enters the UKE budget as a source. UKE is subject to horizontal diffusion. Moreover, energy can be returned from the UKE budget into the resolved scales of the flow via an antidiffusive operator, that is, a spatially smoothed approximate negative Laplacian whose coefficient is calculated locally from the amount of available UKE.
In our scheme, the spatially smoothed Laplacian backscatter operator is paired with a biharmonic viscosity operator which differs from the one used in REF. The viscosity coefficient is scaled by the local absolute velocity; for details on the backscatter scheme and the justification for the choice of the viscosity and backscatter operators, see Juricke et al. (2019) and Appendix A. The fraction of dissipated energy entering the UKE equation is related to the local Rossby number. Compared to Juricke et al. (2019), the scaling parameter R dis which selects the local Rossby numbers beyond which backscatter is reduced has been changed from 1 to 0.5. This reduces the amplitude of the backscatter slightly and was found necessary to ensure model stability. In addition, the time step for BACK1 had to be reduced from 20 to 10 min, also for stability reasons, as otherwise the increased velocities in the backscatter run would lead to CFL violations. There is no special local grid refinement in this setup except for that of the NEMO mesh. This also means that the simulations do not take advantage of the ability of FESOM2 to locally increase resolution where necessary. Futhermore, REF was not tuned for production and may possess larger biases than if it were tuned. The reason is that the 1/4 • mesh is not the one commonly used with FESOM2 but was specifically made for this study to show the results in an eddy-permitting configuration that is also in use by other current state-of-the-art ocean models (e.g., Bernard et al., 2006;Haarsma et al., 2016).
In addition to the main runs, we carried out a set of sensitivity experiments where, starting with the parameter settings of BACK1, we individually perturbed each of a set of main tuning parameters. Motivation and details are discussed in section 4. The runs for this sensitivity study were performed over a 6 year horizon, split off from the main simulation BACK1 starting 1 January 2000. The short integration period is sufficient to get a first impression of the impact of parameter changes to the backscatter run.

Diagnostics
All diagnostics are based on monthly means. Using monthly means is sufficient to highlight changes in simulated SSH variability due to the backscatter parametrization. The first year of simulation, 1979, is excluded to account for the initial spinup of BACK1. This is long enough to remove most of the model drift in the upper ocean velocity field due to the change of the eddy parametrization. However, neither simulation, especially the backscatter run, will be equilibrated at depth, especially with respect to temperature and salinity. This is discussed in more detail in section 3.5.
We compare both simulations against two different data sets. The PHC data set (Steele et al., 2001) is used to diagnose temperature and salinity biases, both for global maps at different depths as well as for time-averaged transects in the Atlantic, Pacific, and Indian Oceans. In all cases, the PHC climatology is compared to the model time average over the Years 1980-2009. The AVISO data (http://www.aviso.altimetry.fr; Ducet et al., 2000;Le Traon et al., 1998) is used to assess the accuracy of SSH variability. Temporal standard deviation of monthly mean data for the time period 1993-2009 is computed for both model and AVISO data and compared to each other.
The time-averaged temperature and salinity transects used for the comparison with PHC data are located in the Atlantic (30 • W), Pacific (130 • W), and Indian Oceans (80 • E). They are computed using linear interpolation from the nearest nodal values. We refrain from showing climatic trends in this study as within the 30 years of integration both simulations still experience strong model drift. The impact of backscatter on trends in the ocean will need to be studied with longer simulations in future work.
In the analysis of the circulation, MKE and eddy kinetic energy (EKE) are defined as the temporal mean over (u 2 m + v 2 m )∕2, with u m , v m the monthly mean zonal and meridional velocities (for MKE) or monthly mean zonal and meridional velocity anomalies (for EKE), respectively. The anomalies are computed with respect to the time-averaged seasonal cycle.
To analyze changes between REF and BACK1, we compute global overturning as well as the overturning in the Atlantic from the monthly mean vertical velocities for the respective regions. Additionally, mixed layer depth (MLD) is used to diagnose changes in the near-surface structure of the ocean. We calculate MLD as the smallest vertical distance to the surface for which the vertical buoyancy derivative is equal to a "local critical buoyancy gradient" following Griffies et al. (2009). We compare the results to a second definition of MLD: the depth for which the vertical density profile differs by 0.125 sigma units compared to the surface density (also following Griffies et al., 2009). The latter tends to be deeper than the former, but both diagnostics agree in the locations of deep mixed layers and changes in MLD due to backscatter. Thus, the analysis in this paper is carried out using the first definition only. We focus on the annual maximum MLD, calculated from monthly means, and averaged over all years. Finally, we analyze changes in sea ice thickness as sea ice responds sensitively to changes in water masses and circulation.

Results of the Main Simulations
In this section, we analyze the performance of REF and BACK1 against the observational estimates and highlight the impact of the backscatter scheme relative to REF with respect to each of the diagnostics outlined in section 2.3.  AVISO and REF. +100% indicates that AVISO variability is twice as large as the simulated variability. −50% means that AVISO variability is 50% less than variability in the respective simulation.

Impact on Eddy Activity: SSH and Horizontal Velocities
We first look at the variability of SSH as a good indicator for eddy activity. Model performance in SSH variability can be readily compared to satellite observational estimates from AVISO. A common bias in ocean models with eddy-permitting resolution is the underestimation of SSH variability.
For REF, this underestimation can be seen almost everywhere (cf. Figure 1c with 1b directly or see Figure 1e for the difference field). In eddy active regions such as the Southern Ocean and western boundary currents, REF lacks variability by a factor of 2 to 3 or even more. Furthermore, the location of eddy variability is also generally biased: The Kuroshio and Gulf Stream experience a coastal intensification, and the missing northwest corner of the Gulf Stream extension in the North Atlantic causes a too easterly flow of warmer waters toward the north. For the Agulhas Current, the Agulhas retroflection is too weak and reaches too far into the Atlantic. The East Australian current is too coastal and also too weak.
With backscatter, the picture changes considerably. BACK1 shows a strongly improved variability throughout the oceans (Figures 1a and 1d). In many regions such as the Southern Ocean, the southern Indian Ocean, the East Australian current, the Agulhas current and retroflection, and large parts of the North Pacific and North Atlantic, variability biases are substantially reduced (by 50% and more; cf. Figure 1d with 1e). In a tropical band between 10 • N and 10 • S, biases do not change much since the Rossby radius, and therefore, eddy effects here are generally well resolved, while increase in SSH variability is especially pronounced in the midlatitudes. However, even in the tropics, the general tendency is an improved representation of SSH variability with BACK1, for example, off the east African coast.
In a few regions, backscatter changes the SSH variability bias from underestimation to overestimation. Examples are close to and west of New Zealand, parts of the Labrador Sea, and close to the coastlines at 30 • N where there is a misplaced coastal intensification of the Gulf Stream and the Kuroshio. However, these regions are substantially fewer and the amplitude of the overestimation is generally smaller than the underestimation in REF. The overestimations of SSH variability generally point to biases in the mean currents that could not be rectified solely by the increased eddy variability of the backscatter scheme. Since the backscatter parametrization intensifies regions where eddy activity is already present, it is not surprising that it cannot fully remove biases in strong, biased mean currents. While increased eddy activity can lead to shifts and improvements in the mean currents, this is not entirely the case for the Kuroshio and Gulf Stream, although some improvement can be seen downstream off the coast for both currents. This is especially true for the northwest corner of the Gulf Stream extension. While the northwest corner was hardly visible in REF, it starts to appear in BACK1 (cf. Figure 1a and 1c), although the too-easterly branch is still the dominant pathway of the current (observe the dipole pattern in the North Atlantic in Figure 1d).
Surprisingly, a strong positive bias of SSH variability around 140 • W and 60 • S in the South Pacific seen in REF is nearly absent in BACK1, even though backscatter generally tends to increase variability. This points to changes in the mean current at that location which shifted eddy activity away from this region, thereby improving the model simulation.
The general improvements in model performance are summarized by ratios of root mean square error (RMSE) between BACK1 and REF in Table 1. BACK1 reduces global SSH variability biases by around 10% compared to REF. Regionally averaged improvements in (parts of) the Southern Ocean are between 30% and 50%, while the improvements over the Gulf Stream area are a modest 5%, and for the Kuroshio region, we actually observe a slight error increase of less than 1%. In the latter two regions, some local biases have changed sign compared to REF. With respect to SSH variability, BACK1 performs much closer to higher-resolution simulations with an even smaller time step (cf. to, e.g., Figure 3 of Sein et al., 2017). However, a direct comparison of different resolutions with and without backscatter regarding various biases will be left for future studies.
In accordance with changes in SSH variability, also the mean state of SSH is globally improved when using backscatter (    unrealistically increased meandering south of Japan that is not in accordance with AVISO and leads to the general deterioration of the area averaged RMSE when comparing BACK1 with REF.
In summary, all global and regional RMSE for mean SSH and SSH variability-with the exception of the Kuroshio region-are improved. Using just the SSH diagnostics as a decision baseline, BACK1 is a superior simulation than REF.
For example, snapshots of the horizontal velocity field at 100 m depth (00:00, 1 January 1999, Figures 3a and 3b) show clearly that the flow field for BACK1 contains much more kinetic energy compared to REF.
Western boundary currents such as the Kuroshio and the Gulf Stream tend to meander and vary over much larger scales and much more intensely in BACK1 than in REF (cf. Figure 3a and 3b). Other strong currents, such as the Norwegian Currents or boundary currents in the Weddell sea, are more localized and intensified in BACK1. This shows that the backscatter is able to destabilize as well as stabilize boundary currents, both possible eddy-feedback mechanisms (e.g., Hughes & Ash, 2001). Eddy intensification is also visible in MKE and EKE (Figures 3c-3f). Both fields increase substantially when backscatter is switched on. This is even more pronounced in velocity snapshots for specific regions: the Agulhas Current (Figures 4a and 4b), the Malvinas Current (Figures 4c and 4d), and the Kuroshio Current (Figures 4e and 4f). In each case, BACK1 has considerably more eddies. Furthermore, backscatter increases flow variability and enables the development of much finer spatial structures. This effect is especially visible away from the strong mean currents. In BACK1, the Agulhas and Malvinas currents shed a substantial number of large-scale eddies. In REF, these are largely absent (also cf. MKE between Figures 3e and 3f). Backscatter causes an eastward shift of the Agulhas retroflection (cf. Figure 1a with 1c, 2a with 2c, and 3c with 3d). As discussed for changes in mean SSH and SSH variability, BACK1 tends to generate larger meanders of the mean flow in the Kuroshio region which are related to the aforementioned overestimation of SSH variability south of Japan (see Figures 1d  and 2g).

Impact on Water Mass Properties: Temperature
We now ask how the intensified variability and changes in pathways of currents in BACK1 affect the vertical water mass structure of the ocean, especially with respect to biases in temperature and salinity.
Generally, REF shows classical temperature biases at different depths ( Figure 5). There is a strong surface and near-surface cold bias in the central North Atlantic surrounded by a strong warm bias, a surface warm bias along the eastern basin boundaries close to Africa and North and South America, a subsurface and deep warm bias in the Southern Ocean, strong biases in the tropics at 100 m depth, a warm bias in the deep Atlantic, and strong dipole biases along the western boundary currents at the surface and near surface related to wrongly simulated current pathways.
With backscatter, some of these biases are substantially reduced. The change in the Gulf Stream extension reduces the North Atlantic cold bias as well as parts of the surrounding warm bias (Figures 5i and 5j). In the Kuroshio region, the effect is not as clear, with an increased bias in the coastal region where backscatter overexcites SSH variability (cf. Figures 5i and 5j with Figure 1). The signal in the surface and near-surface Southern Ocean is mixed, but generally, backscatter tends to reduce biases. This is especially true for the middepth Southern Ocean, where a dipole pattern counteracts the previous warm-cold biases.
Biases in REF at the surface and down to 1,000 m have magnitudes of around 3-5 K. The magnitude of changes due to backscatter are around 1-2 K, most pronounced at 100 m depth. But even in the deep ocean at around 2,000 m, changes due to backscatter can be as large as 1 K, comparable to model biases at that depth.
Generally, global RMSEs in BACK1 are lower than in REF for all four depth levels ( Table 1). The strongest reduction of about 10% is at the surface while at depth improvements only amount to a few percent.
Biases along transects in the Atlantic, Pacific, and Indian Oceans better highlight the impact of backscatter throughout the full depth of the oceans (Figure 6). In all three transects, REF has a subsurface warm bias in the Southern Ocean between around 40 • S and 60 • S, a cool bias down to around 500 m between 40 • S and 20 • S, and-with the exception of the Pacific-a warm bias down to around 1,000 m closer to the equator. There is also a near-surface warm bias closer to the Antarctic continent and a warm-cold-warm bias in the North Atlantic between 30 • N and 60 • N related to biases in the Gulf Stream extension and the subpolar gyre.  In BACK1, some of these biases are substantially reduced. The most evident improvements can be observed in the Indian Ocean. With the exception of a slight increase in a deep cold bias at around 1,000 to 3,000 m depth and a warm bias south of 55 • S, all biases are reduced. In the Atlantic and Pacific Ocean transects, the Southern Ocean warm bias is reduced, most prominently around 40 • S to 50 • S in the Atlantic and around 60 • S in the Pacific. A clear improvement is also visible in the tropical Atlantic and for the warm-cold-warm bias in the North Atlantic, related to the improved representation of the Gulf Stream extension. Generally, in the upper 1,000 m, the changes due to backscatter are dominated by improvements, while below this depth, the balance between improvements and degradations is more mixed.
Nevertheless, the changes due to BACK1 for the three transects correspond to an overall reduction of RMSE. The RMSEs in the Atlantic, Pacific, and Indian Ocean transects are reduced by 22%, 4%, and 31%, respectively.
For all transects, the changes due to backscatter in the Southern Ocean exhibit a clear dipole pattern of increased temperatures further south and reduced temperatures further north (centered between 50 • S and 60 • S in the Atlantic and Indian Oceans and around 60 • S in the Pacific). This dipole pattern reaches from the surface down to about 3,000-3,500 m and is related to the increased eddy activity in the ACC. Eddies are known to flatten steep isopycnals in the ACC, which is one of the main motivations for the GM parametrization in coarse-resolution, non-eddy-resolving ocean simulations (Gent & McWilliams, 1990;Hallberg, 2013). Since in BACK1, eddies are much more vigorous and numerous, the flattening of the isopycnals is more effective, which leads to a change in the ACC fronts and the temperature advection across the ACC toward the Antarctic continent. Heat is more efficiently exchanged between the north and the south, leading to a cooling in the north and a warming in the south.

Journal of Advances in Modeling
The backscatter simulation shows a strong subsurface heating close to the Antarctic continent, on the continental shelf. Unfortunately, this heating increases some of the local warm biases, thus decreases sea ice thickness. We will discuss this in more detail in sections 3.4 and 3.6.

Impact on Water Mass Properties: Salinity
When looking at salinity biases (Figure 7), REF exhibits a surface and near-surface fresh bias with exceptions close to the Antarctic continent, the high northern latitudes, and parts of the tropics, where the water is too salty near the surface. At depths between 500 and 2,000 m, the waters are generally too salty compared to PHC (Figures 7 and 8).
Southern Ocean biases tend to be improved in BACK1, especially near the surface where backscatter generally leads to saltier waters. Biases are also reduced around 100 m depth close to the Antarctic continent where waters are fresher in BACK1. At depths around 1,000 m to about 2,000 m, the changes in biases are not so clear. There are regions where BACK1 performs worse (e.g., the tropical Atlantic at around 1,000 m depth and parts of the Southern Ocean) and other regions where BACK1 performs better than REF (e.g., east of South America).
As for temperature biases, global RMSEs for salinity biases in BACK1 are overall lower than in REF (Table 1). While the global RMSE reduction with BACK1 at the surface is about 12%, at depth, both simulations have very similar RMSE.
Transects of salinity biases in the Atlantic, Pacific, and Indian Oceans show similar results ( Figure 8) with a general improvement above 1,000 m depth, some mixed results between 1,000 and 2,000 m, and only very minor changes below. The most dominant improvement can be seen in the Indian Ocean, where a strong near-surface fresh bias between 20 • S and 40 • S is reduced by the backscatter. However, below this fresh bias, a bias of too much salt is slightly increased in BACK1.
As expected, biases in absolute terms are larger near the surface and so are the changes due to backscatter. Therefore, bias reduction near the surface is stronger in absolute terms than bias increase between 1,000 and 2,000 m depth. Similar to temperature biases, improvements dominate the upper 1,000 m, while signals are mixed below this depth. While temperature and salinity biases are both locally reduced by more than 50% in the Indian Ocean, bias reductions in salinity are slightly less clear in the entire Atlantic transect than they are for temperature. Nevertheless, the changes due to BACK1 for the three transects correspond to an overall reduction of RMSE. The RMSEs in the Atlantic, Pacific, and Indian Ocean transects are reduced by 13%, 11%, and 28%, respectively.

Impact on Near-Surface Stratification: MLD
The analysis of changes in MLD focuses on the annual mean MLD as described in section 2.3. As expected, the MLD is deepest in the Southern Ocean, along the western boundary currents, in the northern North Atlantic, and close to the Antarctic continent. It is mostly shallow in the tropics (Figure 9). The most promi- In the far northern North Atlantic and Labrador Sea, backscatter tends to reduce maximum MLD with exceptions of deepening in a region between Iceland and Svalbard and south of the southern tip of Greenland. It shows that increased eddy activity interacting with the mean flow can lead to restratification and shoaling as well as deepening due to intensified mean currents and changes in local stratification caused by local and global circulation changes.
Another important change in MLD is the (predominant) shoaling of deep convection around the Antarctic continent. As eddy activity is intensified close to the continent, a likely reason is the process of restratification caused by the eddies and the more stable background stratification due to increased heat transport from lower latitudes. This change in MLD is significant because it provides an explanation as to why we observe subsurface heating on the Antarctic continental shelf due to backscatter (see Figures 5 and 6 and discussion in section 3.2) which tends to increase local warm biases. Increased horizontal heat transport and a subsurface trapping of heat by shallower MLDs in a region where surface heat fluxes generally lead to strong cooling increases surface and subsurface heat storage. This effect is not necessarily desirable, as will become evident in section 3.6.

Impact on Vertical Flow Structure: Overturning Stream Functions
To investigate the impact of backscatter on the large-scale overturning circulation, we analyze the differences between BACK1 and REF for the entire global ocean and the Atlantic (Figure 10).  (Kanzow et al., 2010). Also visible are the shallow tropical overturning cells transporting surface waters away from the equator. The canonical picture of a clockwise upper overturning cell and a deep counterclockwise cell is clearly visible on the global scale as well as in the Atlantic.
The common response of the overturning cells to increased eddy activity in BACK1 is a general intensification of the lower cell and a reduction in the strength of the upper cell. Furthermore, the upper cell tends to become shallower, while the lower cell increases in thickness. In this context it should be noted that 31 years of integration for BACK1 is not long enough for the overturning to equilibrate. However, most of the abovementioned responses have to do with changes in isopycnal surfaces in the Southern Ocean and a balance between upper cell thickness and overturning strength that is steered by the eddy activity and resulting changes in isopycnal slopes in the Southern Ocean. A similar behavior has been observed and discussed in detail by Marshall et al. (2017) in non-eddy-resolving model simulations where the mean eddy impact on isopycnal slopes through an eddy thickness diffusivity was parametrized by a version of the GM eddy parametrization (Gent & McWilliams, 1990). Increasing the eddy diffusivity, that is, the strength of unresolved eddies, led to very similar results as the ones observed here. Furthermore, the shift in the ACC that was also observed for temperature in section 3.2 is partly reflected in the changes of the global overturning circulation in the Southern Ocean.

Impact on Sea Ice
Sea ice is a sensitive component of the climate system which reacts strongly both to changes in atmospheric forcing and to changes in the oceanic heat flux. Since the atmospheric forcing in our two simulations is the same, changes in sea ice predominantly reflect changes in heat transport by the ocean as well as changes in vertical ocean heat fluxes.
Annual mean sea ice thickness in the north tends to increase by up to 25 cm with backscatter, especially in the central Arctic ( Figure 11). Only along the ice edge, we can see some reduction in sea ice thickness in some locations. For Antarctic sea ice, the response is very uniform. We see a general reduction in ice thickness. The amplitude is similar to the changes in the Arctic, with a maximum reduction of about 25 cm. However, since Antarctic sea ice is a factor of 4 to 5 thinner than Arctic sea ice, this reduction is substantial. Moreover, only in very few regions, ice thickness does not change or is slightly increased, for example, between the Riisen-Larsen and Cosmonauts seas.
The reason for the general reduction in Antarctic sea ice thickness has been discussed in previous sections: In section 3.2, we noted an increase in temperatures close to the Antarctic continental shelf, leading to an increased local warm bias with backscatter. In section 3.4, we observed a shoaling of the mixed layer on the shelf caused by an increased eddy restratification, trapping more warm water in the subsurface while simultaneously leading to more horizontal heat transport. This leads to a reduction of sea ice thickness through the increase of heat fluxes via increased eddy activity.
While Northern Hemisphere sea ice is within the range of observed thicknesses and covers an area which is just slightly too large in the North Atlantic (cf. to Figure 7 of Danilov et al., 2017), Southern Ocean sea ice is generally underestimated by FESOM2, especially in the summer months (see Figure 7 of Danilov et al., 2017). A further reduction of Antarctic sea ice through backscatter is therefore not a desirable outcome. This is especially true for the summer months (not shown) where sea ice extent is too low and thickness tends to be underestimated as well. However, a well-calibrated balance between Southern and Northern Hemisphere sea ice with respect to observations is a problem which many models face. Sea ice models are very sensitive to oceanic and atmospheric forcing and encompass a large number of parametrizations which need to be tuned. Sea ice parametrizations often show different responses in the north compared to the south as sea ice conditions are very different in both hemispheres (e.g., Hunke et al., 2010).
Antarctic summer sea ice is reduced and too low in BACK1. Nonetheless, it is still present, indicating that the general water mass structure is preserved. We expect to be able to address the performance degradation with respect to Antarctic sea ice through changes in the backscatter scheme. A reduction of the strength of backscatter over the continental shelf as well as a retuning of model parameters in the presence of backscatter are needed, see sections 4 and 5.

UKE Details
UKE is the prognostic variable which tracks the energy dissipation of the resolved dynamics and the backscatter from unresolved to the resolved flow (see equation (A4)). It controls the backscatter amplitude. The mean UKE is generally positive (i.e., backscatter is active) and highlights regions where the backscatter coefficient is largest (Figure 12). This is especially true in the western boundary currents, the Southern Ocean, on shelves, and to some extent also in a band close to the equator. UKE tends to be quite localized, reflecting that strong dissipation is generally also localized near intense eddies or fronts or close to islands or topographic features.
We observe high UKE values on shelves near Antarctica or along the western boundary currents. For some of these locations, it is desirable that excessively dissipated energy is reinjected (e.g., for the East Australian Current). For some shelf regions and also, to some extent, for the Gulf Stream and Kuroshio, the energy available for backscatter might be too large and may contribute to some of the detrimental effects of backscatter discussed in the previous sections.
In the backscatter scheme, the fraction of dissipated energy which enters the UKE budget is governed by the local Rossby number as suggested by Klöwer et al. (2018) and also used in Juricke et al. (2019). The remaining fraction of dissipated energy is discarded to account for physical dissipation. Klöwer et al. (2018) developed their model for the dissipation fraction in the context of western boundary intensification in a double gyre North Atlantic box setup. Thus, it is adjusted to the pathways of dissipation in turbulent western boundary currents but does not accommodate for effects of coastal shelves; neither Klöwer et al. (2018) nor  included complex bottom topography or transient forcing by variable winds or buoyancy fluxes. Thus, the question how the fraction of physically dissipated energy should be modeled to properly include the effects of transient forcing, shelf turbulence, and underresolved viscous boundary layers requires further study.

Model Stability
One of the practical concerns when implementing backscatter is model stability. Parametrized viscosity is not only required to maintain the scaling properties of turbulence but also necessary to keep the model numerically stable. Since backscatter is antidiffusive, model stability could potentially be an issue.
In our setup, the model time step had to be halved to 10 min when backscatter was switched on to keep the model numerically stable. The reason for this was violations of the CFL criterion with the 20 min time step as the flow velocities with backscatter are substantially increased (see Figure 3). Moreover, model variability is much greater with backscatter switched on. As a consequence, stability can sporadically be violated even with a time step of 10 min. This happened on a few occasions during the 31 years of integration. To reduce the occurrence of such incidents, the time step could be further reduced or the backscatter scheme could be recalibrated to be less efficient in the vicinity of very strong flow structures (see also section 4). However, in the cases when backscatter led to model instabilities with BACK1, it was sufficient to restart the same year in which the instability occurred from single-precision data to overcome the instability. Similarly, a short-term (a few months) reduction of the backscatter coefficient R dis or the time step led to similar results. This suggests that the current setup for BACK1 is at the limit of model stability but still stable almost all of the time.

Sensitivity Studies and Potential Tuning
So far, we have used the backscatter scheme developed by Juricke et al. (2019) essentially without tuning; the reduction of the coefficient R dis relative to the earlier study was done for stability reasons, that is, not related to model quality. Thus, the question of how the model with backscatter responds to changes in parameters, with the idea of retuning the entire ocean model in mind, is very important. While we cannot expect the backscatter to solve every problem, it may be possible to reduce or eliminate some of the biases and model degradations noted in sections 3.2, 3.4, and 3.6.
To this end, we carried out a set of short, 6 year sensitivity studies for which we kept the setup of BACK1 with the exception of changes to only one of the parameters in each sensitivity run. The parameters of interest here are the horizontal Redi (isoneutral) diffusion coefficient K h (locally scaled in relation to a reference resolution of 1 • at the equator), the lead closing parameter h 0 which is part of the parametrization set governing sea ice thermodynamics, that is, thermodynamic growth and melt of sea ice, and the amplitude of the backscatter dissipation through R dis (see equation (A7)). In all cases except the last we tested increased and decreased values, that is, K h = 300 m 2 s −1 and K h = 30, 000 m 2 s −1 instead of K h = 3, 000 m 2 s −1 , h 0 = 1.0 m and h 0 = 0.25 m instead of h 0 = 0.5 m. For R dis , we only tested a decreased value of R dis = 0.1 instead of R dis = 0.5, since further increase (i.e., stronger backscatter) would negatively affect model stability and performance.
The rationale for potentially retuning these coefficients is the following. Increased eddy activity motivates a reduction of the diffusion coefficient, as diffusion at eddy-permitting resolution is partially parametrizing unresolved eddy effects. Since some of these eddy effects were previously included in the larger value for diffusion, they are now at least partly resolved, and therefore, it may not be necessary to account for them in the value of the diffusion coefficients anymore. The lead closing parameter in the sea ice model is a strong tuning parameter and inherently uncertain. With changed ocean heat fluxes, it is generally not uncommon to adjust the sea ice model through the lead closing parameter to achieve better results (e.g., Shi & Lohmann, 2017).
In our sensitivity study, we observed the following. Changes in horizontal diffusion led to a considerable reduction of the warming close to Antarctica and of the Antarctic sea ice decrease, with spatial patterns closely resembling the bias increase due to backscatter. Similarly, a reduction of the backscatter amplitude or changes to the lead closing parameter h 0 also had a clear impact on the Antarctic sea ice: A reduction in R dis led to increased sea ice thickness compared to BACK1. The same holds for an increase in the lead closing parameter. A decrease in the lead closing parameter led to a further decrease in ice thickness. The changes with R dis were more pronounced and coherent across the Antarctic sea ice compared to the Arctic sea ice; changes in the lead closing parameter and horizontal diffusion had similar effects in both hemispheres, that is, a reduced value for h 0 or an increased value for K h led to reduced thickness in the north and the south and vice versa. Since the simulations where only carried out for a few years, they were not long enough to come to a quantitative conclusion, but we feel confident that the trends are robust.
We conclude that a tuning of the backscatter coefficient itself might be necessary (see discussion in section 5). But also tuning of the sea ice model or horizontal diffusion can help to reduce some of the biases. More complex formulations for the lead closing parameter are being investigated as h 0 may depend on wind, ocean, and ice velocities and other factors (see, e.g., summary in Shi & Lohmann, 2017). Such new parametrizations may help to focus on bias reductions for the Antarctic sea ice without affecting the Arctic. Interestingly, increasing the Antarctic sea ice thickness through the lead closing parameter also helped to reduce the subsurface warm bias. This points at another potential cause for this bias, namely, that reduced sea ice leads to more vigorous flow variability which enhances the heat transport toward the coast in BACK1. Such a positive feedback could be partly responsible for the subsurface warm bias on the Antarctic shelf and can therefore be counteracted by changes in the sea ice model.
It should be mentioned that these tuning results do not exclude the possibility of reducing specific biases also through retuning of other parameters. Longer and more extensive sensitivity runs will be necessary to clarify this.

Discussion and Outlook
We analyze a kinetic energy backscatter parametrization in a global ocean model, specifically the unstructured-mesh, finite volume ocean model FESOM2 at eddy-permitting resolution. The backscatter parametrization, introduced by Juricke et al. (2019), is an extension of the work by Jansen et al. (2015). Juricke et al. (2019) tested the scheme in a zonal channel where it considerably improved eddy activity and MKE. They anticipated that the scheme may require considerable adjustment before it would be suitable in a global setup. However, it turns out that the "default scheme" selected by Juricke et al. (2019) performs very well also in global simulations. The only necessary adjustment was a slight increase in the ratio of dissipated versus reinjected energy and a halving of the time step to ensure numerical stability with the resulting velocities and more intense mesoscale turbulence.
The backscatter scheme considerably improves the eddy variability nearly everywhere when judged by comparing SSH variability from AVISO data with model output. Especially the Southern Ocean and parts of the western boundary currents are much more energetic and-with the exception of the Kuroshio-represent observed SSH variability more accurately. Improved eddy variability also reduces biases in mean SSH, temperature, and salinity especially in the North Atlantic, the Indian, and Southern Oceans. Increased biases both in SSH standard deviation and mean temperature and salinity are mostly related to an overintensification of eddy activity close to the coastline for the western boundary currents and around the Antarctic shelf. For the latter, an intensified warm bias leads to reduced Antarctic sea ice, which corresponds to an increase in model bias.
Overall, when looking at globally averaged RMSE, backscatter improves SSH variability and mean, and surface temperature and salinity biases each by about 10%. In the Southern Ocean, reductions in area averaged RMSE for SSH variability are between 30% and 50%, depending on region. The regions where backscatter overenergizes the flow largely correspond to strong currents near continental shelves or the western basin boundaries. There are areas where physical dissipation is strong so that the fraction of energy entering the UKE budget should be tuned down. Since the scheme was optimized in a simplified channel setup, this behavior is expected and consistent with our understanding of the energy cycle; further below, we discuss specific changes to the backscatter scheme to improve its behavior in the problematic areas we identified.
We emphasize, however, that the local biases introduced by backscatter are, at worst, comparable to other biases which routinely occur when running ocean models at approximately this resolution.
Follow-up studies will focus more strongly on the effect of variable resolution in combination with backscatter (see Juricke et al., 2019, for results in a channel) for an uncoupled ocean but also for coupled simulations. As FESOM2 permits local grid refinements, we believe that increased grid resolution in eddy active areas in combination with backscatter can be a powerful tool to improve model performance in regions such as the western boundary currents or the high latitudes (Sein et al., 2016(Sein et al., , 2017. In high latitudes, in particular, the Rossby radius is small and currently not well resolved even in simulations with very high resolution (Wekerle et al., 2017).
Another aspect to be investigated in more detail is the ratio of backscatter versus physical dissipation. In the present study, the energy available for backscatter is a fraction of the energy dissipated by the viscous operator. The fraction is estimated by the local Rossby number following Klöwer et al. (2018). The assumption here is that eddy dissipation in balanced flow is, to a large degree, unphysical and should be scattered back. In unbalanced flows, some amount of energy should be dissipated to represent a down-scale subgrid energy cascade. However, as pointed out in sections 3.6 and 4, the local Rossby number heuristic does not work well in some regions, most crucially on the Antarctic continental shelf and other continental shelves. A significant part of transient variability on the shelves is not related to geostrophic turbulence but is forced. The Rossby radius on shelves is also very small, which is why any inverse energy cascade is difficult to resolve, especially with the resolution of the current setup. This is why backscatter may need to be substantially reduced here, since it relies on an at least partly resolved inverse cascade . While this is generally the case in the open ocean, additional pathways to energy dissipation will need to be accounted for close to the coastlines. Another aspect is that in the presence of topography and along western boundary currents, kinetic energy dissipation and generation are not necessarily collocated in space. Finally, spontaneous wave generation (e.g., Chouksey et al., 2018;Shakespeare & Taylor, 2016) and loss of balance (e.g., Molemaker et al., 2005) are additional mechanisms that can lead to dissipation of energy from balanced flow, although the fraction of this energy dissipation is not well constrained. As a consequence, the parametrization for the fraction of dissipation will be improved in future studies.
Furthermore, the current setup does not combine backscatter with the classical GM parametrization for unresolved eddy effects. The GM parameterization is commonly used on non-eddy-permitting 1 • (or coarser) meshes to simulate the release of potential energy by the missing eddy field. In recent years, however, this parametrization has also been used locally in eddy-permitting simulations for regions where they remain non-eddy permitting (e.g., Hallberg, 2013;Sein et al., 2017). In these studies, the GM scheme is scaled down in regions where eddies are partly resolved to transition between the eddy-mean parametrization and a partly or fully resolved eddy field. Two recent conceptual studies in idealized settings by Bachman (2019) and Jansen et al. (2019) combined-using different approaches-GM and backscatter to make best use of both parametrizations. Such ideas will be further investigated in future studies as well.
Performance wise, the backscatter simulation is about 2.5 times more expensive than the simulation without backscatter. Most of the increase in computing costs originates from the halving of the time step. This increase is especially large on our quasi-regular 1/4 • mesh where the increase in mean velocities with backscatter is substantial, necessitating a smaller time step to match the CFL stability criterion. On other FESOM2 meshes with resolution varying in wider limits (e.g., Sein et al., 2017), the time step will be smaller to begin with, since it has to follow the CFL limit of the smallest triangles in the mesh. There, a time step of 10 min is not unusual even for globally relatively coarse meshes. We need to test whether a further time step reduction will be necessary on such meshes. However, even with such an increase in cost, the benefits are substantial. The model bias in eddy activity as indicated by SSH variability is reduced by more than 50% in large parts of the Southern Ocean, with some biases disappearing entirely. This study of backscatter in a global ocean model holds substantial promise for much more accurate simulations of eddy activity in current climate models.