## Abstract

The processes that determine the depth of the Southern Ocean thermocline are considered. In existing conceptual frameworks, the thermocline depth is determined by competition between the mean and eddy heat transport, with a contribution from the interaction with the stratification in the enclosed portion of the ocean. Using numerical simulations, this study examines the equilibration of an idealized circumpolar current with and without topography. The authors find that eddies are much more efficient when topography is present, leading to a shallower thermocline than in the flat case. A simple quasigeostrophic analytical model shows that the topographically induced standing wave increases the effective eddy diffusivity by increasing the local buoyancy gradients and lengthening the buoyancy contours across which the eddies transport heat. In addition to this local heat flux intensification, transient eddy heat fluxes are suppressed away from the topography, especially upstream, indicating that localized topography leads to local (absolute) baroclinic instability and its subsequent finite-amplitude equilibration, which extracts available potential energy very efficiently from the time-mean flow.

## 1. Introduction

Midlatitude gyre flows, confined within closed basins, produce a relatively shallow thermocline. In contrast, the Southern Ocean’s unique geometry permits the Antarctic Circumpolar Current (ACC) to circumnavigate the globe, accompanied by much deeper stratification. Many studies have shown that the stratification generated in the ACC pervades the global ocean below roughly 500-m depth (Toggweiler and Samuels 1995; Gnanadesikan 1999; Wolfe and Cessi 2010; Kamenkovich and Radko 2011; Nikurashin and Vallis 2012; Munday et al. 2013). Therefore, to understand the global deep stratification, it is necessary to understand the dynamics of the ACC and how its equilibrium responds to changes in forcing. There is ample evidence that the wind stress forcing in the Southern Ocean has increased over decadal time scales (Marshall 2003; Toggweiler 2009) and may have been drastically reduced during the Last Glacial Maximum (Toggweiler and Russell 2008), suggesting that the equilibration of the Southern Ocean is of relevance for a wide range of climate problems.

The dominant paradigm for understanding the ACC stratification involves a balance between buoyancy transport by wind-driven upwelling, which tends to steepen the isopycnal slope and deepen the stratification, and baroclinic eddy transport, which reduces the isopycnal slope (Karsten et al. 2002). Thus, the stratification of the ACC is part of the broad fundamental problem in geophysical fluid dynamics and planetary atmospheres, sometimes referred to as “baroclinic equilibration,” of determining the statistical properties of eddy heat transport as a function of large-scale parameters (Green 1970; Stone 1972; Held 1999; Schneider 2006; Jansen and Ferrari 2012). The classical approach of assuming adjustment to a marginally baroclinically unstable state (Stone 1972; Straub 1993) is inappropriate for the ACC, not only because unstable modes always exist in the continuously stratified case, but also because oceanic eddies are generally too small to bring the classical criticality parameter to unity (Jansen and Ferrari 2012): the zonally and time-averaged baroclinic velocity is of the order of 0.2 m s^{−1}, much larger than the phase speed of long baroclinic Rossby waves, which, with a deformation radius of about 15 km, is about 100 times smaller than the observed zonal-mean velocities.

Many eddy-permitting and eddy-resolving ACC modeling studies have examined the sensitivity of the isopycnal slope and associated thermal wind zonal transport on the wind stress forcing (Hallberg and Gnanadesikan 2001; Henning and Vallis 2005; Hallberg and Gnanadesikan 2006; Meredith and Hogg 2006; Hogg et al. 2008; Viebahn and Eden 2010; Farneti et al. 2010; Treguier et al. 2010; Abernathey et al. 2011; Meredith et al. 2012; Morrison and Hogg 2013; Munday et al. 2013). The sensitivity varies somewhat across models, and in all these studies, the eddy efficiency is a key parameter in determining the isopycnal slopes and resulting stratification depth.

Much of what is known about eddy equilibration comes from studies of baroclinic turbulence on a zonally symmetric beta plane (Held and Larichev 1996; Visbeck et al. 1997; Karsten et al. 2002; Thompson and Young 2006; Cessi 2008; Jansen and Ferrari 2012). In contrast, the ACC mesoscale turbulence is highly inhomogeneous in longitude, with pronounced “storm tracks” downstream of major topographic features such as the Kerguelen Plateau or Drake Passage. Most of the cross-frontal eddy exchange of mass and tracers occurs in such storm tracks (Treguier and McWilliams 1990; MacCready and Rhines 2001; Naveira-Garabato et al. 2011; Thompson and Sallée 2012). Accompanying these major topographic features are “stationary waves,” that is, meanders in the time-mean current with some characteristics of standing Rossby waves.^{1} To motivate our study, which is framed in terms of heat transport, in Fig. 1 we plot the vertically integrated divergence of the transient eddy heat flux in the ACC region, as calculated from the eddy-permitting Southern Ocean State Estimate (SOSE; Mazloff et al. 2010). We also plot the net mean (i.e., time averaged) and eddy (i.e., the departure from time average) heat transports across the ACC streamlines. (See section 2 for further details of the calculation.) The figure clearly illustrates that the eddy fluxes are an order of magnitude larger in the vicinity of the Kerguelen Plateau, Macquarie Ridge, East Pacific Rise, and Drake Passage, implying that these few regions make the dominant contribution to the net eddy heat transport across streamlines.

The main goal of this study is to better understand how such localized orography affects the efficiency of eddies in the baroclinic equilibration process. We approach this problem by comparing the equilibration of idealized channel models with a topographic ridge to those with a flat bottom. Although the flat-bottomed case is not a realistic ACC model in itself, the comparison highlights the important role of topography in modifying eddy efficiency. Our paper is organized as follows: In section 2, we introduce an idealized ACC-like problem and describe how the stratification depth depends on the winds and the eddy heat flux across circumpolar contours. In section 3, we present solutions to this problem obtained from a primitive equation, eddy-resolving numerical simulation over a wide range of wind stress magnitudes. To aid in the interpretation of the results, section 4 develops analytical solutions for a two-layer, quasigeostrophic model that demonstrates the mechanism by which standing waves enhance the efficiency of baroclinic equilibration. In section 5, we examine the detailed structure of the cross-stream heat flux in the simulation, revealing the intense localization downstream of topography and the suppression of mixing away from the storm track. In section 6, we discuss the nature of baroclinic instability and wave propagation when topography is present and suggest the importance of locally unstable modes. Conclusions are summarized in section 7.

