Parameterizing Mesoscale Eddy Buoyancy Transport Over Sloping Topography

Most of the ocean's kinetic energy is contained within the mesoscale eddy field. Models that do not resolve these eddies tend to parameterize their impacts such that the parameterized transport of buoyancy and tracers reduces the large‐scale available potential energy and spreads tracers. However, the parameterizations used in the ocean components of current generation Earth System Models rely on an assumption of a flat ocean floor even though observations and high‐resolution modeling show that eddy transport is sensitive to the potential vorticity gradients associated with a sloping seafloor. We show that buoyancy transport coefficient diagnosed from idealized eddy‐resolving simulations is indeed reduced over both prograde and retrograde bottom slopes (topographic wave propagation along or against the mean flow, respectively) and that the reduction can be skilfully captured by a mixing length parameterization by introducing the topographic Rhines scale as a length scale. This modified “GM” parameterization enhances the strength of thermal wind currents over the slopes in coarse‐resolution, non‐eddying, simulations. We find that in realistic global coarse‐resolution simulations the impact of topography is most pronounced at high latitudes, enhancing the mean flow strength and reducing temperature and salinity biases. Reducing the buoyancy transport coefficient further with a mean‐flow dependent eddy efficiency factor, has notable effects also at lower latitudes and leads to reduction of global mean biases.


