## Abstract

Submesoscale processes have been extensively studied in observations and simulations of fronts. Recent idealized simulations show that submesoscale instabilities also occur in baroclinic mesoscale cyclones and anticyclones. The instabilities in the anticyclone grow faster and at coarser grid resolution than in the cyclone. The instabilities lead to larger restratification in the anticyclone than in the cyclone. The instabilities also lead to changes in the mean azimuthal jet around the anticyclone from 2-km resolution, but a similar effect only occurs in the cyclone at 0.25-km resolution. A numerical passive tracer experiment shows that submesoscale instabilities lead to deeper subduction in the interior of anticyclonic than cyclonic eddies because of outcropping isopycnals extending deeper into the thermocline in anticyclones. An energetic analysis suggests that both vertical shear production and vertical buoyancy fluxes are important in anticyclones but primarily vertical buoyancy fluxes occur in cyclones at these resolutions. The energy sources and sinks vary azimuthally around the eddies caused by the asymmetric effects of the Ekman buoyancy flux. Glider transects of a mesoscale anticyclone in the Tasman Sea show that water with low stratification and high oxygen concentrations is found in an anticyclone, in a manner that may be consistent with the model predictions for submesoscale subduction in mesoscale eddies.

## 1. Introduction

Submesoscale processes at fronts in the ocean mixed layer have been the subject of intense study in recent years (e.g., Boccaletti et al. 2007; Callies et al. 2015; Capet et al. 2008a; D’Asaro et al. 2011; Fox-Kemper et al. 2008; Hamlington et al. 2014; Haney et al. 2015; Taylor and Ferrari 2009; Thomas 2005; Thomas et al. 2008, 2013). These studies have identified a range of processes that may be active at mixed layer fronts including mixed layer baroclinic instability, symmetric instability, lateral shear instability, and frontogenesis.

Submesoscale processes at fronts have been shown to have a number of consequences. Mixed layer baroclinic instability leads to restratification of fronts (Boccaletti et al. 2007; Fox-Kemper et al. 2008; Nurser and Zhang 2000). Symmetric instability leads to the viscous dissipation of larger-scale geostrophic flows by inducing a downscale pathway for the kinetic energy of the geostrophic flow (Taylor and Ferrari 2009, 2010; Thomas and Taylor 2010). Beyond the canonical baroclinic and symmetric modes (Stone 1966), mixed modes of instability with characteristics between baroclinic and symmetric instability are also likely to be of importance (Stamper and Taylor 2017). Frontogenesis triggered by a larger-scale strain or by the effect of viscosity acting on the geostrophic flow leads to a sharpening of fronts (Hoskins and Bretherton 1972; McWilliams et al. 2009). Lateral shear instability can lead to the dissipation of the kinetic energy in submesoscale filaments (Gula et al. 2014). The Ekman transport at a front can also affect frontal properties by changing the mixed layer stratification (Thomas 2005) or driving frontogenesis (Thomas and Lee 2005).

The fluxes of tracers such as nutrients between the mixed layer and thermocline at fronts are also thought to be affected by submesoscale processes (Lévy et al. 2012). Smith et al. (2016) show that in a large-eddy simulation that permits both submesoscale instabilities and Langmuir turbulence, vertical fluxes of passive tracers may be inhibited by submesoscale restratification. On the other hand, Smith et al. (2016) find that regions of negative potential vorticity (defined below) in the simulations are areas of countergradient vertical diffusion of passive tracers and suggest that symmetric instability may thus act to counter turbulent mixing.

While research on submesoscale processes has concentrated on fronts, ocean eddies can also have strong lateral buoyancy gradients in the surface mixed layer. In addition, the vertical displacement of isopycnals in eddies means that they have an isopycnal pathway from the surface to the ocean interior. Previous research has focused on processes that can lead to variations in the isopycnal displacements associated with eddies at the eddy scale. Dewar and Flierl (1987), for example, show that the variation in wind stress across an eddy leads to vertical Ekman velocities that attenuate the eddy. McGillicuddy et al. (1998) show that the isopycnal displacements during the formation of cyclonic eddies lead to nutrient transport into the surface layer. Chelton et al. (2004) find that the acceleration of the surface wind stress across an eddy with a warm sea surface temperature anomaly leads to a curl in the wind stress that can generate Ekman pumping inside the eddy. A simulation of a mesoscale cyclone subject to strong surface cooling by Legg et al. (1998) shows that convective, symmetric, and baroclinic instabilities redistribute momentum and buoyancy before ultimately leading to the breakup of the cyclone.

The present work forms part of the Ocean Surface Mixing, Ocean Submesoscale Interaction Study (OSMOSIS). The OSMOSIS study region is away from major frontal systems where much work on submesoscale processes has occurred (e.g., D’Asaro et al. 2011; Thomas et al. 2013) and so aims to understand whether ideas developed at larger fronts are more generally applicable. Observations from the OSMOSIS project show that such open-ocean regions are rich in submesoscale activity, particularly in wintertime (Buckingham et al. 2016; Damerell et al. 2016; Thompson et al. 2016).

An idealized model of a wind-forced open-ocean domain is employed in Brannigan et al. (2015) to examine the seasonal cycle of submesoscale flows. This model is designed as an analog of the OSMOSIS observation site in the North Atlantic where mean flows are weak and the kinetic energy budget is dominated by mesoscale eddies (Buckingham et al. 2016). Brannigan et al. (2015) find that submesoscale processes are not limited to frontal regions between mesoscale vortices but that submesoscale filaments are also found inside mesoscale eddies. Brannigan (2016) investigates the effect of the submesoscale processes inside a mesoscale anticyclone and finds additional mixed layer–thermocline exchange is driven by symmetric instability at submesoscale-permitting resolutions.

The aim of this paper is to examine the causes and effects of submesoscale instabilities in both a mesoscale cyclone and anticyclone in more detail. The structure of this paper is as follows: Section 2 describes the setup of the numerical simulations, methods, and the data-gathering process; section 3 presents the results of the idealized numerical experiments and develops a hypothesis to be tested using the glider observations; and section 4 presents glider observations of an anticyclone from the Tasman Sea, followed by a summary and discussion in section 5.

## 2. Experimental setup

### a. Submesoscale-permitting simulations

An idealized numerical model is used in Brannigan et al. (2015) to consider the seasonal cycle of a submesoscale-permitting flow. The model domain is doubly periodic in the horizontal plane and is an analog of an open-ocean region. Brannigan (2016) uses this simulation once it has spun up to study the development of submesoscale filaments inside a mesoscale anticyclonic eddy. The same experiments are analyzed in further detail in this manuscript.

The spinup of the model is described in Brannigan et al. (2015), and the numerical setup is described in detail in appendix A. To investigate the upwelling of nutrient-rich fluid from the thermocline by submesoscale instabilities, a mesoscale–eddy resolving simulation at 4-km horizontal grid resolution is interpolated to finer resolution grids of 2, 1, 0.5, and 0.25 km at a time when the mixed layer depth is approximately 35 m. The simulations are run for 27 days at all resolutions.

