The Impact of Horizontal Resolution on Energy Transfers in Global Ocean Models

The ocean is a turbulent fluid with processes acting on a variety of spatio-temporal scales. The estimates of energy fluxes between length scales allows us to understand how the mean flow is maintained as well as how mesoscale eddies are formed and dissipated. Here, we quantify the kinetic energy budget in a suite of realistic global ocean models, with varying horizontal resolution and horizontal viscosity. We show that eddy-permitting ocean models have weaker kinetic energy cascades than eddy-resolving models due to discrepancies in the effect of wind forcing, horizontal viscosity, potential to kinetic energy conversion, and nonlinear interactions on the kinetic energy (KE) budget. However, the change in eddy kinetic energy between the eddy-permitting and the eddy-resolving model is not enough to noticeably change the scale where the inverse cascade arrests or the Rhines scale. In addition, we show that the mechanism by which baroclinic flows organise into barotropic flows is weaker at lower resolution, resulting in a more baroclinic flow. Hence, the horizontal resolution impacts the vertical structure of the simulated flow. Our results suggest that the effect of mesoscale eddies can be parameterised by enhancing the potential to kinetic energy conversion, i.e., the horizontal pressure gradients, or enhancing the inverse cascade of kinetic energy.


Introduction
Ocean currents and eddies have systematically less kinetic energy (KE) in low-resolution simulations than high-resolution simulations [1,2].Less energetic currents and eddies can cause biases in e.g., the strength of overturning circulations [3], oceanic heat uptake [4], and mixed-layer depth [5] as well as model drifts [6].To understand what sets the energy levels in ocean models at different horizontal resolutions and what parameterisations may be needed, we must understand how KE is transfered between length scales.A key length scale for the transfer of KE in the ocean is the first baroclinic Rossby radius (hereafter just "Rossby radius"), which is the characteristic length scale at which perturbations in the first baroclinic mode become geostrophically balanced.On scales larger than the Rossby radius, the ocean behaves, to first order, as a quasi two-dimensional turbulent flow due to the strong stratification.A prominent feature of two-dimensional turbulence is the "inverse cascade", where nonlinear interactions transfer KE from smaller to larger spatial scales [7].This is unlike the "forward cascade" in three-dimensional turbulence, where KE is transferred from larger to smaller scales and then dissipated [8].For a two-layer flow, Salmon and Scott, et al. [9,10] showed a simplified view of the transfer of KE as a function of length scale.The baroclinic mode exhibits a forward cascade of potential energy (PE) from large scales to the Rossby radius scale, L R , where there is conversion from PE to KE.Some of the baroclinic KE is then transfered to barotropic KE, often called "barotropisation" [11] and its efficiency depends on on several factors such as stratification and bottom drag [12][13][14][15].Both the barotropic and baroclinic modes exhibit an inverse cascade of KE from the Rossby radius scale to larger scales [10,16], while there is a forward cascade of KE from the Rossby radius toward smaller scales [17].The magnitude of large-scale KE in an ocean model thus depends on the ability of the model to resolve the Rossby radius, typically L R ∼ 10-50 km in midlatitudes, in order to produce the correct transfer of KE and PE between large scales (∼O(10 4 km)) and mesoscale eddies.Non-eddying ocean models, e.g., with 1 • (∼100 km) horizontal resolution, often use an eddy-parameterisation scheme such as [18] (hereafter GM) to account for unresolved eddy processes.However, while such parameterisations imitate tracer fluxes by eddies, they do not include the inverse cascade of KE.As a result, the Gulf Stream, Antarctic Circumpolar Current (ACC), and Kuroshio are generally weak in low-resolution models.Eddy-permitting, ∼0.25 • (∼25 km) horizontal resolution, and eddy-resolving models, ∼1/10 • -1/12 • (∼5-10 km) horizontal resolution, resolve the Rossby radius in midlatitudes and therefore simulate more realistic energy spectra and fluxes than non-eddying models [19,20].However, in eddy-permitting models, the horizontal resolution is close to the Rossby radius in midlatitudes, and thus conversion of PE to KE and the inverse cascade of KE occur near the grid scale where the effects of viscosity often dominates [1,20].Thus, eddy-permitting models generally underestimate the overall KE as well as the magnitude of energy transfers between length scales while eddy-resolving models are more realistic.
As eddy-permitting ocean simulations often struggle to reproduce the strong currents and eddy activity in regions such as the Gulf Stream, recent years have seen much development in parameterising the effect of unresolved scales in eddy-permitting ocean models.One method is to make viscosity scale-aware in order to make the model less viscous while still ensuring numerical stability [21][22][23].Other methods focus on augmenting the weak inverse cascade of KE in eddy-permitting models [19,20,24,25].In addition, horizontal pressure gradients, i.e., the conversion from PE to KE, can be intensified by accounting for unresolved temperature and salinity fluctuations in the nonlinear equation of state [26].A few studies have analysed the spectral properties of KE and PE in global ocean models.Schloesser and Eden [27] studied the North Atlantic and found an inverse cascade of KE and a forward cascade of PE.Arbic, et al. [1] used a very high-resolution model, 1/32 • (∼3 km), and found that both temporal and spatial filtering can introduce an artificial forward cascade of KE at the smallest resolved scales.Chemke, et al. [15] used an ocean state estimate to analyse the barotropic KE budget and understand how barotropisation scales with properties of the flow (e.g., stratification, vertical velocity shear) following the [12] theory.In this paper, we will investigate why varying horizontal resolution changes the amount of KE and if this is true for several regions with strong currents and eddy activity in the world oceans.We will use three simulations with the nucleus for european modelling of the ocean (NEMO) ocean model at non-eddying (1 • ), eddypermitting (1/4 • ), and eddy-resolving resolution (1/12 • ) to diagnose the differences.We present the KE budget in Section 2, and describe the model simulations in Section 3. We will, in Section 4, show that horizontal resolution and the necessary horizontal viscosity exerts a strong control on the KE transfers between length scales, the relative importance of terms in the KE budget, and the partitioning of KE into barotropic and baroclinic modes.We will lastly, in Section 5, summarise our results and discuss possible routes for parameterising unresolved energy transfers.