Introduction
At present, the ocean components of most global climate models are used at resolutions that require parameterizing transport by the oceanic mesoscale (Fox-Kemper et al., 2019).Although coupled simulations with eddying ocean fields are slowly emerging (Chang et al., 2020), mesoscale eddy parameterizations are still likely part of ocean models for another decade.The wave-turbulence duality of mesoscale eddy dynamics can cause very rich transport behavior, involving intermittency, non-local transport by coherent vortices and even upgradient fluxes that energize the mean flow (e.g., Chen et al., 2014;Liu et al., 2023;Yankovsky et al., 2022).Nonetheless, most present-day parameterizations do not include such effects, but have their origins in the works of Gent and Mcwilliams (1990), Gent et al. (1995), andRedi (1982), tackling eddy-induced advection and tracer mixing, respectively.The "GM" advection is specified by an overturning streamfunction which itself is cast in terms of a horizontally down-gradient and vertically up-gradient buoyancy transport-resulting in a reduction of available potential energy."Redi" diffusion, by contrast, mixes both active and passive tracers down-gradient along isopycnals (Gent, 2011).
Much of the current research focuses on how to prescribe flow-dependent eddy transport coefficients, or eddy vertical structure of eddy buoyancy transport certainly appears to be different from that of isopycnal Redi diffusion, including the diffusion of potential vorticity (PV) (see e.g., Abernathey et al., 2013;Bachman et al., 2020;K. S. Smith & Marshall, 2009).Less is known about how GM and Redi diffusivities relate laterally, but since both are sensitive to for example, eddy energy levels, we expect them be spatially correlated.In practice, the diffusivities, whether applied to the GM or Redi scheme, are most often parameterized following mixing length theory.This means that they are constructed from the product of some eddy velocity scale and a mixing length scale.Some work has gone into estimating the eddy velocity scale by implementing a prognostic equation for eddy energy (Bachman, 2019;Eden & Greatbatch, 2008;Jansen et al., 2019;Mak et al., 2018;Marshall et al., 2012), but this is still very much an active and incomplete field of research.So, in the latest iteration of the Climate Model Inter Comparison Project (CMIP; https://explore.es-doc.org/cmip6/models/), the GFDL-CM4.0model (Adcroft et al., 2019) was, to our knowledge, the only one that used a prognostic eddy energy approach to estimate the eddy velocity scale.
The study by Visbeck et al. (1997) therefore continues to influence the practical use of the mixing length approach.Drawing on earlier works by Green (1970) and Stone (1972), they proposed that the velocity scale for the GM diffusivity be based on the product of the growth rate of baroclinic instability in the linearized Eady model (Eady, 1949) and some length scale.Assuming that the mixing length is also set by the same scale, the diffusivity will then scale as the Eady growth rate and the square of the length scale.Visbeck et al. (1997) associated the mixing length with the "width of the baroclinic zone" which they defined as "the width of the region where the local growth rate exceeds 10% of the maximum growth rate of the field."The concept, however, is hard to define in any but the most idealized model geometries, and length scales therefore need to be formed from more rigorous dynamical arguments.
As proposed by Stone (1972), one obvious candidate for length scale is the internal deformation radius, the approximate scale of the fastest unstable growth in the Eady model.Solid observational evidence for the relevance of this length scale has been presented by Stammer (1997) and Eden (2007).However, other relevant scales arise if dynamics beyond the Eady framework are accounted for, most notably bottom friction and potential vorticity (PV) gradients.Jansen et al. (2015), for example, examined the role of bottom friction and the planetary vorticity gradient in a two-layer flat-bottom channel model.They found that bottom friction primarily influences the vertical distribution of eddy energy and that the mixing length in most of their simulations is set by the Rhines scale, that is, the transition scale between nonlinear and linear PV dynamics on the flat-bottom planetary beta plane (Rhines, 1977).More generally, Jansen et al. (2015) found that in order to cover various dynamical regimes, the smaller of several candidate length scales should be chosen, and that inclusion of the Rhines scale amongst these scales is important.In fact, the observational studies of both Stammer (1997) and Eden (2007) specifically pointed to a minimum of the internal deformation radius and the Rhines scale as a best fit for eddy length scales over much of the world ocean.
These principles remain the standard in state-of-the-art models, although development has occurred in later years.As mentioned above, there has been extensive focus on developing prognostic equations for eddy energy.Considerable efforts have also gone into studying effects of horizontal eddy anisotropy (R. D. Smith & Gent, 2004) and the suppression of mixing across strong mean flows (Ferrari & Nikurashin, 2010;Klocker et al., 2012, and references therein).It's worth noting, however, that most of the development up until recently has been guided by observed dynamics in low and mid latitudes.Current parameterizations thus lack any treatment of two aspects that are potentially of huge importance in high latitude oceans, namely the presence of sea ice and the potential vorticity gradients imposed by sloping bottom topography.A sea ice cover can effectively have the same influence as bottom friction on growth of baroclinic instability as well as on dissipation of existing mesoscale and sub-mesoscale eddies (Meneghello et al., 2021).This topic, however, will be left out from the present study.We will instead focus on the dynamical impacts of bottom slopes, that is, continental slopes and mid-ocean ridge systems, whose imprints can be easily seen in observations of both mean currents and mesoscale energy fields, especially at high northern latitudes (Koszalka et al., 2011;Nøst & Isachsen, 2003;Trodahl & Isachsen, 2018).Imprints of topographic PV gradients can also be seen at lower latitudes, for example, in drifter and float paths (Fratantoni, 2001;LaCasce, 2000).
Sloping bottom topography can suppress growth rate and reduce length scales of baroclinic instability (e.g., Blumsack & Gierasch, 1972;Brink, 2012;Isachsen, 2011;Mechoso, 1980) as well as impact finite-amplitude eddy fields (e.g., Bretherton & Haidvogel, 1976;Lacasce & Brink, 2000; K. Stewart et al., 2015;Vallis & Maltrud, 1993;Wang & Stewart, 2018).To this end, new topography-aware parameterizations have started to emerge, both for eddy-induced advection and isopycnal mixing.In particular, Wang and Stewart (2020) and Wei et al. (2022) tested different scaling relations for the GM diffusivity in high-resolution model simulations of flows over idealized continental slopes in re-entrant channels.The two works examined eddy characteristics and fluxes across retrograde and prograde mean currents, respectively, meaning currents that are in the opposite and same direction as topographic waves.Both studies diagnosed the eddy energy from the high-resolution fields and used this to examine traditional mixing length formulations, trying out various choices for mixing length.In addition, they tested the "GEOMETRIC" formulation of Marshall et al. (2012) in which diffusivities are instead constructed from eddy energy and an eddy decorrelation time scale which is set equal to the inverse of the Eady growth rate.In general, the two formulations performed similarly, suggesting that a good knowledge of the eddy energy field is key.However, importantly, both studies also found that empirical prefactors that depend on the topographic slope are needed to reproduce very weak eddy buoyancy fluxes across sloping bottom topography.Wei and Wang (2021) carried on from Wang and Stewart (2020), but focused on the along-isopycnal tracer (Redi) diffusivity in the same channel model-in retrograde flows only.The authors constructed a parameterized Redi diffusivity from (the square root of) the diagnosed eddy kinetic energy and the internal deformation radius, again finding that the actual diffusivity over the slope was suppressed below the original scale estimate.However, instead of testing a set of empirical slope-dependent prefactors, as done by Wang and Stewart (2020) and Wei et al. (2022), this study picked up from Ferrari and Nikurashin (2010) and demonstrated that mean-flow suppression could explain the observed reduction in cross-slope fluxes near the surface, whereas eddy velocity anisotropy contributed to the reduction close to the bottom.
In other words, both sets of studies (see also Brink, 2012Brink, , 2016;;Hetland, 2017) concluded that eddy diffusivities over sloping bottoms are poorly reproduced by traditional open-ocean scaling choices for eddy velocity and eddy length (or decorrelation time), and that additional dynamical impacts of the bottom topography must be brought in.Topographically-induced velocity anisotropy is one obvious factor which could impact both the effective eddy velocity (its orientation relative to the tracer gradient) and effective mixing length or time scale.In addition, mean-flow suppression, caused by eddies propagating relative to the mean flow, may also be reflected in a reduced effective mixing length, as suggested by Ferrari and Nikurashin (2010, their Equation 13).But such interpretations have so far only been applied to Redi diffusion-now also over continental slopes (Wei & Wang, 2021).Whether similar dynamics lie behind the various empirically-fitted suppression factors in the studies of buoyancy diffusion over continental slopes is yet an open question.
The present study will primarily focus on eddy buoyancy transport and thus on GM diffusivities.It is inspired by and builds directly on the results obtained by Wang and Stewart (2020) and Wei et al. (2022), and, as they did, we thus limit the scope to depth-averaged diffusivities.However, as noted, the above works examined prograde and retrograde flows separately and also constructed diffusivities from eddy energy levels diagnosed from very idealized but high-resolution fields.So here we aim to (a) study fluxes and diffusivities over both types of flow situations under one and the same framework, (b) examine how far one can get with parameterizations that do not rely on diagnosing the actual eddy energy field and, finally, (c) expand by assessing impacts both in an idealized setting and in a realistic global ocean simulation.
In the process, we revisit the question of what is the relevant eddy mixing length over continental slopes.The starting point will be the internal deformation radius since this remains a relevant parameter in the Eady problem.In addition, we also consider the topographic Rhines scale, that is, the scale that marks the transition between a linear topographic Rossby wave (rather than planetary Rossby wave) regime and turbulent PV dynamics.The above-mentioned idealized channel studies give conflicting evidence about the relevance of this scale.We are nevertheless inspired by the findings of Stammer (1997), Eden (2007), and Jansen et al. (2015) and bring up this approach here again.To home in on what actually goes on over the slopes, we will also diagnose the eddy velocity anisotropy and, in addition, the phase relationship between velocity and tracer perturbations.This second diagnostic gives additional information about dynamics not reflected in mere scale estimates, at least estimates of eddy velocity.Essentially, no matter how strong the Root Mean Square eddy velocity is, if velocity perturbations are in quadrature with buoyancy perturbations a zero transport results.The analysis done here will indeed show that most of the topographic suppression is reflected in a degraded phase relationship and that velocity anisotropy takes on a secondary role.
The paper is structured as follows: In Section 2 we introduce the modeling tools and various diagnostics and parameterizations used.In Section 3 we begin by diagnosing eddy fields from high-resolution channel simulations that contain both prograde and retrograde flows at the same time.We then see how far mixing-length and GEOMETRIC parameterizations can take us in reproducing the diagnosed depth-averaged GM diffusivity-with and without accounting for effect of anisotropy and phase relations between eddy velocity and tracer perturbations.At the end of this section we examine the impact of a topographically-aware parameterization in a coarseresolution version of the channel model.In Section 4 we finally employ the new parameterization in realistic global ocean simulation.We then take a critical look into some of our parameterization choices and their interpretation in Section 5 before summarizing our findings in Section 6.

Model Setup
We use the Bergen Layered Ocean Model (BLOM), the ocean component of the Norwegian Earth System Model (NorESM; Seland et al., 2020), in an idealized channel configuration as well as in a realistic global setup (both configurations are published in Nummelin (2023b)).BLOM uses 51 isopycnal levels (potential density referenced to 2,000 dbar) with a 2-level bulk mixed layer at the surface.In order to diagnose the various quantities used in the study, we interpolate the outputs locally (in time and space) to height coordinates, except for quantities that are specifically calculated from isopycnal output (see below), in which case we interpolate the outputs to a new density grid so that the bulk mixed layer is properly accounted for.
The channel setup is re-entrant in the zonal (x) direction.The domain is 416 km long (zonally) and 1,024 km (y max ) wide (meridionally).At both sides of the channel there are continental slopes of given width (W) centered at 150 km (Y C ) from the domain edges, stretching 2,000 m (D S ) in the vertical from the shelf break at 250 m depth (D B ) to the bottom of the slope at 2,250 m depth.These parameters then define the bathymetry (H) across the channel (along the y-coordinate): e. in the central basin.
In addition, to trigger instabilities we add 2D random noise with standard deviation of 10 m to the bottom topography.
The model is initialized from rest with constant salinity and a horizontally homogeneous temperature profile.The temperature, which here determines density alone, has a maximum at the surface and decays exponentially toward the bottom.We place the channel in the northern hemisphere, using a constant Coriolis parameter, and then force the flow with a constant westward wind stress.The surface mixed layer is kept shallow by parameterization of submesoscale mixed layer eddies (Fox-Kemper et al., 2008) that counter the vertical mixing induced by the constant wind forcing.See Table 1 for further parameter settings.
We first run the channel model at 2 km horizontal resolution, which is eddy-resolving over the deep central basin and over the slopes (see Table 2 for deformation radius) but only eddy-permitting over the shallow shelves.To investigate the effects of the two bottom slopes on eddy transport and, specifically, on eddy diffusivity, we vary the initial stratification and the width of the continental slope, that is, the slope angle.The various experiments are laid out in Table 2.All simulations are spun-up for 10 years to a semi-equilibrium in which domain averaged eddy kinetic energy is close to constant, and the model fields are then diagnosed over an additional 5-year period (so between years 11-15).We then test and compare various forms of parameterized eddy buoyancy fluxes at noneddying 32 km resolution in the same idealized channel.These are also run for 15 years, with the last 5 years being diagnosed.With the focus on the 5-year means, we are assuming a slowly-varying eddy field and attempting to parameterize its time-mean impact.However, we note that geostrophic turbulence is known to be intermittent with implications for variability in eddy transport and mixing (Busecke & Abernathey, 2019;Huneke et al., 2019;Ong et al., 2023;G. Zhang et al., 2023).Also our simulations suggest that over the northern prograde slope the jet there undergoes periods of alternating strong and weak eddy activity on a timescale of several months (not shown).How the time-mean view of the eddy transport presented here incorporates such variability is left for future studies.
Finally, the impact of the most skillful parameterization is assessed in realistic global simulations.These are nominal 1°resolution global forced ocean-ice experiments which follow the Ocean Model Intercomparison Project, OMIP-II protocol (Tsujino et al., 2020).In these simulations, the mean grid size north of 62°N and south of 64.5°S is approximately 32 km, similar to the coarse resolution channel.We compare simulation with an existing eddy parameterization, which does not include any effects of bottom topography, to simulations with parameterizations that feel the bottom topography through the topographic beta parameter (Figures S1a-S1c in Supporting Information S1; see Section 2.2 for further definitions).Each simulation is 110 year long (2 cycles of 55 long repeat cycle), and we diagnose the results using the last 30 years.At this point there is still a long-term drift in the model (as seen in all models following the OMIP-II protocol; Tsujino et al., 2020), but the general circulation has stabilized.

Diagnostics and Parameterizations
The key parameter of interest is the buoyancy diffusivity, and in this study, we focus exclusively on the depthaveraged diffusivity.We leave the development of depth-varying parameterizations for future studies.A fruitful way forward for this may be to develop a flow-dependent structure function that distributes the depth-averaged diffusivity vertically (see e.g.Bachman et al., 2020;Wei & Wang, 2021).
In the idealized zonal channel simulations, where buoyancy is given by temperature, the cross-channel (i.e., meridional) buoyancy diffusivity is diagnosed from where H is bottom depth, v is meridional velocity, and T is meridional temperature.Angle brackets indicate a zonal (along-channel) mean and primes indicate deviations from such mean.So v′ and T′ are the across-channel velocity and temperature perturbations from the zonal mean.In Equation 1, the flux gradient relation is evaluated at each level, before depth averaging.For analysis of the channel simulations, we also average K diag over time.
In what follows, we make frequent use of depth-averaged variables, which we note with (⋅) .Thus, our eddy kinetic energy density is defined in terms of depth-averaged velocities as and an velocity anisotropy factor is defined as Parameterizing the diffusivity starts with a scale estimate.Two approaches are currently in use, the traditional mixing length and the GEOMETRIC approach.In the first approach, we write .The mixing length may intuitively be thought of as the size of eddies themselves.We examine this possibility below, estimating L from the shape of velocity spectra.Several possibilities exist for this (see e.g., Eden, 2007), but here we chose where v(k) is the Fourier component of the depth-averaged cross-channel velocity at wavenumber k.The expression can be thought of as a kinetic energy-weighted mean wavelength under the spectrum.
In the second approach, the energy-based diffusivity estimate of the GEOMETRIC framework (Mak et al., 2018;Marshall et al., 2012) is constructed as where σ E is the Eady growth rate and E is the total eddy energy.The Eady growth rate is where f is the Coriolis parameter and Ri is the geostrophic Richardson number: Here is the squared buoyancy frequency (g is gravitational acceleration, ρ is density while ρ 0 is a reference density, so that b = gρ/ρ 0 is buoyancy) and is the magnitude of the thermal wind shear.As said, E is the total eddy energy, that is, the sum of the EKE and EPE (eddy potential energy).The former is diagnosed using Equation 2 above and the latter from where η′ is the height of an isopycnal surface above its mean level.The sum is taken over n density surfaces and, as in all of the above, the prime marks deviations from the zonal mean.Note, finally, that the mixing length diffusivity becomes identical to the GEOMETRIC diffusivity if the mixing length in the former is equal to the "Eady scale," EKE 1/ 2 σ 1 E , and if only EKE is used in the definition of the latter (Wei et al., 2022).

Journal of Advances in Modeling Earth Systems
10.1029/2023MS003806 NUMMELIN AND ISACHSEN Since the work on prognostic eddy energy budgets is still a topic of active research, we set out here to parameterize both the eddy velocity and eddy length scale from coarse-resolution variables.Thus, following Visbeck et al. (1997), we write which, moving from proportionality to equality, gives where a 1 is some proportionality constant.Two parameterizations for the eddy length scale are then assessed, namely the WKB-approximation to the internal Rossby deformation radius, and the parameterized version of the topographic Rhines scale, where β T = (|f|/H)|∇H| is the topographic beta parameter.Here we have assumed V par = σ E L T (Eden & Greatbatch, 2008) and introduce a constant tuning factor a T which may, for example, reflect the resolution of the bathymetric data set used. Figure S2 in Supporting Information S1 shows that in the high-resolution channel model the parameterized L T (using a T = 0.1) correspond reasonably well with L T estimate based on diagnosed EKE.Note that Equation 16 and the resulting velocity (Equation 13) and diffusivity (Equation 14) formulations are the same as suggested by Held and Larichev (1996) for potential vorticity diffusivity in beta-plane turbulence, except for their beta being the planetary beta.Here parameterized velocity and length scales are always chosen consistently that is, the parameterized diffusivities will depend on the Eady growth rate and the squared length scale of choice.
Finally, topographic impacts on velocity anisotropy and the phase relationship between velocity and buoyancy perturbations are also diagnosed from the high-resolution channel simulation.Anisotropy is calculated from Equation 3 while the phase relation is assessed from the cosine of the angle between the real and imaginary parts of the cross spectrum: where Ĉo v′,T′) and Qu v′,T′) are the real and imaginary parts of the cross spectrum, respectively (the cospectrum and quadrature spectrum).For analysis, we average θ across all wavenumbers (k) and over time before calculating the cosine.