As for Brannigan (2016), the surface frictional boundary condition is a wind stress calculated relative to a zonal 10-m wind of 6.3 m s^{−1} that gives a mean surface stress of approximately 0.05 N m^{−2}. The development of submesoscale instabilities in the model anticyclones is not sensitive to variations in the orientation of the wind stress (Brannigan 2016). The zonal wind induces a mean Ekman transport to the south in this Northern Hemisphere configuration. A cooling of 75 W m^{−2} is applied to counter submesoscale restratification of the mixed layer (Boccaletti et al. 2007; Couvelard et al. 2015), and so the simulations presented here are best thought of as an early wintertime scenario when the mixed layer is deepening.

A numerical passive tracer release experiment is carried out with the initial condition of a passive tracer concentration *C* = 1 in the surface level of the model and *C* = 0 below. No restoring of the passive tracer is applied.

As in Brannigan et al. (2015), the mixed layer is defined here as the first depth where the temperature is more than 0.1°C below that of the mean temperature in the upper 9 m of the domain. The diagnostics of submesoscale flows through the seasonal cycle in Brannigan et al. (2015) show that this definition captures the different dynamical processes between the mixed layer and thermocline.

### b. Azimuthal and radial velocity

The azimuthal velocity in an eddy is used to understand the effect of the instability on the flow around the eddies and is taken to be , where **u**_{c} = (*u*_{c}, *υ*_{c}, *w*) is the velocity vector in Cartesian coordinates, and *θ* is the angle between the vector from the eddy center to a point and the *x* axis. The radial velocity . The eddy center is defined as the location of the pressure extremum.

### c. Azimuthal averaging

Azimuthal averages as functions of depth and distance *r* from the eddy center are taken around the eddies to examine the effect of varying the model resolution. To do this, model outputs for any function *h*(*x*, *y*, *z*) are averaged over points with radius between *r* − Δ*r* and *r* + Δ*r*, where Δ*r* is twice the horizontal grid spacing:

where *i*, *j*, *k* are the gridpoint *x*, *y*, and *z* indices; *x*_{ij} and *y*_{ij} are the *x* and *y* positions of the grid centers and *A*_{ij} the grid areas; and *z*_{k} is the vertical position of the *k*th grid level.

### d. Turbulent kinetic energy

An energetic analysis of submesoscale instabilities in the eddies is carried out below. Such an analysis requires a decomposition of the flow field into mean and perturbation components. In frontal scenarios, this decomposition can be done by averaging along a relatively straight section of the front (e.g., Capet et al. 2008b; Suzuki et al. 2016), where the mean flow is then the alongfront average, and the perturbation represents departures from that mean. The simplest approach for a vortex is to take the azimuthal average around a given radius of the eddy to be the mean flow (e.g., Smyth and McWilliams 1998). However, this approach is not readily applicable here as the mean mesoscale flow around the eddy is not axisymmetric, as shown in section 3b. In Brannigan (2016), this asymmetry is handled by interpolating a lower-resolution simulation to a finer-resolution grid to be the mean flow. This approach, however, is only valid in the initial days of the simulation.

We use here a more robust decomposition of the flow into mean and perturbation components with a low-pass spatial filter that defines the mean flow as the along-stream average around a portion of the eddy. More precisely, the mean flow is defined at a point by averaging around one-third of the eddy:

where *h* is any field defined around the eddy. The small-scale perturbation field is then defined as .

The derivation of the turbulent kinetic energy budget is set out in appendix B. The focus of the analysis in this paper is the terms that lead to the net production and dissipation of turbulent kinetic energy when averaged across the eddy. The dominant terms—all with units of meters squared per cubic second—are the vertical buoyancy flux with as the azimuthal-averaging operator in Eq. (1a), where the buoyancy anomaly and *g* is the gravitational acceleration, *ρ* is the potential density, and *ρ*_{0} is a reference potential density; the vertical shear production is , where the subscript denotes differentiation; the lateral shear production is ; and the vertical dissipation of turbulent kinetic energy is . The horizontal dissipation of turbulent kinetic energy is somewhat smaller than the vertical dissipation in the mixed layer of the eddies and so is excluded.

### e. Ekman buoyancy flux

The Ekman buoyancy flux is

where *τ* is the surface wind stress vector, **k** is the vertical unit vector, and is evaluated at the surface level. The Ekman buoyancy flux is a source of mean potential energy, where the Ekman transport leads to a steepening of isopycnals and a sink of mean potential energy where it flattens isopycnals.

Thomas (2005) shows that wind stress aligned with the geostrophic shear—known as a downfront wind—extracts potential vorticity from the ocean such that it leaves the ocean unstable to symmetric instability. On the other hand, a wind stress that opposes the geostrophic shear—known as an upfront wind—leads to an input of potential vorticity to the ocean. Downfront winds lead to positive values of the Ekman buoyancy flux, while upfront winds lead to negative values of the Ekman buoyancy flux.

### f. Mixed layer instability parameterization

Additional simulations are carried out at 4-km resolution using a parameterization for mixed layer eddies (Fox-Kemper et al. 2011). This parameterization was developed based on simulations using a baroclinic front with no curvature (Fox-Kemper et al. 2008), and so it is of interest to understand the extent to which it can reduce bias between resolutions when applied in mesoscale eddies where the curvature of the flow is also dynamically important.

The parameterization takes the form of a vector streamfunction **Ψ** defined by

where *C*_{e} is a nondimensional coefficient, Δ is the model grid spacing, *L*_{f} is an estimate of the typical frontal width, *H* is the mixed layer depth, is the mean lateral buoyancy gradient averaged over the depth of the mixed layer, *T* is a time scale of mixed layer instability, and *μ* is a vertical structure function. The length and time scales used here are *L*_{f} = 1 km, and the time scale *T* = 5.8 days. Initial tests with *C*_{e} = 0.07 made little impact on the mixed layer temperature differences between the simulation at 4- and 0.25-km resolution. Therefore, the results presented below use *C*_{e} = 1, and this goes much further in reducing differences in the mixed layer temperature profiles between the 4- and 0.25-km simulation. This higher value of the nondimensional coefficient is in line with that proposed by Bachman and Taylor (2016) in simulations of equilibrated mixed layer instabilities. The results presented below are for the simulation at 4 km with no mixed layer instability parameterization unless otherwise stated. The streamfunction **Ψ** generates an eddy velocity that advects the temperature field and passive tracers.

### g. Convective layer depth

The convective layer depth is the depth to which convective motions driven by surface buoyancy loss dominate the turbulent kinetic energy budget. Below this depth, symmetric instability is expected to be the dominant process in a region of negative potential vorticity (Taylor and Ferrari 2010). A quartic polynomial to predict the convective layer depth (Thomas et al. 2013) is

where *h*_{cld} is the convective layer depth. The depth scale *H*_{pv} is based on a potential vorticity criterion as the deepest depth where *q* < 10^{−9} s^{−3}, while *c* = 14 is an empirical constant, is the convective velocity scale, where *B*_{0} is the surface buoyancy flux; is the friction velocity, with *ρ* as the potential density; Δ*u*_{g} is the bulk difference in the geostrophic velocity across the low potential vorticity layer; and *ϕ* is the angle between the wind vector and the geostrophic shear.

### h. Potential vorticity

The potential vorticity is