## 2. Eddy heat transport and thermocline depth

### a. Description of an idealized problem

The goal of our study is to investigate how transient eddies determine the mean stratification in the presence of isolated topography in a simple context amenable to revealing the underlying physics. To this end, we study a highly idealized problem in which competition between wind and eddies completely determines the stratification. The key ingredients of this problem are as follows:

A zonally reentrant beta-plane channel domain, which permits a zonal current to develop. Using Cartesian coordinates (

*x*,*y*,*z*), the channel dimensions are (*L*_{x},*L*_{y},*H*).Westerly wind stress forcing, which drives an Eulerian-mean overturning. The form of the wind stress is

*τ*=*τ*_{0}sin(*πy*/*L*_{y}), vanishing at the northern and southern boundaries.Surface buoyancy restoring, which maintains a meridional sea surface temperature gradient. For simplicity, we make the gradient linear, with magnitude Δ

*θ*/*L*_{y}*.*A quasi-adiabatic interior, with negligibly weak diapycnal mixing.

For the case with topography, a large topographic obstruction in the abyss.

Given these conditions, the system will equilibrate with sloping isopycnals of the thermocline overlying a weakly stratified abyss. The slope of the isopycnals determines the depth of the thermocline at the northern boundary, which we call *h*. This situation is illustrated schematically in Fig. 2. The sloping isopycnals indicate the presence of available potential energy (APE), which can be transferred to eddy motion.

Our choice of confining the flow in a channel with solid walls, where no flux of heat and no normal flow and no-slip boundary conditions are applied, removes the possibility of having a residual overturning circulation, at least in the low diffusivity limit considered here. Thus, our model ACC approaches the limit of zero residual circulation described by Johnson and Bryden (1989) or Kuo et al. (2005), in which mean and eddy-induced advection cancel completely. The vanishing residual circulation limit is a useful idealization of the complete ACC, in which there is a nonzero residual flow, but where large cancellations between the mean and eddies nevertheless occur (Speer et al. 2000; Hallberg and Gnanadesikan 2006; Volkov et al. 2010). In this way, the stratification in the channel is simply determined by a balance between wind-driven advection of buoyancy by the Ekman circulation, which steepens isopycnals and creates APE, and eddy buoyancy advection, which removes APE (Karsten et al. 2002). More generally, the Southern Ocean stratification is determined by a three-term heat balance involving Ekman and eddy heat transport *and* the residual circulation, which in turn depends on remote processes outside the channel. Gnanadesikan (1999) put forth a model of the global ocean pycnocline with three interacting components: North Atlantic sinking, low-latitude diffusive upwelling, and a Southern Ocean component that involves both Ekman and eddy transport (see also Wolfe and Cessi 2010; Nikurashin and Vallis 2012; Shakespeare and Hogg 2012). Our simplified geometry avoids introducing the additional unknown residual overturning, allowing us to focus purely on the eddy behavior. Our study should be interpreted not as a complete model of the global deep stratification but as a refinement of the Southern Ocean component of Gnanadesikan (1999), which must be coupled with other components to understand the global problem.

### b. Scaling of the thermocline depth

We now discuss the relationship between eddy heat transport efficiency and thermocline depth in this idealized problem. For simplicity, consider first the flat-bottomed, zonally symmetric case. The mean meridional heat transport (MHT, ) across a latitude circle is given, in the Boussinesq approximation, by

where *ρ*_{0} and *c*_{p} are, respectively, the reference density and specific heat of seawater, *υ* is the meridional velocity, *θ* is the potential temperature, and the angle brackets 〈 〉 indicate a zonal and time average (Peixóto and Oort 1992). Since the total vertically integrated mass flux across a latitude circle (or any other circumpolar contour) must vanish, an arbitrary constant *θ*_{0} can be chosen without changing (de Szoeke and Levine 1981). By choosing this constant to be the temperature of the abyss (0°C), we can allow ourselves to ignore the contribution to the heat transport by the deep mean flows, such as the bottom Ekman flow, where *θ* ≃ *θ*_{0}.

In the adiabatic limit, and with no-flux conditions on the solid walls, must tend to zero with the diffusivity. We can use this constraint on the heat transport to derive a scaling for the thermocline depth *h*. We divide the heat transport into a part due to the zonal- and time-averaged meridional velocity _{mean} and a part due to the time-dependent motions _{eddy}. With our choice of *θ*_{0}, _{mean} is simply the heat transported by the upper Ekman layer:

that is, an equatorward heat transport determined solely by externally specified parameters. The eddy heat transport is

where we have denoted with *υ*_{g} ≡ *υ* − 〈*υ*〉 the velocity that departs from the zonal and time average. The subscript *g* indicates that this component of the flow is mostly due to geostrophic motions. In a flat-bottomed case, with no mean zonal pressure gradients or standing eddies, _{eddy} can be due only to transient eddy fluctuations. Furthermore, we expect 〈*υ*_{g}*θ*〉 to be almost zero below the thermocline, where there is little background temperature gradient for the eddies to stir. We therefore define a characteristic value of the eddy heat transport over the thermocline of depth *h* as

The quantity is a key parameter that measures the efficiency of the eddies at transporting heat. Using this definition, we have

The approximate vanishing of the total heat transport, that is, = _{mean} + _{eddy} ≈ 0, leads to a scaling of the thermocline depth near the northern boundary given by

The scaling for *h* is not completely satisfactory since it depends on the unknown eddy heat transport efficiency: *h* is determined by the competition between the restratifying effect of the eddies and the overturning due to the Ekman cell. Another way to interpret (6) is to consider 〈*υ*_{g}*θ*〉 to be specified by the closure theory of Gent and McWilliams (1990) as 〈*υ*_{g}*θ*〉 = −*K*〈*θ*〉_{y}, where the eddy flux is directed along isotherms and is proportional to the meridional buoyancy gradient, ensuring that the eddies act adiabatically to reduce APE. The term *K* is the Gent–McWilliams eddy transfer coefficient. With such a closure, the slope of the thermocline, *s* ≡ −〈*θ*〉_{y}/〈*θ*〉_{z}, is easily obtained (Karsten et al. 2002; Henning and Vallis 2005; Nikurashin and Vallis 2012) and is given by

From the slope *s*, the depth of each isotherm can be calculated given its surface value, and then *h* is defined as the depth of the lowest isotherm on the northern boundary of the channel (Marshall and Radko 2003). While the Gent–McWilliams closure is convenient for models, (6) is more general. Here the goal is to understand how topographically induced asymmetry affects the eddy transport efficiency .

When topography is present, the meridional heat flux across latitude circles (and the corresponding downward flux of momentum) is dominated by standing eddy components, rather than transient eddies (e.g., Treguier and McWilliams 1990; Wolff et al. 1991; Volkov et al. 2010), suggesting a qualitatively different balance. With the zonal and time average denoted by brackets, and the time average denoted by an overbar, the anomalies are defined as

such that *A* = 〈*A*〉(*y*, *z*) + *A*^{†}(*x*, *y*, *z*) + *A*′(*x*, *y*, *z*, *t*). The standing wave is associated with *A*^{†} and the transient component by *A*′. In a flat-bottom simulation with statistical symmetry in the *x* direction, *A*^{†} = 0.

With this decomposition, the MHT is given by

where _{mean} is defined in (2), and the standing and eddy components are given respectively by

The numerical results presented in the next section show that, as the winds are increased, _{SE} becomes increasingly dominant over _{TE}. Such behavior has been documented in highly realistic ocean climate simulations (Dufour et al. 2012; Zika et al. 2013), suggesting that standing eddies are central to the baroclinic equilibration of the Southern Ocean.

A complementary approach, as first demonstrated by de Szoeke and Levine (1981), is to average along a meridional coordinate that follows the time-mean meanders of the ACC, so as to remove the standing eddy contribution, leaving a two-term balance between ageostrophic (i.e., Ekman driven) and transient eddy heat fluxes (Marshall et al. 1993; Hallberg and Gnanadesikan 2001; MacCready and Rhines 2001; Viebahn and Eden 2012).

We follow de Szoeke and Levine (1981) and choose the depth-averaged temperature Θ as our streamwise coordinate:^{2}

The time-mean heat transport across any closed contour Θ_{0} of Θ(*x*, *y*) can be expressed as

where is the unit normal vector to Θ. The overbar indicates a time average. (In the second line, we have used the divergence theorem, yielding an expression that is much easier to evaluate in practice from a numerical model.) We can further decompose ^{Θ} into a time mean and transient eddy components by writing , where the primes indicate the deviations from the time mean, a standard Reynolds decomposition. Analogous to the zonally symmetric case, we respectively denote these two components for the mean (which is dominated by Ekman fluxes; de Szoeke and Levine 1981) and . The cross-stream heat fluxes from SOSE in Fig. 1 were calculated in this way (except that the barotropic transport streamfunction Ψ, rather than Θ, is the streamwise coordinate). Notice that the total MHT transport is a small residual left from a large cancellation between and .

In contrast to the zonal average, the streamwise average shows that the flat-bottomed and zonally asymmetric problems are both governed by the same fundamental balance of eddy and Ekman heat transport. The two perspectives are complementary, as demonstrated by Marshall et al. (1993), who show that the transient eddy flux across a streamline must approximately equal the *sum* of standing and transient eddy fluxes across an equivalent latitude circle. For heat fluxes, this means that . This equivalence is illustrated in Fig. 3, which sketches the mechanism for heat transport by standing eddies. Assuming that eddies flux heat down the mean temperature gradient (which is aligned with the mean streamlines), transient eddy heat flux convergence will occur in northward meanders, and divergence will occur in southward meanders. A water parcel moving along a mean streamline will therefore experience heating in the north and cooling in the south, resulting in a net southward heat transport by the mean flow. To explore this complex interaction between standing and transient eddies, we now turn to numerical simulations of the equilibration process with and without topography.

## 3. Numerical model results

### a. Model configuration

The goal of the model is to realize the system described above (and illustrated in Fig. 2) as simply as possible, while resolving the eddy heat flux. The code solves the hydrostatic Boussinesq equations in Cartesian coordinates on the *β* plane using the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al. 1997a,b). The model grid and numerical parameters are nearly identical to those described in Abernathey et al. (2011), to which the reader is referred for further details. The domain is a box *L*_{x} = 2000 km × *L*_{y} = 2000 km × *H* = 2985 m. The grid spacing is 5 km in the horizontal. There are 40 levels in the vertical, spaced 10 m apart at the surface and increasing to 200 m at depth. Linear bottom drag is applied in the bottom level of the model with a coefficient *r* = 1.1 × 10^{−3} m s^{−1}. With a deformation radius of approximately 15 km, this model adequately resolves the mesoscale dynamics.