Equilibrated Flow Field and Eddy Fluxes
Our setup (see Section 2.1) is very similar to the setup in the series of papers by Wang and Stewart (2018), Wang and Stewart (2020), Wei andWang (2021), andWei et al. (2022), except that we now have continental slopes on both sides of the channel.The forcing is also slightly different as we employ a westward wind stress which, unlike in the previous studies, is kept constant across the channel.The mean ocean state, however, is very similar.Since the channel is in the northern hemisphere, the westward wind stress sets up a northward surface Ekman transport.Thus, Ekman divergence in the south and convergence in the north results in a time-mean sea surface tilt which is in geostrophic balance with a westward mean flow, as shown in the two upper panels of Figure 1.The Ekmandriven overturning circulation in the y-z plane lifts up isopycnals in the south so that they slope with the bathymetry there.Conversely, downwelling in the north sets up isopycnals that slope against the topography.
Despite the simple wind forcing, the total baroclinic velocity field is rather complex (Figure 1a, black dashed lines).In the north there is a strong westward jet over the slope.This jet has a significant thermal wind shear but nonetheless extends all the way to the bottom.Over the southern slope the westward flow is weaker and much more surface-trapped.Lower layers here are almost motionless, so the depth-averaged westward flow takes on a minimum over the slope (Figure 1b).Instead there is a broad and nearly barotropic westward current which has its maximum strength immediately off the seaward side of the continental slope.The north-south asymmetry is clearly not only a result of the stratification being weaker in the south than in the north.Thus, net impacts of mesoscale eddy fluxes must be taken into account.At the most basic level, the tilted isopycnals in both regions are baroclinically unstable, creating an eddy field whose residual mass transport will tend to counter the Ekman-driven overturning circulation.However, because mesoscale eddies also transport momentum, the mean flow field reflects, in part, the integrated effects of eddy momentum and buoyancy fluxes.
Their combined effects can be studied in the Transformed Eulerian Mean version of the zonally-averaged zonal momentum equation: where τ x is the zonal wind stress and is the (quasi-geostrophic) Eliassen-Palm (E-P) flux.It consists of a meridional eddy flux of negative u-momentum and an eddy form stress (this term arises after thickness-weighting).In Equation 18we have neglected small terms describing the transport of zonal mean momentum by the meridional mean flow as well as vertical flux of momentum (see Wang & Stewart, 2018).Note, however, that the eddy form stress term, which is connected to lateral buoyancy transport under the small-slope approximation, may be thought of as a vertical momentum flux.Finally, the Coriolis term contains the residual meridional velocity v*, that is, the equivalent mass transport velocity which accounts for both the Eulerian-mean flow and the mass transport by eddy correlations.
The E-P flux from the high resolution Experiment 3 (Table 2) is shown as arrows in the top panel of Figure 1 and its horizontal component (i.e., 〈v′u′〉) is shown with background shading.So what we see is the direction at which the eddy field transports the westward momentum originally provided by the wind.In general, both in the south and in the north, the downward eddy momentum flux is suppressed over the slopes, in agreement with earlier studies which indicate that baroclinic instability of suppressed over continental slopes.Our estimate of the depth-averaged cross-channel buoyancy diffusivity reflects this signature by being reduced by about two orders of magnitude over the continental slopes (lower panel).What these simulations show, as also seen in the simulations of Wang and Stewart (2018) and Manucharyan and Isachsen (2019), is that eddy motions instead bring zonal momentum laterally across the slopes near the surface and dump it where the ocean bottom flattens off toward the deep basin.There, over the relatively flat bottom, baroclinic instability kicks in to bring the momentum down to the solid ground below.
As lateral eddy momentum fluxes are also clearly important in this and previous simulations, optimal parameterizations will likely need to be build up around down-gradient PV fluxes (see e.g., Wang & Stewart, 2018).However, it is also reasonable to expect that any framework which is successful at reproducing the order-ofmagnitude drop in buoyancy diffusivities seen in Figure 1 will also improve the ocean state in coarse-grained models.So we keep this focus here.Hence, on our way toward a practical parameterization of a GM diffusivity over continental slopes, we begin by examining the length scales and velocity scales associated with the mesoscale eddy field.This approach is motivated by the mixing length argument (Prandtl, 1925), relating diffusivity to an eddy velocity scale and a length scale.However, we will also compare this approach with the energy-based GEOMETRIC framework (Mak et al., 2018;Marshall et al., 2012).

Eddy Length and Velocity Scales
Estimates of eddy length and velocity scales are shown in Figure 2. The length scale is estimated from Equation 6, that is, by calculating a spectral-weighted mean wavelength associated with north-south velocity perturbations.When normalized by its mean value across the channel the length scale shows a near-universal shape across the various model runs (upper left panel).There is a broad maximum over the mid-basin before length scales drop over the continental slopes on both sides.There is, however, a consistent local maximum over mid-slope on the northern (prograde) side, coinciding with the maximum in mean zonal velocity (Figure 1).Scales then flatten out or even increase over the shelf regions.As with other diagnostics below, we will largely ignore shelf values from the discussion due to the model grid not fully resolving the deformation radius there and due to the proximity to model walls.For the eddy velocity scale we show the square root of depth-averaged EKE.When normalized with the across-channel average (upper right panel), the eddy velocity scale in all runs is reduced over the southern slope, save for a slight increase over the upper parts of the slope.In stark contrast, the northern slope is dominated by a large maximum, also that one centered over the upper parts of the slope.The eddy velocity then drops off and flattens out over both shelf regions.
It would seem that forming a diffusivity from the product of these diagnosed length and velocity scales may reproduce the observed reduction over the southern retrograde slope (Figure 1), at least qualitatively.But it should also be clear that this procedure would produce a diffusivity maximum over the northern slope-for which there is absolutely no indication in the model fields.We will return to this issue below but first examine possible scaling approximations to the observed length and velocity scales.
We start by comparing the diagnosed L and V with the classical Stone (1972) prediction.So the diagnosed length scale is normalized by the internal deformation radius L R (Equation 15) and the diagnosed velocity scale is normalized by the product of the Eady growth rate (Equation 8) and the deformation radius, so Leaving out any constant prefactors here, we see that both length scales and velocity scales are well represented by Stone-type scaling in the mid-basin (middle panels).The normalized length scales then drop slightly over the lower parts of both slopes, indicating that the deformation radius overestimates scales there somewhat.Finally, there is a dramatic rise in normalized scales over the upper parts of both slopes as the deformation radius drops toward the shallow shelves.As with length scales, the normalized velocities drop over the lower parts of the slopes before rising again over the upper parts.The normalization brings the EKE peak over the upper parts of the slope down to values similar to those seen over the mid-basin, as the EKE peak there coincides with the region of stronger thermal shear.
Finally, following the suggestion by Eden (2007) and Jansen et al. (2015), we normalize by selecting a smooth minimum of length scales: where L T is the topographic Rhines scale (Equation 16).The results are similar over the central basin since the deformation radius is the smaller of the two scales there (the Rhines scale blows up).But now both normalized length and velocity scales peak over the slopes where the Rhines scale becomes the smaller of the two-and is quite clearly too small to explain the observed fields.As such, consideration of the topographic Rhines scale does not seem to bring any improvement in skill in predicting eddy length scales and velocity scales over the continental slopes.
But before rejecting this scaling choice it is worth noting again that the construction of a diffusivity from the original (non-normalized) length and velocity scale estimates would obviously result in a diffusivity maximum over the central northern slope.Such a maximum is in no way suggested from Figure 1.What may be missing from the story here is a consideration of how eddy velocity anisotropy and the velocity-temperature phase relationship may act to bring diffusivities down over the slopes.So we turn to this issue next.