where *f* = 10^{−4} s^{−1} is the Coriolis frequency, and ∇ is the gradient operator. Horizontal gradients in the vertical velocity are everywhere negligible at these resolutions and so are excluded from the calculation. Division of the Ertel potential vorticity by a varying potential density is neglected because of the Boussinesq assumption.

The vertical component of relative vorticity is , where subscripts denote differentiation. Following Thomas et al. (2013), the components of the potential vorticity are referred to in the text as the baroclinic component and the vertical component .

In the Northern Hemisphere, fluid with positive potential vorticity is stable to gravitational/symmetric/inertial instability, while fluid with negative potential vorticity is unstable (Thomas et al. 2013). The converse applies in the Southern Hemisphere because of the change in sign of *f*. This can be accounted for by considering *fq* as the parameter for determining stability rather than *q* alone (Hoskins 1974). However, we shall only consider the Northern Hemisphere where negative potential vorticity corresponds to instability.

There is some inconsistency in the terminology used in the literature. Hoskins (1974) defines any state where *fq* < 0 as being symmetrically unstable. A state with *fq* < 0 could be due to any or all of unstable stratification, anticyclonic lateral shear, or anticyclonic baroclinic component. In contrast, in the recent oceanic literature, symmetric instability refers specifically to the state where *fq* < 0 because of an anticyclonic baroclinic component (e.g., Thomas et al. 2013). We follow this latter convention.

### i. Glider data

The glider data from the Tasman Sea were first examined by Baird and Ridgway (2012), and a more detailed account of the data gathering and processing can be found there. The glider was deployed in late winter and the deployment continued until early spring 2010. The glider sampled over the upper 1000 m of the water column. The horizontal spacing of the glider profiles varied from about 2 km in weaker currents to up to 16 km when the glider became entrained in faster currents. The data were downloaded from the Integrated Marine Observing System in Australia. (The downloaded data file is IMOS_ANFOG_BCEOSTUV_20100809T003827Z_SG516_FV01_timeseries_END-20101210T104213Z_C-20120117T064827Z.nc.)

Absolute temperature and conservative salinity for the glider are calculated from the in situ data using the Gibbs Seawater Oceanographic Toolbox that uses the Thermodynamic Equation of Seawater—2010 (TEOS-10; IOC et al. 2010; McDougall and Barker 2011).

### j. Remote sensing data

The satellite altimetry data used in conjunction with the glider data are the delayed-time, global, all-satellite gridded product (dt_global_allsat_msla) by AVISO. The surface chlorophyll *a* images for the Tasman Sea are level-two data from the NASA MODIS *Aqua* instrument. (The images used for the Tasman Sea are A2010294041500.L2_LAC_OC.nc and A2010269042000.L2_LAC_OC.nc.)

## 3. Idealized simulations

### a. Initial state

At the time chosen for the interpolation from the 4-km run to the finer grids (Fig. 1a), there is a cold-core cyclone in the north of the domain at (120, 195) km and a warm-core anticyclone in the south at (95, 65) km. The temperature field in the anticyclone is relatively straightforward in that temperature decreases monotonically away from the core of the anticyclone. On the other hand, the temperature increases away from the core of the cyclone for approximately 35 km before decreasing again beyond this. The flow in this outer region of reversed radial temperature gradients nonetheless has cyclonic angular velocity.

The zonal wind stress induces an Ekman buoyancy flux wherever there is a meridional buoyancy gradient. In the anticyclone there are downfront winds and a positive Ekman buoyancy flux in the north of the eddy (Fig. 1b) with upfront winds and a negative Ekman buoyancy flux in the south of the eddy. As the meridional buoyancy gradients are reversed in the cyclone, the relative positions of the positive and negative Ekman buoyancy fluxes are reversed (Fig. 1b). The root-mean-square magnitude of the Ekman buoyancy flux in the anticyclone is 7 × 10^{−8} m^{−2} s^{−3}, which is almost twice as large as that in the cyclone. This difference reflects stronger lateral buoyancy gradients in the anticyclone. In addition, the reversal of the meridional buoyancy gradient in the cyclone means that there is a further region of positive Ekman buoyancy fluxes in the outer northern region of the cyclone near (120, 230) km.

The near-surface potential vorticity *q* (Fig. 1c) has an opposite distribution to the Ekman buoyancy fluxes. In both the cyclone and the anticyclone there is generally negative potential vorticity where there are downfront winds. The magnitude of the potential vorticity differences across the anticyclone is larger than those in the cyclone. The cyclone has negative potential vorticity both in the downfront region in the core of the eddy and in the downfront region on the northern side of the eddy.

The correlation between the Ekman buoyancy flux (Fig. 1b) and the potential vorticity (c) is −0.97. This high correlation shows that symmetric instability is unlikely to be active in redistributing potential vorticity in the 4-km simulation. The magnitude of this correlation drops from −0.5 to −0.3 from the third day of the simulation at 0.25-km resolution (not shown) as symmetric instability becomes active in redistributing potential vorticity (Brannigan 2016).

### b. Anticyclonic eddy

#### 1) Resolution dependence of instability

The degree to which submesoscale instabilities are present in the anticyclone centered at (95, 65) km on day 9 of the simulations can be understood by considering *ζ*^{z} at the surface (Fig. 2). Submesoscale processes lead to anomalies in *ζ*^{z} because of vortex stretching and tilting and advection in the presence of larger-scale vorticity gradients. At 4-km resolution, there is little evidence for submesoscale instabilities occurring in the anticyclone (Fig. 2a). However, the filaments of positive relative vorticity are in evidence in the anticyclone at 2-km resolution (Fig. 2b). With each subsequent increase in resolution (Figs. 2c–e), further filaments with positive relative vorticity values are present. The width of the filaments decreases with each increase in resolution, including the final increase to 250-m resolution. The anomalies in *ζ*^{z} on day 9 primarily have an azimuthal wave structure (Fig. 2e). This structure is in contrast to the anomalies in potential vorticity in days 2–4, where there is a distinct radial wave structure [cf. Figs. 1f and 2b in Brannigan (2016)].

#### 2) Phenomenology of submesoscale processes in the anticyclone

The azimuthal mean flow around the eddy is a broad jet (Fig. 3a) at 4-km resolution where submesoscale processes are not permitted. Even at 4-km resolution, however, this jet is not axisymmetric. Instead the jet is stronger on the northeastern side of the eddy and weaker on the southeastern side. This partly reflects the nonnegligible ellipticity of the anticyclone. The difference in baroclinic structure between the north and south of the eddy—as shown by the difference in near-surface potential vorticity (Fig. 1c)—may also play a role.

The pattern of stronger azimuthal flow on the northeastern side of the eddy is also found at 0.25-km resolution (Fig. 3b). However, at this finer resolution the initial broad jet has been sharpened into a series of narrower and more intense submesoscale jets and occasionally vortices. These submesoscale jets do not necessarily follow pathways of constant radius around the eddy but can instead wrap in toward the core, such as the jet at (115, 68) km in (Fig. 3b). These jets can thus redistribute momentum radially within the eddy. In time, the cyclonic potential vorticity filaments in the eddies can wrap-up to form submesoscale cyclonic vortices [e.g., at (105, 85) km in Fig. 3b] that are transported anticyclonically by the mesoscale eddy.