Kinetic Energy in Physical Space
The 2D KE equation in our ocean model, NEMO ( [28], see Section 3 for further details), can be summarised as where is the horizontal velocity, w is the vertical velocity, p is pressure, and ρ 0 is the reference density.The horizontal gradient operator is denoted ∇ = (∂/∂ x , ∂/∂ y ) and the vertical gradient is ∂/∂ z .The ocean is forced at the surface by wind input, given by the vertical derivatives of zonal and meridional components of wind stress, τ x , τ y .Dissipation at the ocean lower boundary is done using quadratic bottom drag with a coefficient given by Γ = C D u 2 + v 2 + K bg , where K bg is the background turbulent kinetic energy and C D is a non-dimensional parameter.Horizontal and vertical viscosity are set by the respective coefficients A h,m and A v,m .The viscosity operator uses n = 2 for non-eddying simulations and n = 4 eddying simulations, respectively, as is commonly done in ocean simulations [29][30][31].
We identify the terms in Equation (1) as kinetic energy for the horizontal velocity components KE, horizontal nonlinear interactions N, vertical nonlinear interactions N z , wind forcing W, bottom drag D, horizontal pressure work (which is responsible for the conversion between PE and KE) P, horizontal viscosity V, and vertical viscosity V z .We make the approximation that ∂τ x /∂z = τ x /H and ∂τ y /∂z = τ y /H, where H is a reference depth set to either H = 500 m or H = 1000 m depending on the region studied (as explained in Section 3).We calculate the PE to KE conversion using the pressure gradient term above to be consistent with the underlying NEMO calculations (note that it is different from e.g., [17] who used w b , where b is buoyancy).We perform these calculations on results from an idealised double-gyre simulation with parameter settings similar to ORCA025 and compare to the online tendencies from the model, and find only very small differences (see Section 3).When calculating KE budgets, we ignore vertical transfer of momentum, N z , and vertical viscosity, V z .These terms are much smaller than all other terms and do not contribute to the 2D KE budget in real or spectral space.We will also ignore bottom drag, D, when focusing on the upper ocean.

Two-Layer Turbulence
As described in Section 3, we will average the ocean model data onto two layers, ∆z 1 and ∆z 2 , separated at a depth z sep , to partition KE to the barotropic and baroclinic modes.With this approach, the baroclinic mode will encapsulate all baroclinic modes into one.We denote the velocity of the barotropic mode as u bt = (u bt , v bt ) and the velocity for the baroclinic mode as u bc = (u bc , v bc ) where the baroclinic velocities are defined in both the upper layer and lower layer.Similarly to [15], we define the velocity in the two modes as where H tot = ∆z 1 + ∆z 2 is the total ocean depth.Subscripts 1 and 2 denote the upper and lower layers, respectively.The full derivation of the barotropic and baroclinic KE budgets is presented in Appendix A. In this paper, we will mainly focus on the nonlinear interactions and barotropisation in the barotropic and baroclinic KE budgets.These are identified as N bt,bt,bt (barotropic self-interactions) and N bt,bc,bc (barotropisation), N bc,bc,bc (baroclinic self-interactions), and N bc,bc,bt ("catalytic" term).As described in Appendix A, and following [16], we add the "catalytic" baroclinic KE nonlinear interactions, N bc,bc,bt , to the baroclinic self-interaction to get the total baroclinic nonlinear term.In our calculations, the catalytic term has a very similar structure to the baroclinic self-interaction, so adding them enhances the total baroclinic KE flux.We will ignore the "barotropic stirring" term as this is small in the baroclinic KE budget.Scott, et al. [16] argues that the barotropic stirring term in the baroclinic KE budget corresponds to the barotropisation term in the barotropic KE budget.We also note that the baroclinic self-interaction, N bc,bc,bc , is non-zero only when the two layers have different thickness [16].

Kinetic Energy in Spectral Space
In order to study the KE budget as a function of length scales, we do a Fourier transform in the horizontal directions: where k, l are zonal and meridional wavenumbers and curly braces, {}, denotes the Fourier transformed variable.We make similar Fourier transforms of velocity, wind forcing, pressure and their respective gradients to get the time-averaged KE budget in spectral space.We obtain the budget for the KE power spectrum by integrating over wavenumbers where E K = 0.5[{u} * {u} + {v} * {v}], and * denotes complex conjugates.Square brackets, [. . .] = H −1 0 −H . . .dz is a vertical average.We only do a vertical average over the upper ocean, i.e., H = 500 m or H = 1000 m, as defined in Section 3. In the above integral, we assume that the flow is isotropic and proceed to use only the isotropic wavenumber, K = √ k 2 + l 2 .This gives us the budget where Π(K) denotes the integral of a transfer T(k, l), which are given by where u s , v s are surface velocities, u bot , v bot are bottom velocities, [. . .] = H −1 0 −H . . .dz is a vertical average as before and (. . . ) = 1 t 1 −t 0 t 1 t 0 (. . . ) dt is a time average.The time average is taken over the entire simulation length, 1 January 1979 to 31 December 2009.As previously stated, we ignore vertical momentum transfer, N z , and vertical viscosity, V z , as they are relatively small terms.We also ignore bottom drag, D, when studying the upper ocean.Note that we evaluate gradients in (x, y) space before making the Fourier transform, e.g., we calculate u • ∇u, then Fourier transform to get {u • ∇u} and then multiply by {u} * to get the contribution to KE.This is different from e.g., [10], where gradients are evaluated in spectral space from the Fourier transformed velocities.Our method also differs from Aluie, et al. [32], which may be more accurate.
In the case of PE conversion to KE near the Rossby radius wavelength, 2πL R , and a KE flux from 2πL R toward larger scales of O(10 6 m) as shown by [33], we would expect Π N < 0 and Π P > 0 with the extrema between 2πL R and O(10 6 m).Important scales for interpreting the spectra and fluxes of KE are the energy-containing scale, defined as the mean wavenumber weighted by energy, i.e., and the first baroclinic Rossby radius of deformation where ∂z [34].Typically, in ocean mid-latitudes, the energy-containing scale is ∼300 km [34] and the Rossby radius is L R ∼ 10-50 km [34].In our paper, we will calculate the Rossby radius for each grid cell in our eddy-resolving simulation (ORCA0083) in each region, and then present the regional average.As we will show later, K C is approximately the wavenumber at which the KE spectrum, E , as well as Π P and Π N reach their extremes.
We make similar transforms for the baroclinic and barotropic KE and the nonlinear terms.The barotropic and baroclinic KE are where square brackets denote a vertical average over the two layers ∆z 1 and ∆z 2 as before.We integrate as in Equation ( 5) to get power spectra E bt (K) and E bc (K) as well as spectral fluxes.The time-averaged barotropic and baroclinic KE budgets are thus The full derivation of these terms can be found in the Appendix.In this paper, we will focus on barotropic nonlinear interactions, Π bt,bt,bt , baroclinic nonlinear interactions, Π bc,bc,bc + Π bc,bc,bt and barotropisation, Π bt,bc,bc .