The model’s potential temperature equation can be written as

Here *κ*_{h} is a spatially uniform horizontal diffusivity, and *κ*_{υ} is a vertical diffusivity. Advection is performed using a second-order moment scheme (Prather 1986). Explicit diffusivity (*κ*_{h} and *κ*_{υ}) is set to zero, and a detailed analysis has shown that the effective numerical diapycnal diffusivity in this model is weaker than 10^{−5} m^{2} s^{−1}, meaning the interior is almost adiabatic (Hill et al. 2012). However, the K-profile parameterization (KPP) scheme (Large et al. 1994) is employed to simulate the surface mixed layer, where *κ*_{υ} is greatly enhanced. The final term represents the surface forcing, active only in the top model level; *λ* is a temperature relaxation inverse time scale. The surface temperature *θ*_{s} is relaxed to a linear function of latitude of the form *θ** = Δ*θ*(*y*/*L*_{y}). The minimum temperature is 0°C, and we choose a maximum temperature Δ*θ* = 8°C. This leads to a maximum buoyancy contrast of Δ*b* = *gα*Δ*T* = 1.6 × 10^{−2} m s^{−2}. The relaxation time scale *λ*^{−1} is chosen to be 30 days (Haney 1971), which keeps the actual surface temperature very close to the prescribed profile. In practice, we do not achieve a truly vanishing MHT due to the diabatic effects in the surface layer. Our model contains a weak overturning cell in the top 200 m (the extent of the surface diabatic layer), similar to the one noted by Kuo et al. (2005). However, this is a small effect, and we remain very close to the limit of zero net MHT.

In simulations with topography, a Gaussian-shaped ridge is present in the middle of the domain. The motivation for this form was to capture the large meridional obstructions encountered by the ACC along its path, such as Kerguelen Plateau or the Scotia Arc. The depth in this case is given by

We selected *h*_{0} = 1000 m, about one-third of the total depth, and *σ* = 75 km, leading to a steep ridge. However, the topographic length scale *σ* is still large compared to the deformation radius of about 15 km in the middle of the domain. The high-resolution runs presented here are computationally demanding, limiting the choices of parameter exploration. We focus on the wind forcing, comparing the dependence of the stratification and zonal transport, with and without topography, rather than exploring different topographic geometries. Although the form of the topography determines the types of standing waves, broad large-scale topographies such as the one we have selected lead to the most energetic standing waves (Treguier and McWilliams 1990). An additional experiment in a domain twice as long (in the *x* direction) indicates that the domain size is sufficient to properly model the storm track.

The model equilibrates after about 100 yr of spinup. Snapshots of the temperature field from the equilibrated state are shown in Fig. 4. The time-mean isotherms are also superimposed. While both simulations contain mesoscale eddies, the figure illustrates how the flat-bottomed case is statistically symmetric in *x*, while the ridge case contains a standing wave in the time-mean temperature field.

### b. Meridional heat transport and thermocline depth

The heat balance is illustrated in Fig. 5 for both the flat-bottomed and ridge reference experiments (*τ*_{0} = 0.2 N m^{−2}). (Both can be considered streamwise-averaged views, since the Θ contours in the flat-bottomed case are zonally symmetric.) The upper panel demonstrates that _{mean} and _{eddy} are quite similar in both cases, with _{mean} remaining very close to the approximation defined in (2). The term _{eddy} compensates _{mean} almost completely, keeping the net very close to zero. The maximum value of _{mean} is 96 TW. Scaled up to the length of the full ACC, this corresponds to a transport of 1.1 PW, larger than SOSE by a factor of 2. The maximum value of is 82 TW, lower by about 15%.

The bottom panel shows the vertical structure of the eddy heat transport, along with the zonal or streamwise mean isotherms. In the ridge case, the thermocline is shallower, and the heat transport is confined to this shallower layer. As a practical matter, we define the thermocline depth *h* as the first moment of the 〈*θ*〉 profile via the expression

evaluated at the northern boundary, where the thermocline is deepest. This is a standard definition of the thermocline depth (Gnanadesikan 1999; Gnanadesikan et al. 2007; Munday et al. 2013). According to (16), the thermocline depth at the northern boundary is approximately 1200 m in the flat-bottom experiment and 1000 m in the ridge experiment. This means that , the eddy efficiency, is higher with topography present. Accordingly, the thermocline slope is shallower and there is much less APE with topography present: 1.2 PJ compared with 3.0 PJ in the flat-bottomed case. [APE was computed using the definition of Winters et al. (1995).]

The difference in efficiency increases as the wind stress increases. Figure 6 (top-left panel) shows the thermocline depth as a function of wind stress, for both the flat and ridge experiments, for the following values of *τ*_{0}: 0.0125, 0.025, 0.05, 0.1, 0.2, 0.4, and 0.8 N m^{−2}. This range constitutes six successive doublings of the wind stress. It is clear that the dependence of *h* on *τ*_{0} is significantly weaker in the ridge case. The difference is even more pronounced when comparing the APE (bottom-left panel); for the strongest winds, the APE is over 4 times greater without topography. This is because as the winds increase, the geostrophic flow, and the associated temperature transport , become more and more efficient at transporting heat poleward, leading to a weak dependence of *h* on *τ*_{0}. In all cases, increases more slowly than linearly with *τ*_{0}, leading to an increase, albeit weak, of *h* with *τ*_{0}. (We note that even with the weakest wind stress, *h* is much deeper than the mixed layer and is therefore not dependent on the KPP scheme.)

We also ran the ridge experiment with a domain of increased zonal extent (4000 km rather than 2000 km) for the reference value of *τ*_{0} = 0.2 N m^{−2}. One motivation for this experiment was to evaluate whether the relatively short domain was truncating the storm-track region and influencing the equilibration. The results are shown with the star symbol in Fig. 6; *h*, , and APE (an extensive quantity, therefore rescaled by a factor of 2) are nearly identical for the double-length run, confirming that our conclusions are not strongly dependent on the zonal extent of the domain. We shall consider this double-length simulation further in section 6, when we discuss the nature of baroclinic instability in the presence of topography.

It is informative to also consider the heat transport decomposition in terms of zonal averages, which distinguishes between the standing and transient eddy components, defined in (10). Figure 7 shows that, with topography, as the wind forcing increases, it is primarily _{SE} that compensates for the increasing value of _{mean}. Moving from the weakest to strongest winds, the ratio _{TE}/_{SE} decreases, with _{TE} going from completely dominant to completely negligible.

### c. Zonal volume transport

In addition to the stratification, *h* also contributes to the determination of the zonal transport *T*, defined as

where *u*_{b} is the bottom zonal velocity, and the thermal wind balance has been used to obtain the approximate expression in (17). Because 〈*θ*〉_{y} is fixed by the surface forcing, *T* ∝ *h*^{2}, as long as the bottom velocity is negligible. However, in calculations with a flat bottom, 〈*u*_{b}〉 can be very large. This is because the zonal momentum budget requires

where *r* is the coefficient of bottom drag, and *p*_{b} is the bottom pressure. If *h*_{bx} = 0, then the bottom drag is the only means to remove the zonal momentum imparted by the wind, and 〈*u*_{b}〉 is inversely proportional to *r* (and linearly proportional to *τ*) and thus large for weak friction. With topography, the wind stress is balanced by the orographic form stress (Munk and Palmén 1951; Olbers 1998), and 〈*u*_{b}〉 can be much smaller than the baroclinic component in the thermal wind balance. To illustrate this difference, Fig. 6 (bottom-right panel) shows the zonally averaged bottom zonal velocity averaged over the center half of the channel, as a function of *τ*_{0}, with and without the ridge: with the ridge, the bottom velocity is one order of magnitude smaller and largely independent of the bottom drag. The net result is that the zonally averaged zonal flow is much smaller when the ridge is present throughout the water column, and *T* is dominated by the second term on the rhs of (17). The total and thermal wind contributions are plotted in Fig. 8, showing this important difference.

## 4. Eddy heat flux enhancement by standing waves

To understand the enhancement of eddy heat transport in the presence of topography, the resulting decrease of thermocline depth *h*, and the reduction of bottom flow, it is useful to examine a quasigeostrophic (QG) two-layer model, forced by wind stress *τ* and dissipated by bottom drag. This model, in which eddy effects are parameterized by downgradient diffusion of potential vorticity, quantifies and explains how the standing wave contributes to enhancing the efficiency of the eddy equilibration.

The QG model is most easily analyzed by partitioning the flow into three components: a zonal and time average (denoted by angle brackets), a standing wave (denoted by a dagger), and a transient eddy component (denoted by a prime). The streamfunctions *ψ*_{n} and potential vorticities (PVs) *q*_{n} for layers *n* = 1, 2 are given by

The time average (denote by an overbar) of the primed quantities is also zero: . The relation between PV and the streamfunction is

The bottom topography *h*_{b} is the departure from the zonal average of the expression in (15), and . The layer depths are *H*_{1,2}, and *g*′ = Δ*b* is the maximum buoyancy contrast defined in section 3.

For the purpose of comparison with the primitive equations, it is useful to consider the QG dynamics of the waves and eddies to be embedded in the large-scale dynamics of the zonally and time-averaged flow, which we consider to obey the *planetary-scale* equations. The coupled system of the waves and eddies interacting with a planetary-scale flow allows one to consider the stratification parameter, here identified with *H*_{1} to be slowly varying in *y* and *determined as part of the solution*. We identify the *planetary-scale flow* with the time and zonally averaged flow, such that −〈*ψ*_{n}〉_{y} = *U*_{n}, with *U*_{n} slowly varying in *y*, so that 〈*q*_{1,2}〉_{y} = [*β* − *F*_{1,2}(*U*_{2,1} − *U*_{1,2})]. Formally, this requires a multiple-scale expansion that separates the synoptic-scale flow, here identified with and , from the planetary-scale flow 〈*ψ*〉. The details are given in Pedlosky (1987, their chapter 6.24), and for the two-layer model they lead to

In the classical oceanographic treatment of the planetary-scale component of the flow, the rectified contribution of transient eddies and standing waves is neglected. However, this is not appropriate for the periodic geometry considered here. Instead, the planetary-scale flow is determined by the zonal and time average of (23) and (24), which when integrated in *y* with suitable boundary conditions give the following constraints for *U*_{n}:

where *r* is the bottom drag coefficient. The planetary-scale momentum budget and the heat budgets are implicit in (25) and (26).

The standing wave is governed by

The system (25)–(28) can be simplified by adopting the framework considered by Hart (1979) in the barotropic context, where the topography, mean flow, and standing wave are slowly varying in *y*, so that the terms in (27)–(28) that are nonlinear in *ψ*^{†} can be neglected when considering the dynamics of the wave. As in Hart (1979), this approximation requires the amplitude of the wave to be smaller than the amplitude of the zonally averaged flow, so that , but . This ordering allows the explicit calculation of the standing wave and of the mean flow *U*_{1} and *U*_{2} [through (25) and (26)]. Consideration of the coupled QG/planetary-scale dynamics allows us to consider *H*_{1,2} to be also slowly varying in *y*.

The wave response can be easily obtained by parameterizing the time-mean transient eddy PV fluxes as downgradient diffusion of PV, that is,

The parameterization [(29)] ensures that the transient eddies act to damp the standing waves, as found in atmospheric models (Held and Ting 1990). It is also consistent with the diagnosed fluxes in the primitive equation results, although the diffusivity *K* is spatially variable in the simulations (see section 6).

### a. The standing wave response

With (29), and using the PV definition, (27) and (28) become, after neglecting the terms quadratic in and integrating in *x*,

It is clear that the standing wave is forced by the topography *h*_{b}, whose scale *σ* we assume to be larger than the deformation radius . We also assume that *Kσ*/*U*_{2} is small and that the wavelength is of the same order as *σ*. This scale separation permits us, to a first approximation, to neglect the relative vorticity of the wave compared to the vortex stretching terms multiplied by *F*_{1,2}. This limit is the opposite of the atmospheric case, where the deformation radius is often much *larger* than the topographic scales (Held et al. 1985). Although it is possible to solve (30) and (31) exactly, this approximation is more insightful and it gives, to the leading order,

that is, the response is equivalent barotropic. However, there is a correction to , which is 1/(*σ*^{2}*F*_{1}) smaller than the equivalent-barotropic component. The solvability condition is found at this higher order, or more simply by forming the equation for the barotropic mode, *H*_{1} [(30)] + *H*_{2} [(31)], for which relative vorticity cannot be neglected. This leads to

which is the equation for a driven, damped wave disturbance forced by the ridge. In anticipation that *U*_{2} will be much smaller than *U*_{1}, the wavelength of the stationary response is estimated to be . Using values appropriate for the primitive equation computations, the wavelength is of the order of 100 km, that is, of the same order as *σ*. Both these scales are larger than the baroclinic deformation radius, consistent with our approximations.

### b. The zonally averaged momentum and heat balances

To complete the solution and solve for , the values of *U*_{1} and *U*_{2} need to be supplied. These are determined by the zonally averaged PV relations of (25) and (26). The PV fluxes due to the standing wave are given by

Summing the depth-weighted equations *H*_{1} [(25)] + *H*_{2} [(26)], we obtain the total momentum budget

The momentum imparted by the wind is balanced by form stress and bottom drag, analogously to (18). This is one of the two relations used to determine *U*_{1} and *U*_{2}.

The second relation is obtained by considering the heat budget. The standing wave component, proportional to , vanishes to the leading order because of (32), but appears at the next order in (1/*σF*_{1})^{2}. It is possible to calculate , without evaluating the correction to explicitly, by multiplying (30) by and zonally averaging. Using the approximation [(32)], we find that the upper-layer PV flux^{3} is given by

This expression makes it clear that the heat flux carried by the standing wave is mediated by the heat flux by the transient eddies, here proportional to *K*. Using the parameterization [(29)], the combined heat transport by the standing wave and eddies is given by

This shows that the transient eddy heat transport is augmented by the topography and associated standing wave, just as we found in the numerical simulations. As shown in the appendix, this enhancement is best understood by examining the (parameterized) eddy heat transport across time-averaged temperature contours. The enhancement is due to two effects: increased local temperature gradients by the standing wave and increased arclength of the time-mean temperature contours.

### c. Solution of the two-layer model

Substituting (38) in (25), we find the zonally averaged and vertically integrated heat balance, given by

Thus, (39) and (36), together with (33) and the specification of the parameters *f*_{0}, *β*, *r*, *H*, and *τ*(*y*), determine the total flow in both layers.

In the coupled QG/planetary flow dynamics, the interface depth *H*_{1} can also be determined as part of the solution. Via the Margules relation, the term *f*_{0}(*U*_{1} − *U*_{2})/*g*′ corresponds to the *slope* of the interface between the two layers, which we associate with the slope of the lowest buoyancy surface bounding the thermocline in the primitive equation model. We thus define the planetary-scale slope of the isopycnal *s* as

Given this definition, we can rewrite (39) as

which illustrates that, assuming the transient eddy diffusivity (as measured by *K*) to be the same with and without the ridge, the slope of isopycnals is reduced by the presence of the standing wave. Relation (41) can be directly compared to (7) in the primitive equation model. Given the slope, the depth of the interface *H*_{1} is determined, up to a constant, by integrating *s* in *y*.

The slope [(41)] can be further simplified. Multiplying (33) by and zonally averaging, we obtain . We can then eliminate between (41) and (36) to find

The relations (33), (36), and (42) determine , *U*_{2}, and *s* or, equivalently, , *U*_{1}, *U*_{2}, and *H*_{1}. The value of *H*_{1} at *y* = 0 needs to be specified, along with the external forcing and geometrical parameters.

The solutions of the coupled system (33), (36), and (42) are easily obtained numerically. We illustrate the solutions by showing the maximum of *U*_{2} in the domain as a function of the wind stress amplitude as a solid line on the bottom-right panel of Fig. 6. The lines in Fig. 9 show the solutions for *H*_{1}(*y* = *L*) with (dashed) and without (solid) the ridge as a function of *τ*_{0}, assuming *τ* = *τ*_{0} sin(*πy*/*L*_{y}) and *H*_{1}(*y* = 0) = 0, a choice well beyond the range of validity of the QG approximation. In this calculation, we assume the form *K* = *K*_{ref}(*τ*_{0}/*τ*_{ref})^{a}, where *α*, *K*_{ref}, and *τ*_{ref} are constants. The values of the constants are chosen to best fit the primitive equation results of *h* as a function of *τ*_{0} for the flat-bottom case, with *α* ≃ . The QG model cannot predict this relationship, but it is necessary to account for the strong dependence of the overall eddy mixing rate on the wind power input (Cessi 2008; Abernathey et al. 2011). For the same values of *K*, the model predicts a substantially *reduced* value of *H*_{1} and a weaker dependence of *H*_{1} on *τ*_{0} when topography is present. Both results are in qualitative agreement with the simulations, showing how the combination of transient and standing eddies leads to a more efficient overall heat transport. Although the relation between *H*_{1} and *τ*_{0} does not agree *quantitatively* with the eddy-resolving primitive equation computations, the agreement is encouraging given the simplifying assumptions of the QG model (two-layers, constant eddy diffusivity, no mixed layer, etc.). To further compare the QG and primitive equation solution, we contour the total time-averaged interface height in Fig. 10, for the reference case *τ* = 0.2 sin(*πy*/*L*_{y}) (in N m^{−2}). The general pattern is captured by the QG model, but the response downstream is damped compared to the PE model (cf. with the gray contours of Fig. 11).

The QG model also correctly predicts the small amplitude of the bottom flow, so that the bottom drag term gives a small contribution to the momentum budget relative to the form drag and its dependence on the wind (solid line in Fig. 6, lower-right panel). For small bottom drag, *U*_{2} is insensitive to the value of *r*, and thus the slope varies almost linearly with *r*. The strong dependence of some aspects of the solution on *r* is a problematic aspect of the two-layer model.^{4}

Notice that assuming the same value of *K* with or without topography leads to values of *h* that are comparatively *smaller* than those obtained in the eddy-resolving computations in the ridge case. The next section shows that the isolated ridge, while enhancing the eddy diffusivity near the ridge, *suppresses it* in the rest of domain, leading to an average diffusivity that is smaller than in the flat case. This additional response complicates the equilibration process.

## 5. Local cross-stream heat flux

In the flat-bottom case, the eddy statistics are homogeneous in the *x* direction, and the cross-stream heat fluxes (Ekman and transient eddy) are distributed evenly along the front. Prior studies of idealized circumpolar currents with topographic ridges (MacCready and Rhines 2001; Hallberg and Gnanadesikan 2001; Thompson 2010) have shown that eddy thickness fluxes (related to the eddy buoyancy flux) are concentrated in regions near and downstream of topographic ridges. A similar conclusion was reached by Thompson and Sallée (2012) in an analysis of altimetric data; they found that Lagrangian trajectories cross the ACC fronts preferentially in a few locations downstream of major topographic features such as the Drake Passage or Kerguelen Plateau. We find the same result in our simplified eddy-resolving computations: eddy heat fluxes are concentrated near and *downstream* of the ridge, unlike the QG prediction where the eddy fluxes are enhanced directly *over* the ridge, in direct proportion to the local gradient.

To explore the distribution of the eddy heat flux along the front, and its localization downstream of the ridge, we analyze the vertically integrated temperature flux

in the (*x*, *y*) plane. The two separate fields **F**_{mean} and **F**_{eddy} correspond to the components due to the time mean and the time-fluctuating flow, respectively. Dotted with and integrated along the Θ contours, these components correspond to and , the two components of the cross-stream heat transport identified in section 2. The raw flux **F** is generally dominated by the rotational component of the eddy flux, which does not contribute to the integrated cross-stream flux, but can obscure the physics of cross-stream transport (Marshall and Shutts 1981; Illari and Marshall 1983; Jayne and Marotzke 2002; Bishop et al. 2013). This rotational part can be eliminated by solving the elliptic Poisson problem ∇^{2}*ϕ*(*x*, *y*) = **∇** · **F** subject to the boundary condition *ϕ*_{y} = 0 at the northern and southern walls. As shown by Fox-Kemper et al. (2003), this decomposition is not unique; the choice of boundary condition serves to highlight the cross-stream portion of the flux. The divergent component of **F** is then given by **F**^{div}(*x*, *y*) = **∇***ϕ*. We solve for *ϕ* numerically using an algebraic multigrid solver (this python-based solver is freely available online at https://code.google.com/p/pyamg), separating the mean and eddy parts of *ϕ* and **F**^{div}.

We plot in Fig. 11 (left panels) as arrows in the (*x*, *y*) plane, for both flat and ridge experiments. The projection of the flux normal to Θ contours, , is in color. The divergent flux is largely perpendicular to the Θ contours (in gray) and is nearly entirely downgradient (Marshall and Shutts 1981). In comparison with the flat case, the cross-stream flux in the ridge case is much more localized, occurring mostly in the vicinity of the strong meander downstream of the ridge. In fact, close inspection of the arrows in Fig. 11 reveals that is mostly a *zonal* flux across the Θ contours running north–south. These zonal fluxes go in both directions out and away from the trough of the standing wave, consistent with the cartoon in Fig. 3. The term is equal and opposite to the eddy flux, and it is not shown. We note that locally there are both Ekman and geostrophic contributions to , but the geostrophic component vanishes when integrated along each contour (de Szoeke and Levine 1981).

From the divergent eddy heat flux, it is possible to construct a local cross-stream eddy diffusivity. We define this diffusivity as

where *H* is the full depth. This quantity measures the local efficiency of eddies at transporting heat across the Θ contours. The term is plotted in Fig. 11 (right panels). For the flat-bottom experiment, is zonally uniform, peaking around 4000 m^{2} s^{−1} in the northern part of the domain. For the ridge, the region of highest diffusivity is downstream of the ridge, particularly in the trough of the standing meander, where diffusivities exceed 5000 m^{2} s^{−1}, while it is largely suppressed elsewhere. The net result is that the diffusivity is on average *smaller* in the ridge experiment; nevertheless, the efficiency of the eddy transfer, as measured by is clearly larger. This is because the region of strong coincides with a region of large |**∇**Θ|, giving a high correlation and a more effective local cross-stream flux. It is now clear why the QG model, which assumes that the diffusivity away from the ridge is the same with and without topography, leads to a prediction for *h* that is too small: the suppression of eddy diffusivity away from the ridge is not taken into account.

### Dependence of localization on winds

To illustrate the localization of the eddy fluxes as a function of the wind, we show in Fig. 12 (right panel) the cross-stream divergent flux, , as a function of *S* and *τ*_{0}, where *S* is the arclength along a Θ contour. We select the Θ = 0.96° contour, which is near where the maximum flux occurs, and normalize by the maximum value along the contour. In this way, the position of the maximum flux relative to the ridge peak is apparent, rather than its amplitude, which obviously increases as the wind stress is increased (cf. Fig. 5). We also include the double-length experiment with *τ*_{0} = 0.2 N m^{−2}. The fact that the maximum value of *S* increases with *τ*_{0} reflects the previously mentioned increase in the total arclength of Θ contours with the standing wave amplitude.

The maximum cross-stream heat flux is always downstream of the maximum of |**∇**Θ|, whose value along the same contour (also normalized by its maximum along the contour) is shown on the left panel of Fig. 12. It is clear that as the wind and amplitude of the stationary wave increase, the cross-stream heat flux becomes more localized downstream of the ridge, especially in the double-length computation. Moreover, the distributed portion of the cross-stream heat flux away from the ridge, characteristic of the flat-bottom case (cf. Fig. 11), which persists in the low wind regime, essentially vanishes for *τ*_{0} ≥ 0.05, leaving only the localized signal downstream of the ridge (and of |**∇**Θ|).

## 6. Local versus global eddy growth

The localization of the eddy fluxes just downstream of the ridge should be contrasted with their homogeneous nature in the flat-bottom case (cf. Fig. 11). This qualitative difference is associated with fundamentally different propagating properties of the eddies. We illustrate the eddy propagation in Fig. 13, a Hovmoeller diagram that shows the surface temperature anomalies (relative to the time mean) at *y* = 1000 km as a function of *x* and *t*: with the flat bottom (top panel), eddies propagate at a speed that is intermediate between the surface zonally averaged velocity (dashed black line) and the vertically and zonally averaged velocity (black solid line); with the ridge (lower panel), the transient eddies are almost stationary, moving eastward at a speed no larger than the vertically and zonally averaged velocity (black solid line) and much smaller than the surface velocity (dashed black line).

This difference in propagation indicates that in the flat case the eddies are generated and maintained through convective (global) instability, while with the ridge eddies are generated and maintained through absolute (local) instability (Merkine 1977; Pierrehumbert 1984). With topography, the generation of a standing meander locally increases the horizontal buoyancy gradients in the vicinity of the ridge (cf. the gray isotherms in Fig. 11), providing a local source of available potential energy that eddies can release. Furthermore, topography with a zonal slope reduces the stabilizing effect of *β* by orienting the wavenumber of the most unstable mode in a direction with a meridional, as well as a zonal, component (Chen and Kamenkovich 2013). The net result is a localized instability with larger growth rates than in the flat case and suppression of eddy growth away from the ridge.

An additional requirement for local growth is that the mean flow is slow enough to keep the eddies in place as they grow. The ridge reduces the speed of the mean zonal flow by reducing the bottom component relative to the flat case (cf. the bottom-right panel of Fig. 6). The dependence of local instability on the translational properties of the mean flow, as well as on its baroclinicity, is in sharp contrast with the global instability, where the growth rate depends on the shear, the stratification, and *β*, but not on the depth-averaged flow.

In summary, with the ridge, mean buoyancy gradients are locally enhanced and eddies propagate more slowly than in the flat case. Because there is a weak eastward propagation, the maximum in eddy activity is found downstream of the ridge. The difference in absolute versus convective instability is especially apparent in the initial transient development (Fig. 14, upper panel): in the flat case (upper-left panel of Fig. 14), there is a slow development of classical baroclinic instability beginning with “elevator” modes with purely meridional motion and eventually the development of a slow secondary instability in the orthogonal direction (Berloff et al. 2009) leading to finite-amplitude eddies. With the ridge (upper-right panel of Fig. 14), a stationary wave is immediately formed, with eddies quickly reaching finite amplitude in the lee of the ridge.

In the steady, equilibrated state, the source for eddy kinetic energy is given by the conversion term (e.g., Cessi et al. 2006). From the bottom panel of Fig. 14, we can see that this term is highly localized downstream of the ridge, indicating that the local eddy growth seen in the initial transient development remains a feature of the equilibrated state. It is interesting that the maximum is not collocated with the maximum cross-stream heat transport (Fig. 11), which is farther downstream; this reinforces the notion that downstream propagation is responsible for the downstream eddy heat flux maximum.

As the eddies equilibrate in the ridge case, there is an additional positive feedback that enhances their local growth; the poleward heat transport by eddies restratifies the interior, reducing the vertical extent of the zonal shear and thus the baroclinic component of the mean eastward zonal flow, proportional to *h*. This process further slows down the mean flow, reducing the eastward propagation of eddies and allowing continued extraction of mean flow energy into the eddy component.

## 7. Discussion and conclusions

We have explored the equilibration of an idealized baroclinic current with and without a topographic ridge, with the goal of understanding how zonal asymmetry affects the baroclinic equilibration process. In our simplified experiments, in which the interior of the ocean is quasi adiabatic, the thermocline depth is determined by the competition between poleward cross-frontal heat transport by the geostrophic eddies and the equatorward heat transport by the Ekman circulation. We find that, with localized topography, the eddy field accomplishes the same cross-frontal heat transport as in the flat case, but over a shallower layer and a narrower horizontal region. In this sense the geostrophic turbulence is “more efficient” with a ridge. A simple two-layer QG model partially explains the mechanism for this enhancement: the presence of a standing wave leads to 1) a stronger frontal temperature gradient and 2) increased arclength of time-mean temperature contours. Both of these factors allow the same amount of heat to be transported in a shallower layer. By solving for the standing wave amplitude, the QG model makes a quantitative prediction for the enhancement factor based solely on the external parameters.