The temperature anomalies (Fig. 3c) at 4-km resolution also show departures from axisymmetry, whereby the warmest waters in the anticyclone are slightly to the south of the eddy center because of the southward Ekman transport. When submesoscale processes are permitted (Fig. 3d) the temperature field has more small-scale structure. The jet in azimuthal velocity in Fig. 3b at (115, 68) km is a cold filament intruding into the eddy core (Fig. 3d), while the cyclonic eddy at (105, 85) km is a cold anomaly.

The sharp gradients in momentum and temperature around the cold filaments suggest that their dynamics may be similar to those of cold filaments studied in frontal systems. In those frontal systems, the cold filaments undergo rapid frontogenesis caused by ageostrophic secondary circulations driven by flow convergence on the eddy scale and turbulent thermal wind (McWilliams et al. 2009; Gula et al. 2014; McWilliams et al. 2015; Wenegrat and McPhaden 2016), while losing kinetic energy because of lateral shear instabilities (Gula et al. 2014). This is considered further in section 3b(4) below.

#### 3) Active tracer transport in the anticyclone

The net effect on the anticyclone of permitting submesoscale instabilities is considered by comparing azimuthally averaged fields at 4- (both with and without the mixed layer instability parameterization) and 2-km resolution with the simulation at 0.25-km resolution on day 18.

In the mixed layer of the anticyclone, the 4-km simulation (Fig. 4a) is warmer than the 0.25-km simulation in the eddy core but cooler farther out in the eddy around a radius of 30 km. The simulation at 4 km is also colder below the mixed layer. In the 4-km simulation with the mixed layer instability parameterization (Fig. 4d) the warm bias in the core of the eddy is substantially reduced. However, the cold bias farther out in the eddy is reduced to a much smaller degree in the parameterized simulation. The cold bias below the mixed layer along isopycnals that outcrop is unchanged by adding the parameterization that applies only in the mixed layer (Figs. 4a,d). However, in the simulation at 2-km resolution all of the warm and cold biases relative to the simulation at 0.25-km resolution are significantly reduced (Fig. 4g).

There is a substantial difference in the azimuthal velocity profile of the 4-km simulation—both with and without the mixed layer instability parameterization—relative to the 0.25-km simulation (Figs. 4b,e,j) as the azimuthal jet shifts inward at finer resolution. However, the inward shift of the jet does occur at 2-km resolution in the anticyclone (Figs. 4h,j). The difference in the azimuthal velocity profile even with the parameterization shows that the changes in velocity are not due to changes in the buoyancy field but are instead due to instabilities that lead to the redistribution of momentum (Marshall et al. 2012). The relative passive tracer concentrations in the anticyclone are considered in section 3b(5) below.

In summary, the main effects of permitting submesoscale processes in the anticyclone across this range of resolution happen from the change in resolution from 4 to 2 km, with smaller changes when the resolution is further refined from 2 to 0.25 km. The simulation at 2-km resolution is the coarsest resolution where the submesoscale filaments in surface relative vorticity appear (Fig. 2), and so the presence of such filaments may be a good qualitative indicator of active submesoscale processes. The resolution dependence for submesoscale processes in the cyclone is quite different, as shown in section 3c.

#### 4) Energetics of instabilities in the anticyclone

The changes in the buoyancy and azimuthal momentum within the anticyclone at finer resolution [section 3b(3)] show that submesoscale processes have rectified effects on the mesoscale eddy. An energetic analysis shows the spatial structure of energy production and dissipation by submesoscale processes within the anticyclone at 0.25-km resolution (Fig. 5). The analysis is carried out using snapshots output at 12-h intervals over the first 18 days of the simulations. We caution that the sampling frequency means that there is a degree of aliasing in the eddy statistics.

The horizontal structure of energy transfer to turbulent kinetic energy in the mixed layer is shown by averaging in time and depth over the upper 60 m (Figs. 5a–d). The vertical buoyancy flux (Fig. 5a) and vertical shear production (Fig. 5b) are largest on the downwind side of the eddy, where fluid particles have been exposed to downfront winds for the longest period and the lowest values and strongest gradients of potential vorticity are found (Brannigan 2016). The vertical buoyancy flux is large all around the eddy (Fig. 5a), while the largest values of vertical shear production are concentrated on the downwind side of the eddy (Fig. 5b).

The magnitude of lateral shear production is largest on the downwind side of the eddy (Fig. 5c), where it is a net sink of kinetic energy from the turbulent flow. Locally, these negative values are associated with the filaments that penetrate radially inward (Fig. 3b) and lead to the shifts in the mean azimuthal jet. The dissipation of turbulent kinetic energy by vertical viscosity is also largest on the downwind side of the eddy where gradients are strongest. The viscous dissipation of turbulent kinetic energy caused by horizontal dissipation has a similar spatial distribution to the vertical dissipation, though it is weaker (not shown).

The vertical profile of the vertical buoyancy flux (Fig. 5e) shows that it is positive throughout the 40-m-deep mean mixed layer and negative in the entrainment layer below this. The vertical shear production (Fig. 5f) is negative in the upper 15 m and positive deeper in the mixed layer. The mean convective layer depth around the anticyclone is 15 m, consistent with the layer below this being the region of positive vertical shear production (Taylor and Ferrari 2010). The lateral shear production (Fig. 5g) is largest in the upper 25 m. The dissipation of turbulent kinetic energy is largest in the middle of the mean mixed layer, where the vertical viscosity is highest.

Comparing the time- and azimuthally averaged fields over the range of resolutions shows that the major change in the vertical buoyancy flux occurs with the change in resolution from 4 to 2 km (Fig. 6a). The largest difference in vertical shear production, however, occurs with the increase in resolution from 2 to 1 km, and the positive values of vertical shear production increase throughout the lower part of the mixed layer as the resolution is further refined (Fig. 6c). The negative values of the lateral shear production (Fig. 6e) are large at all resolutions finer than 4 km, though with resolution dependence breaking down for the 0.25-km simulation. The overall pattern of larger lateral shear production at the finer resolutions is consistent with the changes in the mean azimuthal jet that are found once the resolution is refined from 4 km (Fig. 4j). The vertical dissipation doubles between 4- and 2-km resolution but then is similar to the 2-km simulation at all the finer resolutions (Fig. 6g).

The energetic analysis suggests that at the finer-resolution simulations the anticyclone is unstable to baroclinic and vertical shear-driven instabilities but is “barotropically” stable.^{1} The increasing values of vertical shear production as the resolution is refined beyond 1 km imply that vertical shear-driven instability is strong in the lower part of the mixed layer and leads to exchange between the mixed layer and thermocline (Brannigan 2016). The positive values of the vertical buoyancy flux imply that baroclinic instability grows through the full depth of the mixed layer. These instabilities grow fastest on the downwind side of the anticyclone, though baroclinic instability occurs the whole way around the eddy.