Model Data
We use model data from realistic configurations of the 3D NEMO ocean model [28] with three different horizontal resolutions.Data are available as five-day averages, which is high enough temporal resolution to resolve most of the mesoscale eddy activity.The model simulations are the same as used by [35].The only differences between our three different configurations are the horizontal resolution, horizontal viscosity, A h,m , tracer diffusion, and time step [28].All configurations use a NEMO 3.6 ocean model with the LIM2 sea-ice model, and are forced at the surface with ERA-Interim [36] wind and buoyancy forcing for 1979-2009.All three configurations all have 75 vertical z * levels where the volume of a column can vary in time.We ignore this variation in layer thickness, but the associated error is on the order of 1%.Bottom topography is represented with partial steps.All configurations use Laplacian iso-neutral tracer diffusion to represent sub-grid scale tracer mixing.Vertical viscosity, A v,m , is calculated from a turbulent kinetic energy (TKE) closure scheme, with a minimum background value given in Table 1.A quadratic bottom drag is applied in all configurations.The time step and horizontal viscosity are different at different horizontal resolution in order to ensure numerical stability.Similarly, the tracer diffusion is increased at coarser resolution so that the resulting fields are smooth.In ORCA025 and ORCA0083, the order of the horizontal viscosity operator is increased in order to make the dissipation more concentrated to the grid scale [21,29].Eddy fluxes of heat and salt are parameterised in ORCA1 using the GM scheme [18], but not at finer resolution since eddies are, at least partially, resolved in ORCA025 and ORCA0083 and thus eddy tracer fluxes are already explicitly simulated.However, in simulations of high-latitude or shelf processes, the GM scheme could be beneficial even at 1/4 • resolution and finer [2,37].The settings used for the ocean model configurations are summarised in Table 1.
To decompose the KE and spectral fluxes of KE into barotropic and baroclinic modes, we average our global ocean model data onto two vertical levels separated at z sep .This allows us to calculate the KE and spectral fluxes in the barotropic mode, leaving the remainder to be the sum of all baroclinic modes.The baroclinic mode will be dominated by the first baroclinic mode, as most of the KE is in the barotropic or the first baroclinic mode [38].The reason for averaging the velocities onto only two levels is to avoid depths intersected by bathymetry, where many grid cells become invalid causing large errors in the spectral transform.).Black dashed lines indicate interface between "upper" and "bottom" ocean in our analysis as described in Section 2. Secondary peak in velocity is part of a warm core eddy that is often found in this region [39].
The KE budget, which is calculated using five-day averages, is not expected to be exactly balanced.The total imbalance in the KE budgets (Equation ( 6)) presented later is of the order 10 −9 m • s −3 in all regions studied.This is partly because we are not including all terms in our calculations, e.g., vertical viscosity, bottom drag, vertical momentum transfer, etc., although they are relatively small terms.The spectral flux of vertical momentum transfer is on the order of 10 −10 , which is an order of magnitude less than the other terms, and vertical viscosity is even smaller.The imbalance in the KE budget is also due to the interpolation we use to calculate the spectral fluxes, as well as the fact that we are only diagnosing the upper ocean, i.e., down to H = 500 m or H = 1000 m, except when we average onto two layers.We have repeated our calculations on an idealised, flat-bottomed, double-gyre configuration of NEMO with parameters similar to the ORCA025 simulation.Comparing our calculated spectral fluxes based on five-day averages with the exact online tendencies from the model, we found some small differences that were mostly due to interpolation (As an example, in physical space, KE contribution from zonal wind forcing is uτ x z , evaluated at the U-point and then interpolated to the T-point.In spectral space, u and τ x z are interpolated to the T-point, then Fourier transformed, and the KE contribution is {u} * {τ x z }.Hence, the two terms will not be exactly equal.).