However, the picture is complicated by differences in transient eddy behavior. Overall, the transient eddy diffusivity is weaker in the presence of the ridge; the cross-frontal eddy flux is concentrated in a narrow storm track and suppressed elsewhere. Also, the localization itself increases as a function of the winds. An explanation for these differences is in the baroclinic instability mechanism generating transient eddies: global eddy growth and equilibration versus a local growth. The geostrophic turbulence of the zonally symmetric, flat-bottomed channel can be viewed as a finite-amplitude equilibration of the classic global (or convective) baroclinic instability problems posed by Charney (1947) or Phillips (1951). Instead, the ridge experiments illustrate the nonlinear equilibration of local (or absolute) instability discussed by Pierrehumbert (1984), where eddy growth is suppressed away from the localized region of enhanced baroclinicity. In the local instability problem, the growth rate of eddies depends not only on the local baroclinic shear, but also on the vertically averaged zonal flow; a fast, vertically averaged mean flow sweeps the disturbances away from the region of baroclinicity before they can extract energy.

These considerations suggest why eddy fluxes in the storm-track region are exceedingly difficult to parameterize using existing frameworks (Hallberg and Gnanadesikan 2001, 2006). The emergence of eddies from local (or absolute) rather than global (convective) instability indicates that any parameterization of the eddy heat transport would have to take into account not just the local baroclinicity (i.e., shear) and stratification, but also the vertically averaged mean velocity, since this is an important parameter for absolute growth (Pierrehumbert 1984). Furthermore, the suppression of divergent eddy heat fluxes away from the ridge, together with the fact that the maximum eddy flux and mean gradient are not collocated, indicates that eddy generation and dissipation are nonlocal in space, as has been noted in western boundary current regions (Wilson and Williams 2004; Grooms et al. 2013). The inability of existing parameterizations to account for local instability and nonlocal eddy life cycles constitutes the main obstacle toward a more complete theory of baroclinic equilibration in the presence of large topography and the more general problem of inhomogenous geostrophic turbulence.