The sink of turbulent kinetic energy caused by lateral shear instability is due to the cold submesoscale filaments that wrap radially inward in the anticyclone such as at (115, 68) km (Figs. 3b,d). In effect, the growing baroclinic waves on the downwind side of the eddy gain energy from positive vertical buoyancy fluxes but are oriented such that they lead to a transfer of energy from turbulent kinetic energy to the mean kinetic energy of the vortex. The energy loss of these cold submesoscale filaments caused by lateral shear instability is similar to that found for cold filaments in simulations of the Gulf Stream (Gula et al. 2014).

This energetic analysis is carried out during the period when the simulation at 0.25-km resolution is adjusting to the change in resolution from 4 to 0.25 km. We carry out the same analysis on an anticyclone at 0.5-km grid resolution during the winter period of year five from the seasonal cycle simulations of Brannigan et al. (2015), where there is no adjustment to changes in resolution (not shown). In this case, a similar spatial pattern holds for the vertical buoyancy fluxes, vertical shear production, and vertical dissipation. The lateral shear production is more variable with large negative values found on the downwind side of the eddy but with large positive and negative values found elsewhere around the eddy. As such the spatial pattern of lateral shear instability (Figs. 5c,g) partially reflects an adjustment to the finer resolution grid.

#### 5) Passive tracer subduction in the anticyclone

At finer resolutions there is subduction of the passive tracer into the thermocline in the anticyclone in the 2- and 0.25-km simulations that does not occur in the simulations at 4-km resolution (Figs. 4c,i). Applying the mixed layer instability parameterization to the 4-km simulation does not affect the magnitude of this difference between the simulations (Fig. 4f) in the thermocline. The difference in effect between the active buoyancy tracer—where the bias between the simulations is reduced—and the passive buoyancy tracer—where the bias is not reduced—is anticipated and arises because the parameterized eddy streamfunction has been determined using buoyancy alone (Bachman and Fox-Kemper 2013).

For the subducted fluid in the thermocline, however, there is a change in the relationship between the tracer concentration and the components of potential vorticity. The regions of high tracer concentration in the anticyclone (Fig. 7a)—shown by the black contours in (Figs. 7b–d)—have low baroclinic potential vorticity (Fig. 7b) only at their boundaries. Instead, the regions of high tracer concentration have low potential vorticity because of strong anticyclonic relative vorticity (Fig. 7c) and weak stratification (Fig. 7d).

These results for a baroclinic eddy are in agreement with the results of Thomas (2008), who shows that parcels subducted from fronts have low baroclinic potential vorticity in the mixed layer but low vertical potential vorticity in the thermocline because of vortex tilting.

Along isopycnals that outcrop to the surface in the anticyclone there is higher tracer concentration in the mixed layer and lower tracer concentration in the thermocline at 4-km resolution (Fig. 8a). The higher values of the passive tracer in the thermocline at finer resolution (negative values in Fig. 8a) are associated with low potential vorticity and a doming of the isopycnals (Fig. 8c), for example, at finer resolution near (78 km, −50 m) in Fig. 8a. As the low potential vorticity lenses are found in the relatively adiabatic thermocline, these structures persist for some time before being eroded by mixed layer deepening. There are also regions of relatively low passive tracer concentration in the mixed layer at 0.25-km resolution, for example, at (82 km, −50 m) or (89 km, −40 m) in Fig. 8a that are associated with upwelling plumes of high potential vorticity (Fig. 8c). Brannigan (2016) concludes that these narrow upwelling plumes are due to entrainment of high potential vorticity fluid from the thermocline by symmetric instability. In light of the further analysis of the tracer fields and energetics here, it is likely more accurate to describe the overall vertical exchange arising as a combination of baroclinic instability and vertical shear-driven instability, though a detailed deconstruction is not possible. The presence of such upwelling plumes would lead to countergradient vertical fluxes for tracers with high concentrations nearer the surface, as suggested by Smith et al. (2016).

The numerical results provide a clear hypothesis for observational data. The simulations predict that wind-driven extraction of potential vorticity at the surface of anticyclonic eddies can lead to submesoscale instabilities and the consequent transport of low potential vorticity fluid into the thermocline. The low potential vorticity fluid should have domed isopycnals, strong anticyclonic relative vorticity, and other indicators of a recent surface origin such as high dissolved oxygen concentrations.

### c. Cyclonic eddy

#### 1) Submesoscale instabilities in a cyclonic eddy

Submesoscale filaments do grow in the core of the cyclone from day 8 in the simulation at 0.25-km resolution (Figs. 9c–h). Filament growth is slower in the cyclone that the anticyclone where similar filaments emerge after just 3 days of the experiment (Brannigan 2016). The filament growth occurs first on the downwind side of the cyclone [near (125, 190) km in Figs. 9c,d], and thereafter filaments are found further around the cyclone (Figs. 9e–h). The initial filaments have a wavenumber primarily in the azimuthal direction, as opposed to the anticyclone where the initial filaments had a wavenumber primarily in the radial direction (Brannigan 2016). Faster growing submesoscale instabilities are also evident in the region with sharp lateral buoyancy gradients in the northern part of the cyclone (Figs. 9a,b).

#### 2) Active and passive tracer transport in the cyclone

As for the anticyclone, the net effect of permitting submesoscale motions in the cyclone is considered by comparing the azimuthally averaged profiles of temperature, momentum, and passive tracer concentration on day 18 (Fig. 10). In contrast to the anticyclone, there is a similar distribution in temperature at resolutions of 4 (Fig. 10a) and 2 km (Fig. 10g), and these are different from the simulation at 0.25-km resolution. These temperature differences with respect to the 0.25-km resolution simulation are somewhat reduced, however, by applying the mixed layer instability parameterization to the 4-km simulation (Fig. 10d).

As for the anticyclone, there is an inward shift of the jet around the cyclone in the inner 10 km of the eddy (Fig. 10j) in the simulation at 0.25-km resolution. This inward shift does not occur in the simulations at 4 km, with and without the mixed layer instability parameterization and the simulation at 2 km (Fig. 10j). The existence of these differences between the simulations is in agreement with the lack of filamentation visible in surface relative vorticity in the cyclone at 4- and 2-km resolution (Figs. 2a,b) compared to 0.25-km resolution (Fig. 2e).

The azimuthally averaged tracer distribution (Fig. 10f) shows that additional tracer transport into the thermocline occurs at 0.25-km resolution. In this case, the bias is increased in the simulation at 4-km resolution with a mixed layer instability parameterization.

Considering vertical sections through the cyclone, there is a much smaller difference in the passive tracer distribution between the simulations at 4- and 0.25-km resolution compared to the anticyclone (Fig. 8b). In addition, regions of negative potential vorticity are more limited in extent with smaller magnitudes of negative values (Fig. 8d). The smaller anomalies with respect to the finer resolution runs are indicative of generally weaker submesoscale instabilities in the cyclone compared to the anticyclone.

#### 3) Energetics of instabilities in the cyclone

An energetic analysis shows the spatial structure of energy production and dissipation by submesoscale processes within the cyclone at 0.25-km resolution. The analysis is carried out using snapshot outputs at 12-h intervals over the first 18 days of the simulations (Fig. 11).