Kinetic Energy (KE) in Physical Space
In regions of strong currents and eddy activity such as the Kuroshio, ACC and the Gulf Stream, the eddy-resolving ORCA0083 simulation has higher surface KE than the lower-resolution simulations (Figure 2).For individual snapshots (not shown), ORCA0083 also has a clear eastward extension of both the Kuroshio and Gulf Stream which is not simulated in ORCA025 or ORCA1.Such an extension can also be seen in observations from satellite altimetry (not shown).A more energetic upper ocean is also found when looking at a latitude-vertical cross section of the Kuroshio current along 145 • E (Figure 1), where the jet has a peak zonal velocity of U ∼ 0.4 m • s −1 in ORCA0083, U ∼ 0.15 m • s −1 in ORCA025 and U ∼ 0.1 m • s −1 in ORCA1.Similar differences are found for the Gulf Stream, ACC and the Malvinas.The zonal velocity is higher and the horizontal density gradients are sharper in ORCA0083 than in ORCA025 and ORCA1, partly because higher-resolution models can resolve sharper gradients, but also because there is less horizontal viscosity and diffusivity that smooths horizontal velocity and density gradients.
We average the model data onto two levels, upper ocean (level 1) and deep ocean (level 2), separated at z sep = 500 m in the Kuroshio region as described in Section 3. As also mentioned in Section 3, we have experimented with different values of z sep and found that, while the absolute values of our calculations can change, the main results remain valid.The separation between the upper and deep ocean is shown in Figure 3 to approximately separate high KE, strong stratification, and clear seasonal variability in the upper ocean from low KE, weak stratification and very little variability at depth.We calculate root-mean-square vertical shear and upper-ocean density anomalies as where subscripts 1, 2 indicate the upper and lower layers, < • • • > is an area average and (. . . ) is a time average over 1979-2009.We find that the mean vertical shear, as well as the density anomalies in the upper layer, are greater in ORCA0083 than in the two other simulations in all regions (Kuroshio shown in Figure 4).The ORCA0083 simulation has larger density anomalies in the upper ocean Kuroshio, which indicates more buoyancy in the upper layer, b 1 = gρ 1 /ρ 1 , and thus more available potential energy (APE) in ORCA0083 than ORCA025 and ORCA1, since APE ∼ b 2 /N 2 ∼ (ρ 1 ) 2 /N 2 [40].The ORCA025 and ORCA0083 simulations have relatively similar profiles of N 2 = g/ρ 0 ∂ρ 1 /∂z, while ORCA1 is different (Figure 3), and this is generally true for all regions studied here.Calculating APE as E A = 0.5(b ) 2 /N 2 [40], we find more APE in ORCA0083 than ORCA025 and ORCA1 in the Kuroshio as well as the other regions marked in Figure 2. Hence, even though the surface buoyancy forcing is the same in the three simulations, there is more APE that may be converted to KE in ORCA0083.All three simulations have a strong annual cycle in upper-ocean density anomalies, ρ 1 , likely due to the annual cycle in solar heating and freshwater forcing (Figure 4b,d).The power spectrum of ρ 1 shows a stronger semi-annual cycle in ORCA1 than ORCA025 and ORCA0083, which can also be seen in the timeseries.All simulations also have a pronounced annual cycle in mean vertical shear of velocity, ∆ z U, (Figure 4a,c) caused mostly by an annual cycle in winds that changes upper-ocean KE (Figure 3), although it is progressively weaker with lower horizontal resolution.We also find low-frequency variability that is stronger or as strong as annual variability in all simulations and all regions.We stress that all three simulations are forced with the same wind and buoyancy forcing and have the same vertical resolution.The low-frequency variability that impacts the Kuroshio may be the Pacific Decadal Oscillation (PDO) [41].Low-frequency variability has also been reported in observations of the Gulf Stream [42], and ocean-only modelling experiments of the ACC [43], which may explain the high spectral power at long time scales for these regions (not shown).The ORCA1 simulation generally has an order of magnitude weaker variability on timescales longer than one month.We also find ORCA025 to have a less distinct annual cycle in the timeseries of ∆ z U (Figure 4a), due to the fact that annual variability has less spectral power than low-frequency variability in that simulation.However, in ORCA0083, the annual and inter-annual variabilities are of similar power and the annual cycle is more distinct in the timeseries.The differences between the three simulations in annual and low-frequency variability is likely due to the fact that the temporal inverse cascade of KE, i.e., the transfer of KE from monthly to inter-annual timescales, is weaker with lower resolution [44,45].