## Acknowledgments

The manuscript was greatly improved thanks to input by Andrew Thompson, Andrew Hogg, and an anonymous reviewer. Support by the Office of Science (BER), U.S. Department of Energy, Grant DE-SC0005100 is gratefully acknowledged.

### APPENDIX

#### An Alternative Derivation of (38)

The expression (38) can also be obtained by considering the heat transport across a vertically integrated time-mean temperature contour . In these coordinates is entirely due to the heat flux by transient eddies, that is,

where is the unit vector normal to the contour. Thus, we have

In the case at hand, we have

Using (32) we find

which shows that the eddy heat flux is enhanced near the ridge due to the locally increased gradients. To recover the expression [(38)] which refers to the zonal average, we need to integrate along contours of constant . This requires calculating the elemental arclength, along a contour, which satisfies so that

Thus, the arclength along the time-averaged temperature is lengthened by the presence of the standing wave. Finally, denoting with *C* the contour at constant temperature, we have that the “streamwise-averaged” heat flux due to the transient eddies is

which is the same expression as (38).

## REFERENCES

*Ocean Circulation: Mechanisms and Impacts, Geophys. Monogr.,*Vol. 173, Amer. Geophys. Union, 19–32.

*Tellus,*

**51A,**59–70, doi:.

*Geophys. Res. Lett.,*

**33,**L16608, doi:.

*J. Geophys. Res.,*

**116,**C09019, doi:.

*Geophysical Fluid Dynamics*. 2nd ed. Springer-Verlag, 710 pp.

*Physics of Climate.*American Institute of Physics, 520 pp.

## Footnotes

This article is included in the Ocean Turbulence Special Collection.

^{1}

Here we use the terms standing eddies and stationary waves interchangeably.

^{2}

Other choices of streamwise coordinate are possible, such as time-mean contours of the Montgomery potential in each isopycnal layer (Marshall et al. 1993; Hallberg and Gnanadesikan 2001; MacCready and Rhines 2001) or the barotropic streamfunction (Viebahn and Eden 2012).

^{3}

The vertically integrated heat flux is given by , and it is proportional to the PV flux.

^{4}

As in previous studies (Charney and DeVore 1979; Hart 1979; Charney and Straus 1980), the solutions of (33), (36), and (42) exhibit multiple equilibria for a large range of parameters. While multiple regimes are interesting in their own right, they were not found in our eddy-resolving computations. In the QG model, multiple equilibria can be eliminated by increasing the bottom drag parameter, and this is what we do.