The vertical buoyancy flux is positive in the core of the eddy and is largest on the downwind side of the eddy (Fig. 11a). In contrast to the anticyclone, however, the magnitude of the vertical shear production and lateral shear production are much weaker than the vertical buoyancy fluxes (Figs. 11b,c). The vertical dissipation (Fig. 11d) is elevated on the downwind side of the eddy, as in the anticyclone.

In the vertical, the buoyancy fluxes are positive through the mixed layer in the cyclone (Fig. 11e). The vertical and lateral shear productions are low at all depths (Figs. 11f,g). The vertical dissipation of turbulent kinetic energy (Fig. 11h) is again highest in the center of the mixed layer.

Comparing the energetics across resolutions, there is a stark contrast with the anticyclone as the mean vertical buoyancy fluxes do not begin to grow until the grid is refined to 0.5 km, and there is a further large change when the grid is refined to 0.25 km (Fig. 6b). Vertical shear production is low compared to the anticyclone at all resolutions (Fig. 6d). Lateral shear production is high only at 0.25 km (Fig. 6f), while the vertical dissipation is similar at all resolutions finer than 4 km (Fig. 6h).

This energetic analysis suggests that restratifying baroclinic instabilities are active in the core of the cyclone. Despite the presence of regions of negative potential vorticity in the cyclone, it appears that symmetric instability is effectively absent from the cyclone in these simulations. This energetic analysis is consistent with the qualitative picture that emerges from the surface relative vorticity in Fig. 2, where submesoscale anomalies in relative vorticity in the cyclone begin to emerge at 1-km resolution (Fig. 11c), and large differences are apparent as the grid is further refined toward 0.25 km (Figs. 11d,e). One consequence of these weaker instabilities is that mixed layer–thermocline exchange is weaker in the cyclone than the cyclone. This cyclone–anticyclone asymmetry is discussed further in section 5.

## 4. Observations

The numerical results set out in Brannigan (2016) and in section 3 also provide a clear hypothesis that can be tested observationally. The model predicts that lenses of low potential vorticity fluid may be found in the thermocline of anticyclonic eddies because of submesoscale instabilities in the mixed layer. An initial test of this hypothesis is performed here using glider observations of an anticyclonic eddy in the Tasman Sea.

The glider deployment took place in the Tasman Sea during the austral winter to spring transition from August to December 2010, as described in Baird and Ridgway (2012). The deployment sampled an anticyclonic mesoscale eddy with diameter of approximately 150 km (Fig. 12a). The nature of the sampling strategy varied over the course of the transect. The glider track (Fig. 12b) shows that when the glider crossed the periphery of the eddy in late September and early October, the sampling strategy can be thought of as a transect. However, when the glider then sampled the core of the eddy, the sampling strategy was quasi Lagrangian (Fig. 12c). The stronger flows in the periphery of the eddy mean that the profiles are separated by up to 16 km in the eddy periphery, whereas profile separations of approximately 2 km are more typical in the eddy core.

The anticyclone has lower chlorophyll *a* concentrations than adjacent regions (Figs. 12b,c) because of the late winter mixed layer being much deeper in the anticyclone compared to around the eddy. This means that restratification (Mahadevan et al. 2010) or suppression of vertical mixing processes (Huisman et al. 1999; Smith et al. 2016; Taylor and Ferrari 2011; Taylor 2016) are likely to be dominant in setting the chlorophyll *a* patterns in the anticyclone rather than the nutrient supply mechanism for an oligotrophic region considered in Brannigan (2016). However, submesoscale filamentary features are present in the interior of the mesoscale anticyclone (Figs. 12b,c).

The in situ data for the portion of the deployment when the glider crossed the periphery of the eddy show that the stratification in the eddy is driven by temperature gradients that are partially compensated by gradients in salinity (near 1400 km in Figs. 13a,c). The density surfaces that are found at depths of hundreds of meters below the core of the anticyclone outcrop around its periphery in this region.

Water with a higher oxygen concentration than the surrounding water on the same density surface is found around 1400 km of the transect at −300 < *z* < −100 m. This water mass is referred to hereinafter as the “high oxygen water.” The density contours show that this water is not in the mixed layer, though its high oxygen values suggest that it must recently have been in the mixed layer. The fluorescence values (Fig. 13d) show that this water also has higher fluorescence than the water around it, consistent with the water having recently been in the near-surface euphotic zone.

In these late winter observations, the fluorescence also provides a useful indicator for the depth to which surface mixing is penetrating. In this case, the high oxygen water is found 200 m below the high fluorescence layer near the surface, showing that the high oxygen water is not subject to surface-driven mixing at this point. Detailed inspection of the data shows that the high oxygen water has stratification an order of magnitude weaker than the surrounding water suggesting that it has low potential vorticity (not shown).

As the glider continued to sample in the core of the anticyclone, repeated observations were made of a lens of water with anomalously high temperature, salinity, and oxygen concentrations between −700 < *z* < −400 m (1600–2000 km in Figs. 13a,b,c). The oxygen concentration in the lens is higher than the oxygen concentration in the deep mixed layer above it, suggesting that the water originated in a shallower mixed layer. There is a clear doming of the isopycnals around the lens, showing that it is a low stratification and presumably low potential vorticity feature.

The combination of high oxygen and low potential vorticity suggests that the water in the lens has a mixed layer origin. Following the isopycnals where the lens is found back along the transect shows that these isopycnals are in the upper 200 m but do not outcrop. Baird and Ridgway (2012) perform a detailed hydrographic analysis of the lenses. They show that the water in the lens does have a local surface origin as the temperature–salinity properties show that water with these properties is only found in late winter in the shallow Bass Strait region between the Australian mainland and Tasmania (for reference in Fig. 12, Bass Strait is near the southwest corner of the images). Water from Bass Strait is known to flow northward along the shelf toward where the eddy was observed (Baird and Ridgway 2012). Therefore, these density surfaces do outcrop in the region. Baird and Ridgway (2012) also observe a number of other anticyclonic eddies approximately 700 km to the south in the direction of the poleward mean flow in this western boundary current region. They find more low stratification lenses in the thermocline of all of the anticyclones that they sample with similar hydrographic properties indicating a Bass Strait origin.