Anisotropy and Phase Relationship
Figure 3 shows the eddy velocity anisotropy A (Equation 3) and the cosine of the phase angle between real and imaginary parts of the v′ and T′ crossspectra (Equation 17).As expected, the eddy velocity field is close to being isotropic in the middle of the basin (upper panel).Values there are around 0.6, implying that cross-channel velocity fluctuations v′ are in fact slightly larger than along-channel fluctuations u′.The eddy fluctuations then become much more anisotropic toward the continental slopes, with A values over the upper parts of the slope close to 0.1 (0.2) in the north (south).This implies that v′ is about 70% (50%) smaller than u′ in the north (south).A notable exception is a peak over the center of the northern slope where v′ is about 50% larger than u′.We have also tested other measures of anisotropy, such as the velocity based measure used by K. Stewart et al. (2015) that takes rotational aspects into account, and the results are similar to those shown here.
The general behavior of increased anisotropy over the slopes, with |v′| < |u′|, will work to reduce the scale-based diffusivity there.But the variations in A from mid-basin values are not great and the mid-slope peak (where |v′| > |u′|) would actually increase the estimates there.So we conclude from this that velocity anisotropy alone can not explain the consistent drop in diffusivity by two orders of magnitude over the slopes seen in Figure 1.The phase relation, however, is able to explain the observed order-of-magnitude drop over the slopes, as the v′ and T′ fields are close to 90°out of phase there (middle panel).Importantly, the low phase agreement over the northern slope largely cancels the local peak in anisotropy.
The lower panel in Figure 3 shows the product of A and cos(θ) (blue lines), an indication of their combined effect.The total suppression is dominated by the information carried in the phase relationship, and velocity anisotropy primarily plays a role near the edges of the two slopes.The suppression over the slopes amounts to more than an order of magnitude, so it is an effect which clearly needs to be parameterized.The slope-dependent prefactors which previous studies have needed to invoke to explain buoyancy diffusion in similar channel simulations are, in effect, attempts at such parameterization (Brink, 2012(Brink, , 2016;;Hetland, 2017; Wang & Stewart, 2020;Wei et al., 2022).However, at this point we temporarily detour from those earlier studies and instead take as a starting point an expression which bears some resemblance to the final form of the meanflow suppression factor proposed by Ferrari and Nikurashin (2010).Thus, we construct an "eddy efficiency" factor as Here, U bc is the large-scale baroclinic flow speed obtained after subtracting the depth-averaged velocity, V is the eddy velocity scale and a 2 is an additional scaling factor which we here take to be constant.The expression does not have a rigorous basis but a simple intuitive interpretation.U bc is directly related to the thermal wind shear and, hence, to the underlying energy source of baroclinic instability (e.g., Sutyrin et al., 2021).Qualitatively, if U bc is large and the flow is baroclinically unstable, one would expect V to be relatively large, giving E eff ∼ O(1) (unless a 2 is very large).But if V remains small despite large U bc , some dynamical constraints must be reducing the efficiency of baroclinic energy conversion, implying E eff ≪ 1.The above interpretation hinges on the parameterized V being an adequate approximation of the actual eddy velocity scale.
We evaluate Equation 21 at each depth but then take the mean over the water column.The large-scale baroclinic flow U bc is extracted directly from the resolved (and zonally-averaged) velocity field, while the eddy velocity is parameterized from Equation 13.The lower panel of Figure 3 shows the resulting efficiency factor, using L min as length scale.The prefactor a 2 has been manually tuned to match the mid-basin values of A ⋅ cos(θ), but it is clear that using L = L T can produce a suppression over the continental slope which is in qualitative agreement with A ⋅ cos(θ) over both slopes for a range of different simulations.Allowing for another prefactor in front of the whole expression (effectively our a 1 parameter) would enable a good quantitative match both over the mid-basin and the slope regions.In contrast, the comparison clearly shows that using L R as length scale does not reproduce the needed behavior over the slopes.We note that several tests with using the thermal wind instead of U bc and with evaluating Equation 21 with depth averaged-quantities (instead of taking the mean of a depth dependent expression) all produce similar results.Here we chose to use U bc due to the ease of implementation at coarse resolution.

Parameterized Diffusivity
Given the above results, we proceed to examine parameterizations of the diagnosed buoyancy diffusivity.The aim is to capture the order-of-magnitude reduction in diffusivities from the mid-basin to the slope regions.The results are shown in Figure 4 where we distinguish between partial parameterizations (panels a-c) and full parameterizations (panels d-f).The partial parameterizations include extensive information about the mesoscale field itself (such as EKE and L S ) which would not be directly available in a coarse-resolution model, whereas the full parameterizations only use information of the zonally and temporally-averaged background buoyancy field, Coriolis frequency, and topographic beta, and are therefore suitable for direct implementation in any existing coarseresolution model.The one exception is an estimate (VII) which combines diagnosed EKE with a parameterized length scale.Panels a and d are from one single simulation, showing both the actual depth-averaged diffusivity diagnosed (black line) and the various approximations (distinguished by Roman numerals and color).Panels b-c and e-f then show statistics over both slope regions collected over the whole range of simulations.
A first thing to notice from the partial parameterizations is that the mixing length (I) and GEOMETRIC (II) approaches behave nearly identically.This suggests that (a) EKE and EPE are proportional to each other, as found in the simulations of Wei et al. (2022), and (b) that our diagnosed eddy scale, L S , reflects the "Eady scale" Jansen et al., 2015;Kong & Jansen, 2021;Larichev & Held, 1995).As also noted by Wang and Stewart (2020), both approaches thus give reduced diffusivities over the southern retrograde slope (squared Pearson correlation coefficient r 2 is 0.49 and 0.66, for I and II, respectively; note that r = cov(x, y)/(σ x σ y )).But the reduction is still underestimated by up to one order of magnitude and reflected as a large relative error.More importantly, over the prograde slope in the north, both approaches result in a serious qualitative mismatch, as the high EKE and EPE levels there (EKE seen in Figure 2; EPE not shown) produce a non-existing diffusivity peak over mid-slope.Although σ E also peaks over the northern slope (not shown), this is not enough to pull down K GEOM there.We note that the peak in K GEOM over the southern shelf might be spurious as the deformation radius is not well resolved there, which is why we have left the shelves out of this analysis.
The observed discrepancies, particularly the qualitative mismatch over the northern slope, confirms that scaling arguments alone are unable to reproduce the diagnosed diffusivities-even with knowledge of eddy energy levels and eddy sizes.It is worth noting that this is in line with previous studies (e.g., Wang & Stewart, 2020;Wei et al., 2022) who found that such scaling estimates needed to be multiplied by slope-dependent prefactors to align with diagnosed diffusivities.Here, we instead explicitly examine the role of velocity anisotropy and the velocitytemperature perturbation phase relationship.Accounting for the diagnosed velocity anisotropy, so that ̅̅̅̅̅̅̅̅̅ ̅ EKE √ will be replaced with v′ (III), improves the mixing length estimate slightly but not nearly enough.Multiplying the two estimates by A cos θ, however, largely removes the diffusivity peak in the north and even produces a clear suppression over the slope-for both the mixing length and GEOMETRIC estimate (IV and V).The values are still higher than the observed diffusivity (mean absolute relative error stays above 100%) but r 2 increases to 0.71 and 0.77 for (IV) and (V), respectively.Over the retrograde slope in the south the match is even closer, with r 2 reaching 0.92 and 0.85 for (IV) and (V), respectively.
Guided by the observed agreement between the mixing length and GEOMETRIC estimates above, we focus on the former approach when examining how well full parameterizations can do.So we assume that a diffusivity can be written as the Eady growth rate times the square of a length scale.Including our efficiency factor, the effective diffusivity becomes  2) whereas the right panels summarize the statistical comparison between the diagnosed and parameterized diffusivities across all experiments (b-c and e-f; relative error is defined as |(K par K diag )|/K par , where K par is one of the parameterization (I)-(IX); r 2 is based on linear regression using 200 points across all cases that is repeated 5,000 times).Boxes and whiskers come in pairs, with the one on the left (right) corresponding to the southern (northern) slope.Linear regressions are done over the slope regions only (gray shading; similar to Figure 1).In panels a and d, the estimates I-IX are scaled by constant a 1 optimized to match the mid-basin diffusivity, and in addition we use a T = 0.5 (for VII-IX), and a 2 = 2 (for IX).The a 1 values are 0.1, 1.17, 0.09, 0.31, 3.77, 0.25, 0.66, 0.25, and 0.25 for estimates I-IX, respectively.
where K 0 = σ E L 2 is the scaling estimate of diffusivity before considering the efficiency factor and where, as discussed above, we have a choice to make for the length scale.For the constant parameters, we chose a 1 = 0.25 so that mid-basin values of VI closely match with the diagnosed diffusivities across all cases, and we use the same a 1 for (VIII) and (IX), whereas for (VII) we use a 1 = 0.66.Coefficients a 2 and a T are then tuned manually such that the correlation between the diagnosed diffusivities and estimates (VIII) and (IX) over the slopes (Figure 4f) is maximized across all cases.
We start by looking at K 0 first.Using the traditional Stone (1972) expression where the length scale is taken to be the internal deformation radius everywhere, seriously overestimates diffusivities over both slope regions (VI).The estimate, in fact, bears some resemblance with both the mixing length and GEOMETRIC estimates based on diagnosed eddy quantities (I and II), but with even larger discrepancies over the slope regions.
Switching temporarily to a hybrid estimate, with a diagnosed velocity scale but a mixing length set to be the smooth minimum of the deformation radius and the parameterized topographic Rhines scale (VII; a T = 0.5), improves the skill considerably, giving in suppressed diffusivities over both slopes.We note that over the slopes, where L T is selected, the diffusivity estimate is essentially a topographic version of the Rhines-based estimate tested by Jansen et al. (2015).They too found that this scaling (with planetary rather than topographic beta) reproduced diagnosed diffusivities well in an two-layer channel model.Interestingly, in our simulations the skill of this hybrid parameterization is better than both estimates using full eddy quantities (I and II).Our interpretation is that our parameterized L T is a better estimate of the eddy mixing length than the diagnosed eddy size L S (which, from the close correspondence between I and II must reflect the Eady scale L E ).
Finally, using the full parameterization with Eady growth rate and the minimum length scale squared (VIII) leads to further improvements over the slopes where the parameterized diffusivities now drop by nearly two orders of magnitude, approaching the behavior of the diagnosed diffusivities.Multiplying this estimate with the parameterized efficiency factor (IX) improves the match somewhat over both slopes, especially in terms of the mean absolute relative error, which drops from 215% (70%) to 116% (53%) while r 2 increases from 0.73 (0.62) to 0.79 (0.64) over the southern (northern) slope.We note that across all parameterizations the correlation is higher over the southern slope than the northern slope, whereas the relative error shows the opposite pattern.This indicates that if one would further tune a 1 , it would be possible to further reduce the relative error metric, especially over the southern slope.However, here, for simplicity, we chose a 1 such that the diagnosed and parameterized diffusivities match in the central basin.We also note (not shown) that scaling by E eff improves the hybrid estimate (VI) such that its performance becomes similar to estimate (IX).
A summary of these findings is in order.The parameterized eddy velocity scale and eddy mixing length, whether based on L R or L T , do not reproduce their diagnosed counterparts (

̅̅̅̅̅̅̅̅̅ ̅ EKE
√ and L S ) over the continental slopes, as seen from Figure 2.But the failure of diagnosed eddy energy and L S in predicting a diffusivity over the slopes, particularly over the prograde slope in the north, also indicates that these eddy quantities do not give the full story.In particular, the eddy scale-which also appears to be related to EKE via the Eady length scale L E = EKE 1/ 2 σ 1 E -is not a good predictor of the effective mixing length.Accounting for the phase relationship between eddy velocity and buoyancy perturbations and, to a lesser degree, the eddy anisotropy brings the estimates much closer to the actual diffusivity.Our attempts at full parameterizations then clearly shows that L T is a much better choice than L R over both continental slopes.And the comparison between estimates VII and VIII even suggests that V par,T ∝ σ 2 E / β T is a slightly better predictor for the effective eddy mixing velocity over the slopes than the square root of eddy kinetic energy.Thus, a further consequence appears to be that the need for an explicit suppression factor (our E eff ) for the fully parameterized diffusivity becomes smaller.We leave further examination of this topic for later and instead carry on to see what effects the parameterized expression (IX) will have when used to actually operate in coarse-grained, non-eddying, simulations.