Spectral Fluxes in the Upper Ocean
We are interested in the impact of horizontal resolution on the KE balance of oceanic flows and will therefore focus on the KE in four regions of high eddy activity and one region with low eddy activity.The five regions are marked in Figure 2: the ACC in the Pacific sector, the Malvinas region in the South Atlantic, the Kuroshio extension, the Gulf Stream extension and, for comparison, a region in the East Pacific where eddy activity is low.The ACC region is the same as studied by [10], which allows for a comparison with their results.We will use only the upper ocean, defined in Section 3, and investigate the Fourier transformed KE budget defined in Equation (6).
We calculate the spectral KE density, E (K), from Equation ( 5) for all regions (Figure 5).We find that the eddy-resolving ORCA0083 simulation has more KE than the other two simulations at all wavenumbers.We also find that the energy-containing scale is shifted from ∼320 km to ∼280 km from ORCA025 to ORCA0083 in the Kuroshio, and observe similar shifts for the other regions.However, if the energy-containing scale is defined by the peak in E , as in [15,46], then we find no change in the energy-containing scale between ORCA025 and ORCA0083.The non-eddying ORCA1 run has at least an order of magnitude less KE than the other two simulations, which is true whether the eddy-induced velocities from the GM scheme are added to the Eulerian velocities or not.We also find that ORCA025 and ORCA0083 have a clear peak near the 300 km scale, while ORCA1 has no peak.Hence, the KE spectra in ORCA025 and ORCA0083 peak near their respective energy-containing scales.As we will show later in this section, KE fluxes converge to the energy-containing scale of ∼200-300 km in all regions of eddy activity studied here.The KE spectra (Figure 5) are somewhat steeper than the E ∼ k −3 slope predicted in the 2D turbulence inertial range when stirring and dissipation are isolated to large and smaller scales respectively.The slope is steeper for ORCA025 than ORCA0083.As we will show later, the spectral fluxes due to wind forcing and PE to KE conversion occur over a range of wavenumbers, and thus there is no clear inertial range where a k −3 slope could be expected.The two eddying configurations, ORCA025 and ORCA0083, clearly show negative spectral fluxes for nonlinear interactions, Π N < 0, corresponding to KE transfer from smaller to larger scales, i.e., an inverse cascade, with a peak near the energy-containing scale of ∼300 km (Figure 6).An inverse cascade corresponds to either eddies merging into larger eddies, or eddies rectifying a jet such as the Kuroshio, as found in previous results for global ocean models, e.g., [1,27].We find positive values of Π P with a peak near the energy-containing scale, hence opposite the inverse cascade of KE (Figure 6).
The scale at which Π N and Π P peak is larger than the Rossby radius of deformation, L R , shown in Figure 5 for ORCA0083.
However, the fastest growing mode can be up to a factor ∼ 2 larger than the Rossby radius wavelength, 2πL d [19,33].Taking L R ∼ 20 km, which is the average Rossby radius in our ACC region, gives ∼ 2 • 2π • 20 km ≈ 250 km which is close to the peak of the inverse cascade, Π N , and PE to KE conversion, Π P .The positive lobe of PE to KE conversion corresponds to adding KE to scales smaller than the energy-containing scale and removing KE from the larger scales.
The energy-containing scale of ∼300 km is significantly larger than the first baroclinic Rossby radius (Figure 6), and larger than the horizontal resolution of ORCA1.However, since ORCA1 does not resolve or parameterise the process by which KE fluxes to larger scales, the high KE at the energy-containing scale can not be simulated by ORCA1.
The inverse cascade of KE, Π N , is stronger in ORCA0083, weaker in ORCA025, and nearly zero in ORCA1 in all regions, which is in agreement with [1,20] who found that increasing viscosity or using coarser horizontal resolution in an idealised model dissipates KE that would otherwise be transferred to larger scales.Here, it is confirmed to also be the case in state-of-the-art global ocean models with realistic bathymetry and forcing.Similarly to the inverse cascade, and in agreement with [47], the conversion of PE to KE, Π P , attain larger values in ORCA0083, than in ORCA025 or ORCA1.More conversion of PE to KE in ORCA0083 than ORCA025 and ORCA1 implies stronger horizontal pressure gradients and thus a stronger geostrophic flow, which can be seen in a cross section of the Kuroshio (Figure 1) and the other regions studied (not shown), where the zonal flow and the meridional density gradients are stronger at higher resolution.As previously noted, stronger horizontal pressure gradients can be explained by the higher horizontal resolution as well as the lower iso-neutral tracer diffusion.In combination, this means that, in ORCA025, buoyancy anomalies will be subjected to more tracer diffusion and at somewhat larger scales than in ORCA0083.Hence, even though all three simulations have the same surface buoyancy forcing, higher horizontal resolution makes the ocean model more efficient in converting PE to KE.
We find that the spectral transfer for wind forcing, T W = ∂Π W /∂K, is the same at scales O(10 6 m), i.e., Π W has the same slope, for all three simulations, but they differ near the energy-containing eddy scales.Wind forcing is a source of energy at larger scales but, for ORCA025 and ORCA0083, wind forcing is also a sink of KE at near the energy-containing scale.This is due to the fact that zonal wind stress accelerates a zonal jet, but induces a torque that opposes the vorticity of the eddies so that eddies lose KE [48].Eddies on scales near the energy-containing scale in ORCA0083 lose KE to wind forcing and viscosity in roughly equal portions, while, in ORCA025, viscosity dominates over wind (Figure 6).Hence, the manner by which eddies lose KE can be very different depending on the horizontal resolution of the model.It is not clear if this impacts the shape or the lifespan of eddies.Overall, the net KE input from wind is largest in ORCA1 for all regions, except in the East Pacific where there is almost no eddy activity and thus the wind forcing is identical in all three simulations.We note that the KE budgets in Figure 6 do not close perfectly.As stated in Section 3, vertical viscosity and vertical momentum transfer are omitted due to being small but still introduce small errors when attempting to close the budget.However the largest errors in the budget are due to interpolation errors and fluxes of KE on the lateral boundaries of the regions studied.In all simulations, there is a non-zero flux of KE in and out of all the regions studied.We have repeated our calculations using the output of an idealised double-gyre simulation with settings similar to ORCA025, and found that the spectral fluxes sum to nearly zero on the boundaries.
In geostrophic turbulence, the inverse cascade is predicted to terminate near the Rhines scale where most of the KE is contained within large-scale Rossby waves.The Rhines scale is defined as L β = 2 < u > /β [49], where < u > is the root-mean-square zonal velocity and β = ∂ f /∂y.The Rhines scale scales with the eddy KE, and is thus smaller in ORCA025 than in ORCA0083, implying that the inverse cascade can act over a wider range of scales in ORCA0083 than in ORCA025.The Rhines scale has also been found to coincide with the peak of the KE spectrum [49,50], i.e., the peak of the KE spectrum shifts toward smaller scales at lower resolution [19] However, we find neither a shift in the peak of the KE spectrum (Figure 5) nor a shift in the scale at which the inverse cascade arrests (Figure 6) from ORCA025 to ORCA0083 in any of the regions studied.We note that the change in Rhines scale from ORCA025 to ORCA0083 in all regions is only ∼5-10 km, which is too small to be resolved in the simulations.Therefore while the eddy resolving model has more eddy KE than the eddy permitting simulation, this change in eddy KE is not large enough to drastically impact the Rhines scale or the arrest of the inverse cascade.This explains why the peak of the KE spectrum and the inverse cascade do not shift from ORCA025 to ORCA0083.We find two orders of magnitude difference in KE between ORCA1 and ORCA025, and thus a much larger difference in Rhines scale between ORCA1 and ORCA025 than between ORCA025 and ORCA0083.However, ORCA1 has no clear peak of the KE spectrum and no inverse cascade, so we can not relate the change in Rhines scale between ORCA1 and ORCA025 to a shift in scale where the inverse cascade arrests.
The non-eddying simulation, ORCA1, generally shows almost no spectral flux from nonlinear interactions or PE to KE conversion except at the largest scales.We include the eddy-induced velocities in our analysis, which are responsible for some conversion of PE to KE by flattening isopycnals to generate eddy velocities.This could be the positive Π P seen at large scales in ORCA1 (Figure 6).It is clear, however, that the eddy-induced velocities do not result in an inverse KE cascade, which may be the main reason for the low KE.