While the observational dataset does not allow definitive conclusions, we conjecture that the high oxygen filament and lenses have been created by subduction of mixed layer fluid because of submesoscale instabilities in the mixed layer such as baroclinic instability (Spall 1995) or symmetric instability (Brannigan 2016). There are of course alternative hypotheses for the generation of the high oxygen water. For example, frontogenetic or geostrophic adjustment processes (Hoskins and Bretherton 1972; Shakespeare and Taylor 2013) can lead to tracer subduction (Macvean and Woods 1980). A further hypothesis is that the filament arose due to the straining of tracer fields along isopycnals by the mesoscale quasigeostrophic-type flow (Smith and Ferrari 2009). However, the slope of the high oxygen water is much shallower than the *N*/*f* slope predicted by quasigeostrophic theory.

## 5. Discussion

The results presented in this paper show that mesoscale baroclinic eddies are unstable to a range of submesoscale instabilities. In this range of resolutions, the production of turbulent kinetic energy is stronger in anticyclones than cyclones. The strongest eddy production occurs in the downwind quadrant of both cyclones and anticyclones. The submesoscale instabilities also lead to a radial redistribution of momentum compared to the simulation at 4 km. This radial redistribution of momentum occurs from 2-km resolution in the anticyclone but only with 0.25-km resolution in the cyclone. The idealized simulations show that these submesoscale instabilities can lead to the subduction of low potential vorticity fluid in the thermocline of anticyclonic eddies with much less subduction in cyclonic eddies. Limited glider observations from the Tasman Sea show the presence of such low potential vorticity fluid in the thermocline of mesoscale anticyclones.

The submesoscale instabilities are different in the cyclones and anticyclones of the idealized simulations across this range of resolutions. In the anticyclone, the submesoscale flows are energized by vertical buoyancy fluxes and vertical shear production. However, in the cyclone only vertical buoyancy fluxes are a significant source of turbulent kinetic energy. Furthermore, in the anticyclone the simulations at 2 and 0.25 km are similar in many respects, whereas in the cyclone the submesoscale processes are much stronger at 0.25-km resolution compared to the coarser-resolution simulations.

The differences between the cyclone and anticyclone may be because the relative vorticity and angular velocity are stabilizing in the cyclone compared to the anticyclone. To illustrate this we derive in appendix C an analytical expression for the linear growth rate *σ* of symmetric disturbances in a baroclinic vortex:

where growing modes occur for real, positive values of *σ*, *M*^{4} = |∇_{h}*b*|^{2}, *N*^{2} = *b*_{z}, Ω = *V*/*r*_{0} is the angular velocity at the reference radius *r*_{0}, *k* is the radial wavenumber, and *m* is the vertical wavenumber. This expression generalizes two known limit cases. For vortices with no lateral buoyancy gradients or stratification where *M*^{2} = *N*^{2} = 0 the generalized Rayleigh criterion is recovered where (Kloosterziel and van Heijst 1991), and growing modes occur when . On the other hand, if we neglect curvature by letting *r*_{0} → ∞, we recover the corresponding expression for baroclinic fronts

derived by Bachman and Taylor (2014). The expression for *σ* in Eq. (7) shows that the fastest growing mode is expected to be aligned with the slope of the isopycnals, as for a baroclinic front (Bachman and Taylor 2014). The expression also shows that for a cyclone and anticyclone with similar balanced Richardson number, the growth rate will be larger in an anticyclonic vortex because of the second term on the right-hand side of Eq. (7). Equation (7) shows that for a given vertical and lateral buoyancy gradient, the fastest-growing mode moves to a longer wavelength as the mixed layer deepens. As such, at a given grid resolution these symmetric modes may be permitted during the period of deepest mixed layers but not during periods of shallower mixed layers. In these simulations, the lateral buoyancy gradient is about 20% stronger in the anticyclone with similar vertical stratification and mixed layer depth. This means that symmetric instability has a longer wavelength in the anticyclone and so the instability can grow faster over this range of resolutions.

Unfortunately, Eq. (7) only accounts for the growth of symmetric modes, while the results of this study indicate that nonsymmetric baroclinic instabilities are a dominant component of the instabilities in both cyclones and anticyclones. We have not found any comparable expression to Eq. (7) for these nonsymmetric modes as these appear to be subject to higher-order dynamics. Lahaye and Zeitlin (2015) study the stability of a two-layer baroclinic anticyclone and find that it is unstable to symmetric-, baroclinic-, and barotropic-type instabilities. They note that solving numerically for the linear stability of baroclinic vortices with continuous stratification is much more challenging, as the problem becomes nonseparable in the vertical and radial directions. The results here indicate a further complication compared to the large body of research of vortex stability (e.g., Kloosterziel and van Heijst 1991; Billant and Gallaire 2005; Lahaye and Zeitlin 2015; Lazar et al. 2013; Smyth and McWilliams 1998) in that the stability of the vortices varies in the azimuthal direction. This variation occurs because of the opposing effect of the Ekman buoyancy flux on either side of a vortex. A full understanding of the stability of vortices will require accounting for the effect of this azimuthal variation on various modes.

The results presented here show that submesoscale instabilities are a pathway to dissipate the kinetic energy of the mesoscale eddy field. The mean dissipation of turbulent kinetic energy in the mixed layer of the anticyclone is approximately 5 × 10^{−9} m^{2} s^{−3}, while the mean kinetic energy in the mixed layer is of order 0.1 m^{2} s^{−2}. This gives a notional spindown time scale of about 200 days caused by the dissipation of kinetic energy in the mixed layer by submesoscale instability. If this is an accurate assessment of the magnitude of kinetic energy dissipation by submesoscale instabilities, it would mean that submesoscale instabilities play a relatively small role in the kinetic energy budget of mesoscale eddies, particularly eddies that extend deeper with only a fraction of their volume in the mixed layer. There is also likely to be a strong seasonal variation to the dissipative effect of submesoscale instabilities, as submesoscale instabilities are likely to be more prevalent and so more dissipative in wintertime conditions (Brannigan et al. 2015; Callies et al. 2015; Mensa et al. 2013; Thompson et al. 2016; Buckingham et al. 2016).

The glider and remote sensing observations show that submesoscale processes are occurring inside a mesoscale eddy. The analysis indicates that drawing more definitive conclusions on the dominant processes in eddies likely requires multiplatform in situ observations. The results from the idealized simulations suggest that indicators of submesoscale processes such as increased viscous dissipation (D’Asaro et al. 2011) are most likely to be observed on the downwind side of the eddy rather than on the downfront side of the eddy.

Classical papers on subduction in the North Atlantic find subduction driven by frictional (also referred to as mechanical) fluxes of potential vorticity out of the ocean to be negligible in the large scale (Marshall et al. 1993). In the light of the results of Thomas (2005), Maze et al. (2013) consider the potential role of frictional fluxes of potential vorticity in creating “mode water” in an eddy-permitting simulation of the North Atlantic. The authors there also find that the frictional fluxes averaged over a large area of outcropping isopycnals are an order of magnitude lower than the diabatic fluxes. We note that frictional flux averaged over a vortex would average to zero in a steady state, yet we find that a net subduction of the passive tracer occurs (Fig. 4, right panels). This net subduction arises because submesoscale instabilities lead to a flux of potential vorticity across the base of the mixed layer (Thomas 2005; Taylor and Ferrari 2009; Brannigan 2016) where negative potential vorticity occurs, but no such advective response develops where positive potential vorticity occurs. This highlights the difference between the net fluxes of potential vorticity into an isopycnal layer considered by Maze et al. (2013) and the net transport of tracer along the isopycnal. An estimate of the global tracer transport caused by this process will be considered in future work.

## Acknowledgments