Performance in a Coarse-Resolution Channel Simulation
Before moving to realistic global simulations, we test the proposed parameterizations in a coarse-resolution channel setup.This setup has the same geometry as the high-resolution channel setup, but differs in resolution (from 2 to 32 km) and in that the GM-Redi parameterization scheme is activated.The model is forced and run similarly to the high-resolution setup, but some parameter settings, such as timestep and viscosity, are necessarily modified.Acknowledging that the channel setup is quite specific, and in anticipation of the global simulations to be studied next, we have kept tuning of the coarse-resolution channel setup to a minimum.
Figure 5 shows parameterized buoyancy diffusivities and the time-mean density field in three of the equilibrated simulations that had wide continental slopes but differing initial stratification.We chose to use the wide continental slope case because that is best resolved at 32 km resolution.We also show the corresponding diagnosed quantities from the corresponding high-resolution simulations for comparison (black lines).As in Figure 4, we show the three versions of the parameterized diffusivity (and corresponding density structure): one using the internal deformation scale (with a 1 = 8), one using the minimum between internal deformation scale and topographic Rhines scale (with a 1 = 8 and a T = 0.1) and, finally, one using the minimum and also applying the parameterized eddy efficiency factor (with a 1 = 32, a 2 = 1, and a T = 0.1).Choosing a T = 0.1 is consistent with the high-resolution model diagnostics as it provides a reasonable fit in terms of the order of magnitude between the full and parameterized versions of the topographic Rhines scale (see Figure S2 in Supporting Information S1).Together, constants a 1 and a T scale the diffusivity estimate by 0.08 (without E eff ) and 0.32 (with E eff ) over the sloping regions.These appear similar to the scaling coefficients used in Figure 4 (∼0.25,see figure caption) as well as previously reported values with mixing length theory: 0.33 by Wang and Stewart (2018) and 0.17 by Wei et al. (2022), both focusing on buoyancy diffusivity, as well as 0.2 by Wei and Wang (2021), although they focused on along isopycnal diffusivity.The tuning of a 1 was done to approximately match the mid-basin parameterized diffusivity with the corresponding diagnosed diffusivity in the high-resolution simulations.And yet, Figure 5 shows that the parameterized diffusivities have a clear north-south gradient in the magnitude whereas their diagnosed counterpart does not (panels a-c).The discrepancy is caused by a stronger difference in stratification between north and south at coarse resolution (see Figure S3 in Supporting Information S1) which directly impacts the internal deformation radius used by the parameterization over the central flat region.Despite this, we push forward and review the performance of the different parameterizations over the continental slopes.The deformation scale-based parameterization (Figures 5a-5c, orange line) clearly does worst, producing local diffusivity maxima over both slopes, as also seen in Figure 4 (blue line).Thus, there are only very weak lateral density fronts, or thermal wind shears, over the continental slopes.Essentially, the high parameterized diffusivities effectively wash out any density front there (see Figures 5d-5f for density structure and Figure S3 in Supporting Information S1 for the density gradients).This, it should be remembered, is exactly the effect one wishes to reduce with a slope-sensitive parameterization.
The run using a parameterization which selects the minimum of the two length scales does much better over both continental slopes where the topographic Rhines scale kicks in.With suppressed diffusivities, the density front which is set up by the topographic PV gradient is no longer washed out completely, especially in the north (see Figures 5g-5l for density structure and Figure S3 in Supporting Information S1 for the density gradients).The result is an enhanced thermal wind shear over the northern slope, albeit with a lower absolute strength than in the high-resolution simulation (about a factor two lower).In the south, where the vertical stratification is much weaker, the parameterization is not able to set up a strong thermal wind shear, although the location and strength of the surface density front has improved.Further scaling by the eddy efficiency E eff (Figure 5, green) enhances the diffusivity reduction slightly in the north, but not necessarily in the south.Therefore, the feedback to the resolved fields strengthens the baroclinic jet in the north further, but not in the south.
Clearly, there remains significant discrepancies between the low-resolution and high-resolution fields.Some may be due to higher levels of implicit (numerical) diffusion in the coarse-resolution simulations, but others likely reflect limitations in the parameterization.More extensive tuning, for example, based on matching the baroclinic transport over the slopes, would likely bring the coarse resolution simulations closer to the high-resolution fields.But the aim here has primarily been a qualitative examination, and even with a minimum of tuning the behavior over the slopes is robust, as the response to the different parameterizations is consistent across all cases.All in all, the above results are encouraging in that the Rhines-based parameterization is able to reduce wash-out of the density fronts over the slopes-particularly over the prograde slope in the north.However, although the channel setup is a reasonable test bed for development, it is extremely idealized and lacks multiple features from the real world (e.g., variable Coriolis parameter, uneven topography and complex atmospheric forcing).Therefore, we also test the slope-aware parameterization in the realistic global domain next.

Eddy Parameterization Adjustments
We carry out a control simulation and 5 different perturbation experiments.For simplicity, the focus will mostly be on a comparison between the control simulation and two of the perturbation experiments.All of these simulations operate with 2D diffusivities based on the depth-averaged Eady growth rate and a square length scale, as in Equation 14.The control run selects a length scale from the minimum of the internal deformation radius and the planetary Rhines scale (not the topographic Rhines scale).Then, in two distinct "topo" runs we (a) introduce the topographic Rhines scale as an additional length scale and (b) also turn on the eddy efficiency factor E eff .The OMIP "topo" runs then differ slightly from the coarse-resolution channel setup in the choice of constant scaling factors.The constant factor a 1 which scales the overall diffusivity magnitude is set to 3 and factor a 2 used in E eff is set to 1.In addition, we adjust L T by trying a T = 1 and a T = 0.5.We view these constants as tuning factors specific to one particular setup; for example, a T is impacted both by the quality of the parameterized Rhines scale (Figure S2 in Supporting Information S1) as well as the bottom topography data set.Here we calculated the topographic beta parameter based on the model's bathymetry.An improved approach might be to take a high-resolution bathymetry product and low-pass filter it up to the mesoscale-e.g. to the deformation radius or slightly above -before calculating the slopes.Further analysis is left for future studies, but in Figure S1 of Supporting Information S1 we show that the model bathymetry-based β T is a reasonable fit to high resolution bathymetry-based β T (using 15 s resolution bathymetry by Sandwell et al. (2022)) after filtering to the local deformation radius.Finally, in all runs the diffusivity magnitude is scaled down with a resolution function (Hallberg, 2013) when the deformation radius is resolved by the model grid.
For simplicity, the along-isopycnal (Redi) tracer diffusivity is set to be the same as the GM diffusivity.We have assessed the impact of this choice in a set of additional experiment included in Figure S4 of Supporting Information S1 and summarized the findings in Section 4.2.
To put the OMIP experiments in some context, it should be mentioned that the model settings for the control run are similar to the NorESM model version used in the latest Climate Model Intercomparison Project (CMIP6) except for some aspects of the GM diffusivity formulation.The CMIP6 version of the model included a mixing length formulation where the length scale was selected as the minimum of the internal deformation radius and the planetary Rhines scale-as in our control simulation.However, the local Eady growth rate was then evaluated at each model level, rendering a 3D profile for both eddy driven advection (GM) and for along isopycnal mixing (Redi).Finally, the scaling-based diffusivity was adjusted by a zonal velocity-dependent mean flow suppression following Ferrari and Nikurashin (2010), and as in the experiments here, a resolution function (Hallberg, 2013) was also used.
The lack of vertical structure of the 2D parameterization proposed here, turned out to be a clear deficiency in the global domain as our initial simulations showed an unrealistically strong sensitivity to bottom slopes in the low and mid-latitude deep ocean.For example, large reductions in the parameterized diffusivity across mid-ocean ridges were not seen in eddy-permitting studies that diagnosed eddy diffusivity in the global domain (e.g., Bachman et al., 2020).Therefore, to reduce the topographic impact on eddy fluxes in strongly stratified low and mid-latitude regions, we added an ad hoc "limiter" of topographic effects-based on the assumption that if the resolved flow does not feel the bottom then it is unlikely that mesoscale eddies would do so either.Specifically, the topographic Rhines scale is scaled by cos(α) 10 which rapidly increases the topographic Rhines scale when the angle α between the resolved depth-averaged flow and the bottom slope tangent vector deviate by more than ∼30°that is, when the resolved depth-averaged flow is not aligned with the bottom slope.