Baroclinic and Barotropic Kinetic Energy (KE)
As described in Section 3, we average the model data onto two vertical layers, representing the upper and deep ocean, to limit the effect of bathymetry.We calculate the barotropic spectral flux, Π bt,bt,bt , the baroclinic spectral flux (with the added catalytic flux), Π bc,bc,bc + Π bc,bc,bt , and the spectral flux due to barotropisation, Π bt,bc,bc [16] as described in Appendix A. We find a distinct inverse cascade in both the barotropic and baroclinic modes (Figure 7), similar to [16], in all regions, although it is very weak in the East Pacific.The two terms Π bc,bc,bt and Π bc,bc,bc have very similar structure, thus adding the catalytic flux to the baroclinic self-interaction enchances the total baroclinic KE flux.In both ORCA025 and ORCA0083, barotropisation injects KE over the 150-600 km range, and the barotropic inverse cascade transfers energy from the 150-300 km range to the 300-600 km range (Figure 8).This implies a KE sink over this range, which can partly be explained by bottom drag (Figure 8), and wind forcing (Figure 6).The scales over which barotropisation acts, as well as the scale at which the barotropic KE spectrum peaks, are unchanged from ORCA025 to ORCA0083.It is difficult to conclude whether the abovementioned scales are insensitive to horizontal resolution from studying just two eddying simulations.

Wavelength [km]
P bc,bc,bc + P bc,bc,bt ORCA1 P bc,bc,bc + P bc,bc,bt ORCA025 P bc,bc,bc + P bc,bc,bt ORCA0083 P bt,bt,bt ORCA1 P bt,bc,bc ORCA1 P bt,bt,bt ORCA025 P bt,bc,bc ORCA025 P bt,bt,bt ORCA0083 P bt,bc,bc ORCA0083 Figure 7. Spectral flux for baroclinic nonlinear interactions , Π bc,bc,bc + Π bc,bc,bt (dashed lines, left panels), barotropic nonlinear interactions, Π bt,bt,bt (solid lines, right panels) and barotropisation, Π bt,bc,bc (dotted lines, right panels) for ORCA1 (blue), ORCA025 (red), and ORCA0083 (green).Note the different scales for the East Pacific region.As in Figure 6, note also that spectral fluxes have units [m • s −3 ] since they are the integrals of spectral transfers, which have units of KE tendency, [m 2 • s −3 ].The baroclinic inverse cascade does not approach zero at the smallest wavenumber, implying that the energy-containing scales drive baroclinic flow on scales larger than the regions studied here.The results show that the baroclinic inverse cascade is about ∼7-10 times larger than barotropic inverse cascade in both ORCA025 and ORCA0083.This is likely due to the strong bottom friction [51] which, in idealised models, results in a baroclinic inverse cascade much larger than the barotropic inverse cascade [16].Furthermore, we find that the East Pacific region has inverse cascades of barotropic and baroclinic KE, but of at least one order of magnitude smaller than in the other regions due to the low eddy activity.As is the case for the KE cascades in the upper ocean, the ORCA025 simulation clearly has weaker inverse cascades of barotropic and baroclinic KE than ORCA0083.It is clear that the ORCA1 simulation has virtually no KE cascades, as shown for the upper ocean (Figure 6).We stress that with our two-layer approach, all baroclinic modes are encapsulated into one.The inverse cascade of baroclinic KE may thus be a residual of an inverse cascade in the first baroclinic mode and forward cascades in the higher modes [10,11], and the partitioning of KE between baroclinic modes can differ between simulations.For realistic ocean stratifications, KE is efficiently transferred from higher baroclinic modes to the first baroclinic mode, and subsequently to the barotropic mode [13,52].It is therefore likely that the barotropisation mostly represents the transfer between the first baroclinic mode and the barotropic mode.
We find clear differences in the strength of barotropisation between our simulations, with stronger barotropisation in ORCA0083 than ORCA025 and very small barotropisation in ORCA1.The differences in barotropisation could be due to differences in bottom friction, stratification, or barotropic eddy KE, which have been show to impact the barotropisation [11][12][13][14][15][16]46].Comparing ORCA025 to ORCA0083 in each region, we find that the mean stratification, N 2 , in the upper ocean (Figure 3) is similar in both simulations, and can thus not explain the differences in barotropisation between the simulations.Differences in bottom drag in the barotropic mode (shown for the Kuroshio in Figure 8), can explain some of the stronger barotropisation at higher resolution [15].Furthermore, some of the differences in barotropisation can also be explained by the fact that horizontal viscosity dominates over barotropisation and the barotropic inverse cascade in ORCA025, but not in ORCA0083.While the relative effects of bottom drag and horizontal viscosity on the barotropisation are unclear from our results, we note that [15] found bottom drag to have a low effect on barotropisation.Hence, we speculate that horizontal viscosity is the dominant factor here.To isolate the impact of bottom drag, one could run multiple ORCA025 and ORCA0083 simulations with varying bottom drag coefficient, C D , but this is beyond the scope of this paper.
For each of the regions, we calculate the area-averaged barotropic and baroclinic KE, and then the fraction of the total KE that is in the barotropic mode (Figure 9).For the western boundary currents, the ORCA0083 simulation agrees with observations [38], where it was reported that the Kuroshio and Gulf Stream have about 40% and 50% of the total KE in the barotropic mode, respectively.We find that the ACC and Malvinas have ∼60% and ∼80% of KE in the barotropic mode, respectively, but it is difficult to compare to observations as there is very poor observational coverage for these regions in [38], and we are unaware of any other observations of this.We find that in all regions except the East Pacific, the flow is more barotropic in ORCA0083 than in ORCA025.The weaker barotropisation in ORCA025 compared to ORCA0083 (Figure 7) suggests that the lack of barotropic energy in ORCA025 is due to less barotropisation, which can be explained by baroclinic eddies that are becoming more barotropic being dissipated faster in ORCA025.In the ACC, Malvinas and Gulf Stream, the ORCA1 simulation is even less barotropic than the other simulations, while, in the Kuroshio, it is strongly barotropic.Thus, horizontal resolution can have an impact on the vertical structure of the simulated flow by changing the partition of barotropic and baroclinic KE, which can be seen as a bias in e.g., the vertical shear.