This work forms part of the OSMOSIS project funded by the Natural Environment Research Council grants NE/I019921/1 (Oxford University), NE/I01993X/1 (National Oceanography Centre), NE/I019999/1 (University of Southampton) and NE/I019905/1 (University of East Anglia). L. Brannigan also received funding from the Wenner-Grenn Foundation. The authors thank Andrew Thompson, James McWilliams, and two anonymous reviewers for their insightful comments on the manuscript. The model output is available on request from the authors. The surface chlorophyll *a* images are produced by the NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group; (2014) MODIS-*Aqua* Ocean Color Data; NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group (https://oceancolor.gsfc.nasa.gov/data/aqua/).

### APPENDIX A

#### Numerical Configuration

The MITgcm (Marshall et al. 1997) is employed in hydrostatic mode with a doubly periodic domain of horizontal length of 256 km and a depth of 3700 m. There are 200 vertical levels with vertical spacing increasing from 3 m in the upper 90 m to approximately 30 m at depth. Further experiments in nonhydrostatic mode show no difference from the hydrostatic simulations at these resolutions.

A linear equation of state in temperature is employed with a thermal expansion coefficient of *α* = 2 × 10^{−4} K^{−1}. Horizontal viscous and diffusive coefficients are biharmonic (Griffies and Hallberg 2000) with the Smagorinsky (1963) scheme used to set the horizontal viscous coefficient to preserve the small-scale structures of interest (Ilicak et al. 2012). The biharmonic Smagorinsky viscous coefficient is 3, the Laplacian vertical viscous and diffusive coefficients are 4 × 10^{−5} m^{2} s^{−1}, and the biharmonic horizontal diffusive coefficient is 10^{5} m^{4} s^{−1}. A quadratic bottom drag with nondimensional drag coefficient of 3 × 10^{−3} is imposed to dissipate momentum. The flux-limited version of the Prather (1986) scheme is used to advect temperature and the passive tracer. Hill et al. (2012) show that the diffusion with the Prather scheme is of similar order to that estimated from dye release in the ocean interior. Momentum advection is based on the model’s default centered second-order scheme. Surface boundary layer processes are parameterized with the K-profile parameterization (KPP; Large et al. 1994).

### APPENDIX B

#### Turbulent Kinetic Energy Budget in Polar Coordinates

##### a. Preliminaries on cylindrical coordinates

We convert from Cartesian coordinates (*x*, *y*, *z*) to cylindrical coordinates (*r*, *θ*, *z*). In terms of the horizontal Cartesian coordinates,

where the inverse tangent is appropriately defined for the particular quadrant.

##### b. Boussinesq equations in cylindrical coordinates

The radial momentum balance is

where *u* is the radial velocity; *υ* is the azimuthal velocity; *π* = *p*/*ρ* is the dynamic pressure, where *p* is pressure; *A*_{υ} is the vertical viscous coefficient set by the KPP scheme; and *A*_{h} is the biharmonic Smagorinsky viscosity.

The azimuthal momentum balance is

The vertical momentum balance in hydrostatic balance is

where *g* is the gravitational acceleration, and *ρ* is the potential density.

We define the material derivative as

##### c. Large-scale mean flow

The large-scale mean flow is defined using the spatial filter in Eq. (2) in the main text. As the vortices in the simulations are slightly elliptic, this large-scale mean flow has both azimuthal and radial components that vary in *r*, *θ*, and *z*.

The mean azimuthal flow is anticipated to be in gradient wind and hydrostatic balance with

and, to approximately satisfy the balance relationship,

Following Smyth and McWilliams (1998), to simplify the equations we define

as the angular velocity of the mean azimuthal flow and

as the relative vorticity of the mean azimuthal flow.

The perturbation radial velocity is *u*′ = *u* − *U* and similarly for the other variables.

##### d. Linear perturbation equations around the mean flow

We linearize the equations around the large-scale mean flow. This linearization means that the triple correlation terms will be neglected in the turbulent kinetic energy budget. As the triple correlation terms lead primarily to transport of turbulent kinetic energy rather than its production or dissipation they are not the focus of this study. The radial momentum balance of the perturbed flow is

The azimuthal momentum balance of the perturbed flow is

The vertical momentum balance of the perturbed flow is

where the buoyancy anomaly , where is the large-scale potential density field.

We define a modified material derivative for this linearized flow

and so the radial momentum balance is

while the azimuthal momentum balance is

where Eq. (B9) has been used.

##### e. Linearized turbulent kinetic energy budget

We take the dot product of (*u*′, *υ*′, *w*′) with the respective momentum balances to form the turbulent kinetic energy budget. For clarity these products are taken for each balance before they are combined:

and

Letting *K* = 0.5(*u*′^{2} + *υ*′^{2}) be the kinetic energy per unit mass of the perturbations, combining the components and taking the azimuthal average gives

where we use the equality to rewrite the lateral shear production term. The azimuthal averages include the gradients of the large-scale flow components, as these terms vary slowly in the azimuthal direction.

The terms in the first row on the right-hand side of (B17) are the production of turbulent kinetic energy from the along-flow gradient of the mean flow, the terms in the second row are the production of turbulent kinetic energy from the shear normal to the mean flow, the terms in the third row are the production of turbulent kinetic energy from the vertical shear of the mean flow, and the term in the fourth row is the production of turbulent kinetic energy from vertical buoyancy fluxes. The terms in the fifth row are the pressure work, the terms in the sixth row are the viscous transport of turbulent kinetic energy, and the terms in the last row are the negative definite dissipation of turbulent kinetic energy by vertical and horizontal viscosity.

### APPENDIX C

#### WKB Normal-Mode Solution

We look for solutions that vary in the radial direction only and so neglect any terms with ∂/∂*θ* from the linearized equations in appendix B [(B10) and (B11)]. The analysis is similar to that carried out by Bachman and Taylor (2014) for a baroclinic front. The region of interest is taken to be a radius *r*_{0}, which can be thought of as representing the conditions in the part of the vortex where lateral gradients in buoyancy and momentum are strongest. The model uses a WKB approach. As such, the modes considered are assumed to have high radial and vertical wavenumber and are driven by the local dynamics without boundary conditions.

We assume an axisymmetric and inviscid flow for simplicity and so the system of equations is

We define and and substitute these into (C1d). We define the absolute angular velocity *A* = (2Ω + *f*). In cyclogeostrophic and hydrostatic balance, we have the cyclogeostrophic thermal wind relation ∂*V*/∂*z* = *A*^{−1}*M*^{2}, and we substitute this into (C1b).

We look for normal-mode solutions of the form . We omit the ^ and so arrive at

From the continuity equation (C2e), we have

Inserting this expression for *w* into Eqs. (C2b) and (C2d) and using the hydrostatic relation (C2c) in (C2d), we get

To eliminate *υ*, we multiply Eq. (C4a) by *σ* and Eq. (C4b) by *A*. We also use Eq. (C4c) in (C4a) to eliminate *π*:

which we solve as an expression for *σ*:

Following Bachman and Taylor (2014), we can rearrange Eq. (C6) as

The final two terms drop out for modes aligned with the isopycnals, where *k*/*m* = *M*^{2}/*N*^{2}.

## REFERENCES

*Ocean Modeling in an Eddying Regime*,

*Geophys. Monogr.*, Vol. 177, Amer. Geophys. Union, 17–38.

## Footnotes

*Publisher’s Note:* This article was revised on 20 December 2017 to include a missing funding acknowledgement in the Acknowledgments section.

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

Barotropic in this case refers to the instability that arises due to an inflection point in the lateral momentum gradients that can occur in its simplest form in a barotropic domain with no density variations. We stress that the term barotropic here does not mean that the instability in this case has a barotropic structure in the vertical; indeed, the instability has a highly baroclinic structure (Fig. 5g).