Model Response in the Global Domain
As expected, introducing the topographic Rhines scale leads to locally reduced diffusivities over sloping topography, as shown in Figure 6 (top row).The effect is enhanced at high latitudes with a ∼50% reduction over Arctic and Antarctic continental slopes.Bringing in the eddy efficiency E eff (see Figure S6 and Text S6 in Supporting Information S1 for an estimate of the E eff pattern) leads to additional and more severe diffusivity reduction globally, also away from topographic features (bottom row).This is also in agreement with recent studies (admittedly focusing on Redi mixing) that found that, at large scales, the scaling by mean-flow dependent suppression has the largest impact on diffusivity (Holmes et al., 2022;Stanley et al., 2020;W. Zhang & Wolfe, 2022).Note that in the tropics, the diffusivity is limited by the grid resolution function (Hallberg, 2013), that is, the diffusivity is reduced when the grid size is smaller than the local deformation radius.Therefore, the large relative reduction in tropical diffusivity is small in absolute terms and less important there as transport is dominated by the resolved flow.Finally, we note that a comparison between the top and the bottom rows in Figure 6 shows that in multiple continental slope regions, especially in the Arctic and around Antarctica, the eddy efficiency simply enhances the response seen with the topographic Rhine scale.Indeed, the diffusivity reduction due to introducing the topographic Rhines scale and due to eddy efficiency are close to linearly additive (not shown).
As the impact of eddy efficiency on diffusivity is more broad, its impact on flow speed, temperature, and salinity is also more widespread than the impact of the topographic Rhines scale alone.Table 3 collects bias reductions (relative to the control case) across 5 different experiments while Figures 7-9 show the spatial patterns for subsurface (100-200 m) current speed and temperature, as well as zonal-mean temperature and zonally-integrated overturning streamfunction anomalies for the two "topo" experiments that are in focus here (the overturning volume streamfunction is diagnosed by dividing the online-calculated overturning mass streamfunction with a constant reference density, ρ 0 = 1,000 kg m 3 ).We show results for the subsurface response since the surface response in these forced simulations is strongly forced by the non-responsive atmosphere.Both the topographic Rhines scale alone and its combination with eddy efficiency increase the mean kinetic energy of the resolved flow globally (at 100-200 m depth, by 2.7% and 10.5%, respectively).This increase is especially noticeable over sloping bathymetry where the two impacts contribute approximately equally to the overall increase (over slopes where β T > 5 ⋅ 10 10 m 1 s 1 kinetic energy at 100-200 m depth increases by 9.1% and 20.8%, respectively).The two modifications also warm the ocean below the global thermocline and cool the surface, reducing the overall temperature bias at depth.But they increase the temperature bias at the thermocline (Table 3; Figures 9a and 9c).
Figure 6.Anomalies from the control case in parameterized (depth-averaged) GM diffusivity due to implementation of (top row) the topographic Rhines scale and (bottom row) eddy efficiency in addition to the topographic Rhines scale.Red contours show the 1,000 m 2 s 1 isoline for diffusivity in the control case and light gray contours show areas in the tropics where the grid size is smaller than the internal deformation radius and therefore the resolution function (Hallberg, 2013) reducing the GM coefficient is in effect.Note.The observational data sets are the WOA 2018 climatologies for temperature (Locarnini et al., 2018) and salinity (Zweng et al., 2018).The experiment names correspond to parameterizations in Figure 4 as follows: "L T " and "0.5 ⋅ L T " correspond to (VIII), with a t = 1 and a t = 0.5, respectively; "L T and E eff " and "0.5 ⋅ L T and E eff " correspond to (IX), with a t = 1 and a t = 0.5, respectively.Experiment E eff does not have a counterpart in Figure 4, but uses parameterization like (IX), with the exception that L T is not considered as a length scale.
Overall, the mean overturning response in the "topo" runs is characterized by a positive (cyclonic) anomaly which implies that the Atlantic overturning cell and the Deacon cell in the Southern Ocean strengthen, whereas the Antarctic Bottom Water cell and the shallow surface overturning cells within the subtropical and subpolar gyres weaken.These changes generally reduce biases.The simulated strength of the Atlantic overturning at 26°N is 15.5 Sv in the control simulation, 17 Sv when topographic Rhines scale is considered, and 18 Sv with the addition of eddy efficiency, whereas the observational estimate from the RAPID array (∼26°N) is 17 ± 3.3 Sv (Frajka-Williams et al., 2019).The Antarctic bottom water cell at 32°S weakens from 26.0 Sv in the control simulation to 23.5 Sv with topographic Rhines scale and 20.3 Sv with addition of eddy efficiency, whereas inverse modeling suggest 20.9 ± 6.7 Sv (Lumpkin & Speer, 2007).The Deacon cell strengthens from 13.2 Sv in the control simulations to 15.2 Sv with the topographic Rhines scale and 18.4 Sv when eddy efficiency is considered, whereas previous modeling estimates (Döös et al., 2008) and observational estimates (Speer et al., 2000) suggest a strength of 20 and 20-25 Sv, respectively.
Some more specific impacts of the topographic Rhines scale and eddy efficiency are a poleward shift and strengthening of the boundary and slope currents, with E eff generally speeding up the boundary currents at locations where observations show the core of the currents (Figure 7, observed currents in black contours).Changes in the net volume transports in most key passages remain small (Table 4), but the results show a strengthening of the ACC (Drake Passage transport; reduced bias), a general enhancement of water exchange between the Arctic and mid-latitudes (opposing influence on the bias in different straits), and strengthening of the Gulf Stream (Florida-Bahamas strait transport, reduced bias).The spinup of the ACC is a direct consequence of reduced diffusivities, allowing for stronger thermal wind currents.In the northern North Atlantic, the current speed response is directly reflected in the temperature response as the Atlantic Water warms up along its path from the Nordic Seas to the Arctic (Figure 8, reduced bias).Despite the speed-up of the Gulf Stream off the North American coast, its observed turning around Grand Banks off Newfoundland is not reproduced.Due to this deficiency, the cold bias off Newfoundland strengthens (Figure 8).This cold bias is a long standing issue in coarse resolution ocean models (Tsujino et al., 2020) and reducing the diffusivity along the current path or along the shelf break clearly does not mitigate the bias.We speculate that, similar to the southern retrograde slope in the channel configuration and recent results on the Gulf Stream reported by Uchida et al. (2022), the eddy momentum flux convergence that is not included in the parameterization plays a crucial role in determining the current path.
The overall overturning response leads to increasing heat transport toward the northern hemisphere (Figure 10).The northern hemisphere subtropical peak in northward heat transport in the Atlantic basin (globally) is 0.83 PW (1.07 PW) in the control simulation, 0.91 PW (1.15 PW) when topographic Rhines scale is considered, and 1.00 PW (1.26 PW) with the addition of eddy efficiency, whereas Trenberth et al. (2019) estimate approximately 1.1 PW (1.6 PW).Breaking down the zonally-integrated impacts into resolved and parameterized eddy components illustrates how the reduced eddy mass transport across the ACC (Figures 10a and 10b) also leads to less southward heat transport (Figures 10c-10e) and therefore a cooling of the Southern Ocean surface, but also warming over the continental slopes (Figure 8).Both these effects reduce the bias in the model.Note that the heat transport response is dominated by eddy-driven advection with a smaller contribution due to the eddy diffusion (Figures 10d and  10e).In contrast to the Southern Ocean, in the northern mid-latitudes the overall northward mass and heat transport increase as the mean overturning spins up (Figures 10a and 10c; Figure 9 right panels) and the eddy contributions actually weaken (Figures 10b, 10d, and 10e).
We note again that in these OMIP experiments we have taken the eddy driven advection ("GM") and isopycnal mixing ("Redi") coefficients to be the same, making the origin of the response ambiguous.However, with a set of additional experiments (Figure S4 in Supporting Information S1) we show that the circulation response is mostly due to GM and, to a large extent, can be constructed as a sum of experiments where GM and Redi are changed separately (i.e., the response is linearly additive).The detailed temperature response, especially the thermocline bias, is sensitive to the treatment of the along-isopycnal mixing, but here the combined response is not linearly additive.The non-linear temperature response suggest that when developing and tuning a model system the changes to GM and Redi parameterizations should be made simultaneously.For temperature, black contours show the control case bias relative to the WOA observations in 0.25°C intervals (dashed for negative, solid for positive, the thick solid curve shows the zero contour).Therefore, whenever solid (dashed) contours surround blue (red) areas the bias to the observations is reduced.For the MOC, the contours show the control case MOC at 5 Sv intervals with the thick solid curve indicating the 0 Sv contour.Therefore solid (dashed) contours surrounding red (blue) indicates intensifying overturning.Similar to Figures 7 and 8, a student's t-test is applied here, but gray dots are now shown as all values are significant at 5% level.2020), for pure observational estimates see Koenig et al. (2014) and Donohue et al. (2016); and Florida-Bahamas Strait transport come from Larsen and Sanford (1985).

Discussion
Our study has focused on a relatively small range of parameterization choices, essentially (a) re-examining the topographic Rhines scale as a relevant mixing length and (b) checking the importance of an additional suppression factor which we have called the eddy efficiency E eff .The studies by Wang and Stewart (2020) and Wei et al. (2022) did a more comprehensive sweep over possible parameterization choices but did not analyze prograde and retrograde bottom slopes under one and the same framework, which has been the intention here.Also, to the best of our knowledge, the current OMIP simulations constitute the first assessment of the impacts of a topographically-aware GM parameterization in realistic global ocean models.As such, this work should be taken as a pragmatic investigation into what can be achieved with simple parameterization approaches applied to existing models that do not contain a prognostic eddy energy equation (which in itself requires parameterization choices).As with all parameterizations, the options examined here are far from perfect, and below we discuss some shortcomings and unresolved questions.