Discussion and Conclusions
We have presented the energy budget for state-of-the-art global ocean simulations of varying resolution, which allowed us to study the impact of horizontal resolution on the flow in models with realistic bathymetry, wind forcing, and bottom drag.To our knowledge, a spectral KE budget for state-of-the-art ocean models has not been presented before.The eddy-permitting simulation, ORCA025, showed a weaker inverse cascade of KE than the eddy-resolving simulation, ORCA0083 (Figure 6), which was caused partly by the lower horizontal resolution, and partly by the necessary higher viscosity.There were also clear differences in the effect of wind forcing and conversion of PE to KE between ORCA025 and ORCA0083 (Figure 6).As a result, in ORCA025 currents were too weak (Figure 1), the midlatitude jets did not extend far enough eastward (Figure 2), and viscosity had a dominating effect on the KE budget.
We found no shift in the scale at which the KE spectrum peaks from ORCA025 to ORCA0083, and also that the inverse cascade acted over the same range of scales in both simulations.However, in idealised experiments, Jansen, et al. [19] found that the peak of the KE spectrum shifts to larger scales with higher resolution, which implies that the inverse cascade arrests at larger scales as well.We explain this disagreement by the fact that the Rhines scale, which predicts where the inverse cascade arrests and the KE spectrum peaks, changed by only ∼5-10 km from ORCA025 to ORCA0083.Hence, the shift in Rhines scale was too small to be resolved by either the ORCA025 or the ORCA0083 simulations.The small differences in the Rhines stem from the small differences in eddy KE between the eddy permitting and eddy resolving simulations.We similarly found no shift in the scales over which the barotropisation and the barotropic inverse cascade act from ORCA025 to ORCA0083.Studying how these macroturbulent scales change with horizontal resolution is interesting and would require a larger suite of eddying simulations than the two (ORCA025 and ORCA0083) presented here.
Several studies [12,13,15] have investigated the sensitivity of barotropisation to e.g., stratification, bottom drag and Coriolis parameter.However, to our knowledge, none of these studies investigated the sensitivity of barotropisation to horizontal resolution or horizontal viscosity.We found that differences in stratification, N 2 (Figure 3), could not explain the large differences in barotropisation.Comparing ORCA025 to ORCA0083, we found that the stronger bottom drag at higher resolution could explain at least part of the stronger barotropisation, although bottom drag was noted to have a low effect on barotropisation in [15].The differences in horizontal viscosity and resolution may explain some of the differences in barotropisation as well, as horizontal viscosity has a much larger impact on the barotropic KE budget in ORCA025 than in ORCA0083 (Figure 8).Further sensitivity studies are needed to fully understand the impact of viscosity and bottom drag on barotropisation in realistic global ocean models.
Partitioning the KE into barotropic and a baroclinic mode, we found that the ORCA025 simulation had a more baroclinic flow than ORCA0083, since part of the KE that should be transfered from the baroclinic mode to the barotropic mode is dissipated by the higher viscosity, although bottom drag was lower.This suggests that increasing the horizontal resolution or reducing horizontal viscosity (while maintaining numerical stability) can reduce biases in the simulated vertical structure, although the vertical resolution is important as well.When the flow is too baroclinic, it will impact the relative importance of terms such as buoyancy forcing over Reynolds stresses in the eddy energy budget.Bottom drag, which primarily affects the barotropic mode [16] and is the primary way of dissipation in the ocean [51], was underestimated in ORCA025 due to the lower barotropic KE, thus changing the relative importance of drag and viscosity on KE dissipation.Furthermore, the vertical structure of the flow, to some extent, set the level at which the effective horizontal mixing is largest in regions like the ACC or western boundary currents [53][54][55].Thus, simulating a realistic vertical structure is of vital importance for the horizontal mixing of tracers in the model.
Our findings show that the some key processes in mesoscale eddies can be different in eddy-resolving and eddy-permitting ocean models.The eddy-resolving ORCA0083 had a significantly larger sink of oceanic eddy KE due to wind forcing than the eddy-permitting ORCA025.An eddy KE sink implies a transfer of momentum from the ocean to the atmosphere, which in turn implies that, as horizontal resolution of the ocean in coupled climate models increases, the momentum budget for the atmospheric boundary layer will change.The shapes of eddies may also be different, as high viscosity could prevent eddies from becoming horizontally anisotropic [56] and the lack of barotropisation leads to too much vertical shear.
For ORCA025 and ORCA0083, our spectral KE budget results agree with the energy cycle in physical space diagnosed by [57] who showed a route of energy where mean KE and mean PE can generate eddy PE, which is then converted to baroclinic eddy KE and lastly barotropic eddy KE.This is also mostly in line with [58].In the Lorenz energy cycle by [58], winds and dissipation are the main source and sink for eddy KE, respectively.However, our results show that conversion from eddy PE is the main source of eddy KE, while both dissipation and winds act as sinks, similarly to [48,59].In the framework of a Lorenz energy cycle, our results show that all parts of the transfer from mean PE to eddy PE and subsequently to eddy KE are weaker at eddy-permitting resolution and that their underestimation is balanced by higher rates of dissipation of KE and diffusion of PE.
We suggest that we can improve the eddy-permitting simulation by decreasing viscosity, parameterising PE to KE conversion or parameterising sub-grid transfers of KE.Lowering the horizontal viscosity parameter, A h,m , can lead to some improvements of the simulated currents but also leads to more numerical noise [60].A more promising way to improve the simulations is to use a flow-dependent viscosity [21,22] where the dissipation is set by the deformation rate or vorticity gradients.Such methods generally lead to improvements by having overall less dissipation while still suppressing numerical noise [23,61], which could result in a stronger inverse cascade of KE and more barotropisation.
It is possible to improve eddy-permitting simulations by amplifying the PE to KE conversion, i.e., strengthening the horizontal pressure gradients (Figure 6).Such an approach was taken by [26], who used stochastic processes to represent unresolved density fluctuations in a coarse-resolution model.Brankart, et al. [26] found that the parameterisation resulted in more ageostrophic motions, more intense vertical velocities, and improvements in many regions with intense eddy activity.Using a coupled climate model, Ref. [62] found that the parameterisation by [26] reduced some of the biases in low-resolution climate models as well as increased the interannual variability of the overturning circulation.Horizontal pressure gradients can also be strengthened by strengthening horizontal gradients of sea surface height, which would require increasing the divergence of the barotropic velocities, e.g., using a stochastic parameterisation to represent unresolved divergent motions.
Despite the fact that the KE input from wind forcing at large scales was very similar in all simulations, the ORCA0083 simulation had more KE at all wavenumbers due to a stronger inverse cascade (Figure 5).A way to improve eddy-permitting simulations like ORCA025 is to parameterise or enhance the inverse cascade of KE, which can be done in a few different ways.The "backscattering" parameterisation by [19] aims to calculate the KE lost due to viscosity and re-inject it at larger scales in order to parameterise the inverse cascade.The "α model" by [63,64] uses filtered velocities where Rossby radius is shifted to larger scales, thus making it more clearly resolved.Work has also been done using the "Rivlin-Ericksen" stress [20,25] where subgrid eddy stresses are represented by a stress tensor that can strengthen velocity gradients and thus amplify eddies.Implementing the above parameterisations in ocean models [20,64,65] has shown that they enhance the inverse cascade and allow an eddy-permitting model to have a KE spectrum more similar to that of an eddy-resolving model.
Future work should involve calculating KE budgets for ocean models of even higher horizontal resolution that include sub-mesoscales and internal waves, e.g., 1/48 • -MITgcm simulation [66], to diagnose to what extent the ORCA0083 simulation underestimates KE of ocean currents and eddies.We would also propose to perform sensitivity experiments with state-of-the-art global ocean models to disentangle the roles of horizontal resolution, bottom drag and horizontal viscosity in controlling barotropisation.This would help in understanding what sets the scales of the inverse cascade and the peak of the KE spectra in our simulations, which could be studied using a larger suite of simulations of varying horizontal resolution.In the future, we also hope to investigate the KE budgets in global ocean models that utilise some of the parameterisations mentioned above, e.g., scale-aware viscosity or augmented energy transfers, to see how well they may reduce biases in the inverse cascade of KE or barotropisation for eddy-permitting models.Work is currently ongoing to implement and test new parameterisations of subgrid energy transfers for state-of-the-art ocean models.
where vertical nonlinear interactions and vertical viscosity are ignored as there is no vertical component in the barotropic mode.We multiply by baroclinic velocities, u bc , v bc , in the upper and lower layers and integrate vertically to get the baroclinic KE budget as  The first remaining nonlinear interaction in the baroclinic KE budget is "barotropic stirring", N bc,bt,bc , which represents baroclinic energy advected by the barotropic mode and corresponds to the barotropisation in the barotropic KE budget [16].The second term is the "catalytic" term [16], N bc,bc,bt , which involves the barotropic mode and redistributes KE within the baroclinic mode.The last term involves only the baroclinic mode and is the baroclinic self-interaction, N bc,bc,bc .where the barotropic forcing, F bt , includes the surface pressure gradient term and wind forcing, while barotropic dissipation, D bt , includes bottom drag and horizontal viscosity from Equation (A1).Baroclinic forcing, F bc , includes the hydrostatic pressure gradient term, wind forcing and vertical momentum transfer, while dissipation, D bc , includes bottom drag, horizontal viscosity and vertical viscosity from Equation (A2).Following [16], we add the baroclinic self-interaction to the "catalytic" term to get the total baroclinic nonlinear interactions, N bc,bc,bc + N bc,bc,bt .We find that both are very similar in sign and magnitude.In our calculations, barotropic stirring is the smallest term in the baroclinic KE budget, and we therefore ignore it.In spectral space, the barotropic and baroclinic KE budgets are found in the same way as the full KE budget (Equation ( 6)): where, as in the main text, curly braces {} denote Fourier transformed variables, and square brackets, [. . .] denotes a vertical average.