The Relevance of the Topographic Rhines Scale
Earlier idealized model studies have given conflicting evidence for the relevance of the topographic Rhines scale.Jansen et al. (2019) and Kong and Jansen (2021) reported that using a generalized Rhines scale which accounts for both planetary and topographic beta in their eddy parameterization of flows in an idealized ACC-like domain improved their model skill.More in line with our work here, the idealized channel studies of Wang and Stewart (2020) and Wei et al. (2022) found the topographic Rhines scale to be a useful choice over retrograde slopesbut not over prograde slopes.This conclusion was drawn, however, after an empirical slope-dependent prefactor was applied in the retrograde case but not in the prograde case.Both studies also constructed diffusivities from diagnosed depth-averaged EKE.In other words, they set the eddy velocity scale to be V =

√ and then defined
, that is, using the actual definition of the topographic Rhines scale.However, here we find that over both prograde and retrograde slopes, a full parameterization using Equation 14with L T = a T σ E /β T , better reproduces diagnosed diffusivities from our high-resolution simulations than partial a parameterization using Equation 5 with L T = ̅̅̅̅̅̅̅̅̅̅̅ V /β T √ (see Figure 4).Although not analyzed in detail, our hypothesis is that the full parameterization produces better results because it leads to a β 2 T dependence for the overall diffusivity, instead of the β 1/ 2 T dependence when using Equation 5 and the actual definition of L T .The different power dependence is important because β T varies by several orders of magnitude across the slopes whereas EKE and σ E vary by less than one order of magnitude (see e.g., Figure 2).

The Interpretation of E eff
The above results also suggest that some of the discrepancy between the pure scaling-based diffusivity and the diagnosed diffusivity is contained in the suppression factor which ours and other studies have pointed to.Although we don't try to identify the underlying dynamics here, the suppression is reflected in an imperfect phase relationship between eddy velocity and buoyancy perturbations, possibly causing larger suppression over prograde slopes.
It is worth nothing that our E eff may be related to the topographic Eady problem of Blumsack and Gierasch (1972).This connection becomes apparent if we evaluate the 2D version of Equation 21.We begin by setting U bc = U tw , where U tw is the top-to-bottom thermal wind shear (a 2D quantity).Then we first consider the slope region where the topographic Rhines scale will be the relevant length scale.So, here, V = σ 2 E/ β T , where σ E is now the depthaveraged (2D) Eady growth rate.Noting that in the Eady model, where both N 2 and ∂U g /∂z are constant, σ E = 0.3 ⋅ U tw /L R .This allows us to rewrite Equation 21as where a 3 is a modified tuning factor.Here δ = β T L 2 R / U tw is the slope parameter of Blumsack and Gierasch (1972) which measures the ratio between topographic and isopycnal slopes.Equation 23 is further supported by Figure S5 of Supporting Information S1 showing reasonable correspondence (up to a constant) between U 2 bc/ V 2 and δ 2 .This expression is interesting not only because it brings in the controlling parameter of the modified Eady problem but also for its similarity to the slope-dependent prefactor used by Wang and Stewart (2020) over retrograde slopes in the parameter regime where the bottom slope is not much larger than the isopycnal slope.Their prefactor F MLT (from their table 3) has the topographic delta parameter to the power of one in the denominator, in contrast to our squared power.But we suggest that the impact of sampling errors in the empirical fitting be studied in future studies before the correspondence is rejected.We also note that the similar studies of prograde fronts by Brink (2016) and Wei et al. (2022) found best fits using similar expressions but using topographic Burger number Bu in place of the delta parameter, where the two are related via Bu = (σ E/ f )δ.The latter study concluded that scalings using δ instead of Bu where not successful over prograde slopes.But, again, a comparison with our results is not straightforward since their diffusivities were constructed using diagnosed EKE while ours are fully parameterized.The relationship between δ-based and Bu-based formulations is also an obvious topic for future work.
Note, finally, that over the flat regions where the deformation radius will act as the relevant length scale, the 2D version of our efficiency factor becomes constant, in agreement with the behavior seen in Figure 3.In fact, the 2D version of E eff was able to qualitatively reproduce the observed eddy efficiency behavior in the idealized channel simulations, with some changes required for the tuning constants (not shown).We nonetheless chose to use the 3D version in the realistic OMIP simulations in anticipation of a more complex hydrography and flow field where the various assumptions of the Eady model can be expected to hold to an even lesser degree than in the channel model.Interior thickness PV gradients, for example, are expected to be small in systems that are only forced by Ekman pumping, as our channel model is (see e.g., Manucharyan & Stewart, 2022;Meneghello et al., 2021).In a real ocean, where for example, thermohaline forcing can produce interior PV gradients, the suppression of eddy efficiency will inevitably be governed by additional non-dimensional parameters beyond Blumsack and Gierasch (1972) δ (or, alternatively, the topographic Burger number).Such 3D effects, caused by thermohaline forcing in addition to wind stress, may also be the underlying reason for why E eff had a much bigger impact in the OMIP simulations than it did in the channel.

Summary and Conclusions
Efforts to include topographic effects into mesoscale eddy parameterizations are warranted, especially at high latitudes where observations show that hydrographic fronts are typically locked to topography.The very existence of such fronts along continental slopes and submarine ridges imply not merely topographic steering of large-scale currents but also suppression of lateral mixing across topography.Yet, despite all the observational evidence, as well as solid theoretical arguments for example, reduced growth rates and length scales of baroclinic instability over sloping topography, most eddy parameterizations still do not account for any bathymetric influence.
Here we have re-examined the relevance of the topographic Rhines scale in the mixing length approach to parameterizing the Gent-McWilliams buoyancy diffusivity which is used for eddy advection.Constructing diffusivities using the Eady growth rate and a parameterized version of the topographic Rhines scale reproduces an observed order-of-magnitude reduction in diffusivity over continental slopes in idealized channel simulations.
The simulations and analysis cover both prograde and retrograde continental slopes, representing mean flows in the same and opposite direction to topographic waves, respectively.Although differing in detail, both the and parameterized mixing suppression are of similar order of magnitude on both sides.The skill of the parameterization is enhanced further, at least over the prograde slope, when the diffusivity is multiplied by an eddy efficiency factor E eff that is sensitive to the strength of the mean-flow vertical shear relative to the parameterized eddy velocity scale.Finally, we find that selecting a smooth minimum of the topographic Rhines scale and the internal deformation radius for length scale gives good skill over the entire idealized channel domain.
The parameterization is then tested in a realistic global ocean simulation.Comparison with a simulation where topographic effects on diffusivities are not included suggests that the topography-aware parameterization enhances the sharpness of hydrographic fronts and, as such, strengthens the thermal wind shear in boundary currents.The improvement is particularly noticeable at high latitudes, but we also observe large impacts throughout the world ocean.The globally-averaged temperature and salinity bias reductions are in the range O(1%)-O(10%), with largest reductions seen in Southern Ocean temperatures and in Atlantic Water temperatures in the Arctic.However, existing low-latitude thermocline biases tend to increase.
The complex pattern of bias changes seen is not uncommon in a realistic global model, as bias reduction is very much a tuning exercise involving a range of free parameters associated with different parameterizations (e.g., eddy transport, vertical mixing and air-sea-ice fluxes).Our parameterization also has free parameters and, as is common, we found that the different model configurations, specifically different resolutions, might require different values for these.But we did not attempt a rigorous tuning, especially not for the dynamically complex OMIP simulations.Simply put, the focus at this stage has not been on a well-tuned realistic global simulation, but rather on illustrating possible impacts of a topography-aware eddy parameterization.
The suggested parameterization is clearly incomplete.The relatively large difference in importance of the efficiency factor E eff between the coarse-resolution channel simulations and the realistic OMIP simulations is one indication of this.A second one is the fact that we had to use an ad hoc limiter when applying this in the OMIP simulations.One key reason why a limiter had to be used is likely that we have been ignoring any vertical structure in eddy velocities and, ultimately, diffusivities.Fundamentally, the kinematic interaction with the bottom involves eddy bottom velocities, and a number of observations as well as theoretical arguments have indicated that these are often significantly smaller than surface or even depth-averaged eddy velocities (see e.g., de La Lama et al., 2016;Killworth, 1992;Lacasce, 2017;Wunsch, 1997).The topographic impact, under such considerations, would probably be smaller than if estimated with depth-averaged quantities.Future work clearly needs to be put on such vertical structure, for example, by taking an equivalent barotropic structure as a starting point (Killworth, 1992).We also observe that in our coarse-resolution channel simulations the flow remains too baroclinic, similar to the results presented by Kjellsson and Zanna (2017) and Yankovsky et al. (2022).Although addition of vertical structure to the buoyancy diffusivity might mitigate the issue, feeding the mean flow with vertically-distributed eddy energy (e.g., via a backscatter-type parameterizations) might be needed to resolve it (Yankovsky et al., 2022).
Another key topic which we have entirely neglected in this study is the impact of bottom roughness or corrugations on fluxes-and how such impact may be asymmetric with respect to the flow direction.As demonstrated by Wang and Stewart (2020), bottom roughness along a retrograde topographic slope can set up additional eddy buoyancy transport and, thus, form stresses due to arrested topographic waves.The dynamics governing such fluxes are likely distinct from those captured by our parameterizations here for smooth topography.The relevant eddy length scale, for example, is probably not the same as for transient eddies (Khani et al., 2019), and even coarse resolution models might be able to reproduce some of the largest standing meanders (Kong & Jansen, 2021).The application of standing Rossby wave theory (e.g., Abernathey & Cessi, 2014;A. L. Stewart et al., 2023) appears to give promising results on the planetary beta plane with a flat but rough bottom.A natural next step may therefore be to examine such ideas to the "topographic beta" problem, using for example, the idealized two-slope model used here.
Yet another issue ignored here is the role of lateral eddy momentum fluxes over continental slopes.As shown in Figure 1 and also highlighted in earlier studies (e.g., Manucharyan & Isachsen, 2019;Wang & Stewart, 2018), such fluxes bring wind momentum off the slopes to relatively flat regions where baroclinic instability kicks in to transfer the momentum to the ground below.The lateral momentum flux may be up-gradient in places and form eddy-driven jets, as seen offshore of the retrograde slope in our idealized simulations (Figure 1).As with eddy form stress, lateral momentum fluxes also appear to be impacted by corrugated bottoms, being associated with the formation of prograde jets near the bottom (Wang & Stewart, 2020).This last effect is again probably related to the formation of arrested topographic waves, as discussed by for example, Haidvogel and Brink (1986), as well as being linked to down-gradient PV diffusion in the finite-amplitude limit (Bretherton & Haidvogel, 1976;Vallis & Maltrud, 1993).
Finally, it's worth remembering that eddy transport, even of buoyancy, may be anisotropic.So what really needs to be parameterized is a diffusion tensor rather than a single scalar.Bachman et al. (2020) discussed such anisotropy of the Redi diffusion tensor and showed that at global scale the direction of the major axis of the tensor is well correlated with the mean flow direction and the minor axis is well correlated with the gradient of Ertel PV.In addition, Nummelin et al. (2021, Appendix A) suggested that the Ferrari and Nikurashin (2010) type of mean-flow suppression indeed suppresses the across-flow Redi mixing, but that the inverse of the same factor enhances mixing in the along-flow direction.It remains unclear whether our eddy efficiency factor-here primarily applied to buoyancy mixing-and the other empirical scaling factors (e.g., Wang & Stewart, 2020;Wei et al., 2022) act similarly (i.e., relate to tensor anisotropy) or if they indeed suppress the overall tensor magnitude.In other words, it remains a research question whether the mean flow and topography merely direct the eddy transport or if they impact the overall magnitude of the eddy transport.Nevertheless, if the tensor major axis is correlated with the mean flow (as suggested by Bachman et al. (2020)) -and if that mean flow transport dominates over eddy transport-then the focus on the minor axis is likely justified.
Even if important questions remain, and despite its many shortcomings, the relatively simple parameterization investigated here at least reduces an excessive washing out of hydrographic fronts over submarine ridges and continental slopes in ocean climate models-a known problem with eddy parameterizations that are insensitive of bathymetry.One of several important consequences of such adjustment is likely a more accurate representation of oceanic heat transport across Antarctic and Greenland continental slopes and onward to the great ice sheets whose melt rates depend intimately on such transport.On the shallow continental shelves, tides and other ageostrophic processes which we have neglected entirely here will also contribute.However, getting fluxes right across the strong fronts along the continental slopes is no less important.For this and other reasons, further scrutiny of all of the above unresolved issues related to mesoscale eddy transport and their impacts in both regional and global realistic simulations are much needed.

Figure 1 .
Figure 1.Cross section of zonally and temporally-averaged (a) mean zonal velocity (dashed black contours), mean density (dotted gray contours), E-P flux (gray arrows) and its horizontal component (meridional eddy flux of negative u-momentum; shading, note that the units have been scaled by 10 3 and here the bar denotes time-mean), (b) vertically-averaged zonal velocity and (c) vertically-averaged meridional buoyancy (temperature) diffusivity.In panels b and c the blue lines show results from the various experiments listed in Table 2.The black line is Experiment 3 and corresponds to the case shown in panel (a).Gray shading shows the location of the slope regions in the different simulations (where 300 m < H < 2,250 m).For some of the simulations the diffusivity lines are broken because of negative diffusivities that are not shown on the log scale.

Figure 2 .
Figure 2. Diagnosed length scales (panels on the left) and velocity scales (panels on the right) for all experiments.All measures are zonal and time averages that have been normalized.The top row (panels a and b) are normalized by the basinmean values (denoted by square brackets).Length scales in panels (c and e) are normalized by the deformation radius (L R ) and by the minimum of deformation radius and topographic Rhines scale (L min ), respectively.In panels d and f we normalize by the parameterized velocity scale, using length scales from c and e, respectively.Colors and line styles as in Figure 1.Gray shadings indicate the slope regions (similar to Figure 1) and vertical lines indicate the location of maxima in depth-averaged velocity in each experiment.

Figure 3 .
Figure 3. Measures of anisotropy and phase angle relationships: (a) eddy velocity anisotropy (A), (b) cosine of the phase angle between T′ and v′ and (c) the product of (a) and (b), as well as the parameterized eddy efficiency factors E eff (brown when using deformation radius, pink when using the topographic Rhines scale).In panel c we use a 2 = 10 for E eff .To match the midbasin values of E eff with A cos(θ), we scale E R eff with 0.35 and E T eff with 0.32.Colors and line styles as in Figure 1, and gray shadings indicate the slope regions.

Figure 4 .
Figure 4. (a-c) Partly-parameterized and (d-f) fully-parameterized across-slope buoyancy diffusivities.The left panels show across-basin profiles for Experiment 3 (Table2) whereas the right panels summarize the statistical comparison between the diagnosed and parameterized diffusivities across all experiments (b-c and e-f; relative error is defined as |(K par K diag )|/K par , where K par is one of the parameterization (I)-(IX); r 2 is based on linear regression using 200 points across all cases that is repeated 5,000 times).Boxes and whiskers come in pairs, with the one on the left (right) corresponding to the southern (northern) slope.Linear regressions are done over the slope regions only (gray shading; similar to Figure1).In panels a and d, the estimates I-IX are scaled by constant a 1 optimized to match the mid-basin diffusivity, and in addition we use a T = 0.5 (for VII-IX), and a 2 = 2 (for IX).The a 1 values are 0.1, 1.17, 0.09, 0.31, 3.77, 0.25, 0.66, 0.25, and 0.25 for estimates I-IX, respectively.In panels (b and e), the dashed gray lines correspond to 25%, 50%, 100% absolute relative error.

Figure 5 .
Figure 5. Buoyancy diffusivity (top row) and potential density anomaly (referenced to 0 dbar) in the upper 1,000 m of the coarse-resolution channel simulation compared to the high-resolution simulation (black line; coarse-grained to the coarse resolution grid).The different columns show experiments with different stratification such that the initial conditions are the same as for Exp 3, 6, and 9 in the left, middle and right columns, respectively (the initial stratification decreases to the right).For panels (dl), title indicates the parameterization.Bottom topography is indicated with a gray line.

Figure 7 .
Figure 7. Flow speed anomalies from the control case at 100-200 m depth due to implementation of: (top row) the topographic Rhines scale and (bottom row) eddy efficiency in addition to the topographic Rhines scale.Black contours show the 0.25 m s 1 isolines for observational estimate of the quasi-geostrophic current speed (Buongiorno Nardelli, 2020) in the same 100-200 m depth interval.Gray dots mark grid cells where 30 year mean is not significantly different from the control case at 5% significance level (student's t-test).

Figure 8 .
Figure 8. Temperature anomalies from the control case between 100 and 200 m depth due to implementation of: (top row) the topographic Rhines scale and (bottom row)eddy efficiency in addition to the topographic Rhines scale.Black contours show the ±1°C (solid/dashed) isoline for the control case bias relative to the WOA observations.Therefore, whenever solid (dashed) contours surround blue (red) areas the bias is reduced.Gray dots mark grid cells where 30 year mean is not significantly different from the control case at 5% significance level (student's t-test).

Figure 9 .
Figure9.Zonal-mean temperature anomalies (left panels) and global meridional overturning stream function (MOC) anomalies (right panels), relative to the control simulation, due to implementation of: (top row) the topographic Rhines scale and (bottom row) eddy efficiency in addition to the topographic Rhines scale.For temperature, black contours show the control case bias relative to the WOA observations in 0.25°C intervals (dashed for negative, solid for positive, the thick solid curve shows the zero contour).Therefore, whenever solid (dashed) contours surround blue (red) areas the bias to the observations is reduced.For the MOC, the contours show the control case MOC at 5 Sv intervals with the thick solid curve indicating the 0 Sv contour.Therefore solid (dashed) contours surrounding red (blue) indicates intensifying overturning.Similar to Figures7 and 8, a student's t-test is applied here, but gray dots are now shown as all values are significant at 5% level.

Figure 10 .
Figure 10.(a and b) Resolved and eddy contributions to the global meridional overturning circulation (MOC) and (c-e) to the global northward ocean heat transport (OHT).For the MOC we show the maximum (solid) and minimum (dashed) below 500 m to avoid the shallow surface overturning cells.For the OHT we show both advective and diffusive eddy contributions (panels d and e, respectively).

Table 1
Bergen Layered Ocean Model Model Constants for the Channel Simulations

Table 3
CORE-II Hydrography Bias (Root Mean Square Error) Reduction Compared to the Bias of the Control Case

Table 4
Observed and Simulated Current Transport in Selected Straits Note.The various perturbation experiments show percentage changes relative to the control case.The references for the observational values are as follows: Arctic Ocean gateway transports come from de Boer et al. (2018) with the original citations being Ingvaldsen et al. (2004) for Barents Sea Opening, Beszczynska-Möller et al. (2015) for Fram Strait, Curry et al. (2014) for Davis Strait (CAA), and Woodgate (2018); Woodgate et al. (2015) for Bering Strait; ACC transport come from Xu et al. (