Figure 1 .
Figure 1.Cross section of time-averaged (1979-2009) zonal velocity in the Kuroshio from ORCA1 (a), ORCA025 (b) and ORCA0083 (c).Black solid lines show isopycnals (labels in kg • m −3).Black dashed lines indicate interface between "upper" and "bottom" ocean in our analysis as described in Section 2. Secondary peak in velocity is part of a warm core eddy that is often found in this region[39].

Figure 2 .Figure 3 .
Figure 2. Time-averaged (1979-2009) kinetic energy at the surface in ORCA1, ORCA025, and ORCA0083.Overall, the currents are more diffuse at lower resolution.Note that the colour scale is nonlinear.

Figure 4 .
Figure 4. (a) timeseries of root-mean-square (RMS) vertical shear in ORCA1 (blue), ORCA025 (red), and ORCA0083 (green); (b) RMS upper-ocean density anomalies for the Kuroshio region.Horizontal dashed lines indicate time average; (c,d) power spectral density in units m 2 • s −2 and (kg • m −3 ) 4 of the two time series in the top and middle figures respectively.Horizontal axis is 2π/t, where t is time in days.Vertical dotted line shows annual frequency.

Figure 5 .
Figure5.Kinetic energy spectrum, E K , as a function of wavenumber, for ORCA1 (blue), ORCA025 (red), and ORCA0083 (green).Vertical lines show energy-containing scale, K C .Black vertical line shows first baroclinic Rossby radius, L R , in ORCA0083, although it is nearly identical in all simulations.Numbers show the result of a fit to determine p for E ∼ k p between the energy-containing scale and the highest wavenumber.

Figure 6 .
Figure 6.Spectral flux for nonlinear interactions, Π N , and pressure gradient, Π P , (Equation (1)) for ORCA1 (blue), ORCA025 (red), and ORCA0083 (green).In the left figure, advection, Π N , is solid, and pressure gradient, Π P , is dashed.In the right figure, wind forcing, Π W , is solid, viscosity, Π V , is dotted, and bottom friction.Note the different vertical scale for the East Pacific region.Note also that spectral fluxes have units [m • s −3 ] since they are the integrals of spectral transfers, which have units of kinetic energy (KE) tendency, [m 2 • s −3 ].

Table 1 .
Parameters used for the ocean model simulations used in this study.The vertical coordinate z * is explained in the text.Note that the simulations have different forms of horizontal viscosity.