Cross-isobath transport of Upper Circumpolar Deep Water (UCDW) provides a major source of heat to the Antarctic continental shelves. Adaptive sampling with a Slocum glider revealed that the UCDW regularly intrudes onto the western Antarctic Peninsula shelf within mesoscale eddies, and a linear stability analysis of the shelf-break current upstream confirmed eddy length scales and vertical structure are consistent with the baroclinic instability of the current. The properties of the most unstable mode are insensitive to current orientation but sensitive to bottom slope in accordance with modified Eady theory. Once on the shelf, the eddies’ core properties mix with ambient shelf water to form modified CDW (mCDW). Concurrent shipboard CTD and ADCP data are used to diagnose the responsible mixing processes and highlight the importance of thermohaline intrusions. The genesis mechanism of the interleaving layers cannot be confirmed, however a simple analytic model suggests the upper limit contribution of advection by internal waves cannot account for the observed temperature variance unless the cross-eddy temperature gradient is large. Data-adaptive sampling of an eddy with the glider revealed it lost heat across two isopycnals and a fixed radius at a rate of 7 × 109 J s−1 over 3.9 days. This rate is corroborated by a diffusion model initialized with the eddy’s initial hydrographic properties and informed by the heat fluxes parameterized from the shipboard data. The results suggest eddies predominately lose heat laterally and downward, which preserves subsurface heat for melting of marine-terminating glaciers.
The West Antarctic Peninsula (WAP) is bordered by the Antarctic Circumpolar Current (ACC), which flows along the continental slope. Below the surface layer, the ACC advects a warm mass of Circumpolar Deep Water (CDW) with significant heat content relative to the in situ freezing point. The CDW spans a range of properties, though Gordon (1971) distinguishes between an upper (UCDW) and lower (LCDW) variety defined by temperature and salinity maxima, respectively. The southern boundary of the ACC is defined as the southernmost presence of UCDW (Orsi et al. 1995) and, unique compared to the rest of Antarctica, near the WAP the ACC flows immediately adjacent to the shelf break, making UCDW readily available to the shelf (Fig. 1).
The WAP is undergoing rapid climate change and the ocean is a primary heat source, particularly in winter when there is no direct radiative forcing. The marginal seas of West Antarctica are warming (Schmidtko et al. 2014) and the increase in heat content along the WAP margin has been attributed to a warming of the UCDW Tmax (Martinson et al. 2008). Cook et al. (2016) confirm an oceanic role in glacier retreat along the WAP by demonstrating an asymmetry in glacial advance/retreat: southern marine-terminating glaciers that have access to warm subpycnocline waters are retreating whereas northern glaciers under the influence of much colder Bransfield Strait waters are not. More and/or warmer CDW has also left its imprint on the atmosphere. The northern portion of the WAP is undergoing rapid winter warming (Turner et al. 2013), presumably related to lighter sea ice cover venting ocean heat to the atmosphere. The reduced sea ice cover, in turn, may be related to changes in the winds and their effect on UCDW delivery and mixing across the pycnocline (Dinniman et al. 2012).
The myriad consequences of UCDW heat have made the exchange of CDW with the WAP shelf an active area of research (Klinck 1998; Smith et al. 1999; Klinck et al. 2004; Dinniman and Klinck 2004; Moffat et al. 2009; Dinniman et al. 2011, 2012; Martinson and McKee 2012; Spence et al. 2014, 2017; Graham et al. 2016; Couto et al. 2017), which is summarized in a recent review by Moffat and Meredith (2018). Various processes have been implicated in driving the cross-isobath transport of CDW onto the WAP shelf, each of which may dominate on different time and/or space scales. The importance of the mesoscale has been argued for by theoretical means (Stewart and Thompson 2015), demonstrated in high-resolution numerical models (St-Laurent et al. 2013; Graham et al. 2016; Stewart et al. 2018), and observed (Moffat et al. 2009; Martinson and McKee 2012; Couto et al. 2017). Warm-core, subpycnocline, primarily anticyclonic eddy-like features have been found within several tens of kilometers from the shelf break, particularly in the vicinity of Marguerite Trough, and are steered along isobaths. The hydrographic properties of the eddies indicate an injection of UCDW as it is found on the continental slope and their length scale is as large as or slightly larger than the first baroclinic Rossby radius, which near the shelf break is about 5 km. Decorrelation lengths of physical and geochemical scalars on the WAP shelf are also about 5 km (Eveleth et al. 2017), suggesting mesoscale eddies may dominate tracer stirring. The eddies’ large heat content relative to surrounding waters indicates a potentially large onshore heat flux. For example, in numerical models, cumulative onshore heat transport is reduced by as much as 50% when model grid spacing is increased from 1 to 2 km (St-Laurent et al. 2013).
Still, very little is known about the origin of the eddies. Their core hydrographic properties and stratification suggest an origin along the continental slope and their observation along isobaths suggests they are, at least in part, advected within the mean flow. Consistent with the former, an idealized numerical model with a continental slope straddled by a jet similar to that observed along the WAP suggests that the jet soon becomes unstable and a Rossby wave containing alternating warm anticyclones and cold cyclones emerges (St-Laurent et al. 2013). Additionally, it is important to understand the processes that attenuate the mesoscale variability and work toward setting the larger-scale shelf stratification. Once UCDW intrudes onto the shelf it mixes to become modified CDW (mCDW), but attempts to understand UCDW transformation have been based on shelf-integrated budgets and have been process independent (Klinck 1998; Smith et al. 1999; Klinck et al. 2004).
In this study we seek to both understand the origin of mesoscale eddies observed on the WAP shelf and then identify and quantify the major processes responsible for their heat loss. Our focus is on the southern grid region of the Palmer Long Term Ecological Research project (Pal LTER; Smith et al. 1995) and in particular the vicinity of Marguerite Trough as this is a region of frequent eddy activity, is close in proximity to the shelf-break jet that we hypothesize generates the eddies, and contains well-defined bathymetric pathways by which to steer the eddies. We primarily use data from a novel Slocum glider survey designed to sample a known pathway for CDW exchange, documenting gradients along the axis of advection for any eddies encountered along the way, and then identify and track an eddy in real time through data-adaptive sampling. This allows, for the first time, high-spatial and temporal resolution transects directly through some of these eddies and a real-time quantification of the attenuation of their core properties. We supplement this with shipboard CTD and ADCP data used to quantify the mixing processes in the WAP environment around and within the eddies and to diagnose the instability that generates the eddies. Ultimately, we simulate the evolution of the eddy as documented by the glider with a simple diffusion model informed by our parameterized mixing processes.
2. Data and observations
The principal dataset used in this study is a set of temperature and salinity profiles collected by a Teledyne–Webb Slocum glider (Schofield et al. 2007) equipped with an unpumped Sea-Bird CTD. Slocum gliders traverse the water column (to 1000 m) in a sawtooth pattern by changing their buoyancy, traveling with average horizontal speed of 0.35 m s−1. There is a hysteresis effect apparent in the up versus down traces that we correct by applying the thermal lag correction of Garau et al. (2011). We also empirically correct for a salinity bias by regressing glider-recorded salinity against ship-recorded salinity, each averaged in subpycnocline temperature bins, for five stations that were occupied by both platforms (maximum temporal separation of 12 days). To facilitate analysis, all corrected up and down traces are binned into 1-m profiles whose latitude, longitude, and time coordinates are taken as the average of all intraprofile samples.
The glider was deployed at Palmer Station at the head of Palmer Deep Canyon (PD) and traveled downgrid before reaching grid station 400.100 (see Fig. 2 for grid convention and physical setting), at which point it began the first of two phases of sampling. In the first phase, the glider flew against the mean current within the anticyclonic cell extending out of Marguerite Trough in order to sample any eddies being advected across the shelf and to identify lateral gradients in their properties. In the second phase, after reaching station 290.115 within Marguerite Trough, it waited to come across an eddy in order to track it in real time. The sampling strategy and basic observations of each of these stages are described below.
a. Glider survey: Advective Path
1) Sampling strategy
The Advective Path (highlighted in Fig. 2) refers to the anticyclonic cell extending out of Marguerite Trough that carries upwelled UCDW to the northern portion of the grid. This path was confirmed to carry UCDW by Martinson and McKee (2012), and its width and central location were inferred from the depth-averaged currents and CDW dye transport in Dinniman et al. 2011 (their currents superimposed on our Fig. 2). Because the width of the current is about the same as the expected diameter of the eddies we were confident that by flying straight through it we would be able to cross any existing eddies traversing the shelf. Nominal glider profile spacing is 1 km, which should afford several profiles through each eddy.
2) Eddy characteristics
As the glider flew upstream (downgrid), it encountered five eddies imbedded within the subpycnocline Tmax layer. We define a mCDW temperature profile as the along-isopycnal average of all profiles whose Tmax is <1.55°C and then compute heat content per unit area Q of the eddy profiles relative to this mCDW profile by differencing the two and integrating between σshal = 1027.64 kg m−3 and σdeep = 1027.76 kg m−3. We use neighboring mCDW as our reference profile because we are interested in how the eddies mix with their surroundings. We choose these isopycnals since 1) they encompass a positive temperature anomaly within eddies, 2) σshal is relatively stable within the permanent pycnocline, and 3) σdeep is near, but does not intersect, the bottom. Series of temperature, salinity, heat content per unit area, and geostrophic current at 280 dbar relative to the bottom are shown in Fig. 3, with the eddies labeled A–E (geographic locations of eddies A–E shown in Fig. 2). Fundamental statistics for each eddy developed and discussed below are given in Table 1.
To quantify the dimensions, gradients, and heat content of these eddies, we obtain an analytical representation for the eddy in terms of a Gaussian function that is fit to the data via nonlinear regression. Specifically, we identify intraeddy profiles via positive heat content and assign them an along-eddy coordinate χ along an axis determined via orthogonal regression of profile grid line against grid station. This allows inversion of the following model
for core heat content per unit area Qmax, center χ0, and radius R. Implicit in this model is the assumption of axisymmetry, which we maintain throughout the paper. Our model is similar to those used by other authors (Couto et al. 2017), however we fit the function to integrated heat content as opposed to along-isopycnal temperature anomalies since integration smooths over some of the locally high temperature variance that can reduce the quality of fit. The total eddy heat content is given by the radial and azimuthal integrals of Q:
It is important to note that radii (and heat content) determined from this method are underestimates because the glider may not be going straight through the eddy (they are half chord lengths) but also because the glider is deliberately flying against the mean flow which is advecting the eddy through the glider during sampling. Therefore, we provide an additional measure of eddy size as half the along-track distance spanned by intraeddy profiles (RTD), which tends to be slightly larger.
The cores of the eddies contain as much as 5.5 × 108 J m−2 relative to the reference profile. In T–S space, the properties of these eddies are consistent with UCDW. Particularly so, eddy C contains water disjoint from neighboring water, consistent with UCDW as found upstream on the continental slope and indicating warm-core isolation. The upstream eddies A–C are narrower than their downstream counterparts D–E, having radii R ~3 km (RTD ~5 km) compared to R ~7 km (RTD ~9 km). Eddies D–E are broader than the upstream eddies but do not have a substantially smaller heat content per unit area and in fact tend to have more total heat. The large discrepancy in upstream versus downstream total heat content might be more reflective of a sampling bias (offset trajectory through eddy center and stronger head current encountered through eddies A–C) as opposed to an actual difference in heat content. The injection of UCDW into the shelf water column means that the eddies are associated with a downward deflection of isopycnals and anticyclonic shear (Fig. 3). The deflection of isopycnals is largest in the weakly stratified portion of the water column at and below the Tmax, which leads to a cross-track geostrophic velocity signal. To quantify this, we evaluate the geostrophic current at 280 dbar relative to the bottom. We choose this depth as it is the approximate depth of the moored current meters used by Martinson and McKee (2012) upstream where they found the largest eddy signal. A composite over all of the eddies (Fig. 4) reveals a well-defined signature of anticyclonic rotation centered about a warm core, similar to the composites of Moffat et al. (2009).
We measure temperature gradients at the eddy boundaries via a least squares approach. We define by regressing a straight line to the portion of the temperature profile spanning 10 m above and below the depth of σshal and then averaging the slopes laterally across all profiles within . We define similarly, considering the region 40 m above and below the depth of σdeep. To estimate , we fit a straight line on each isopycnal to the temperature value at the profile nearest to and up to two values on either side of that profile and then average the slopes between σshal and σdeep and across both hemispheres. Owing to the nature of glider sampling, for the same reasons radius estimates are underestimates, lateral gradient estimates are overestimates. Vertical gradients are similar for all eddies while lateral gradients are larger for the three upstream eddies.
b. Glider survey: Tracking Stage
1) Sampling strategy
We conducted transects along a “fence” perpendicular to the eastern wall of Marguerite Trough (Fig. 2, blue line) to find an eddy and then track it via real-time adaptive sampling. Given the local Rossby radius (~5 km), the mean flow speed, an assumption that eddies are advected within the mean flow, and the lateral deviation of the mean flow, we determined that the fence should be 15 km wide and that it should be surveyed back to back in a 24 h period, which is within the operational constraints of the glider. During sampling, profiles were inspected for a Tmax ≥ 1.8°C as evidence of a UCDW-core eddy and, if found, trajectories were forecast by integrating along both the vector of 24-h mean glider depth-averaged currents and along a streamline fit to the Dinniman et al. (2011) model time–depth-averaged currents, both of which generally agreed within one Rossby radius. We encountered an eddy shortly after initiating sampling and crossed it five times over 4 days. The trajectory followed the eastern wall of Marguerite Trough, consistent with the eddy-like features observed by Moffat et al. (2009).
2) Eddy characteristics
We define a local reference profile in the same manner as before and subtract it in order to obtain Q. Analogous to the Advective Path, sections of data are shown in Fig. 5. However, unlike in that stage, we are not confident that we consistently crossed the eddy via a chord length. Connecting spatial end-members of threshold Tmax along a line perpendicular to isobaths leads to an estimated mean diameter of 8.5 km, similar to the eddies observed on the Advective Path. There is, in general, a decrease in heat content per unit area over time. Anticyclonic shear is less apparent beyond the first few crossings, and the vertical structure of temperature and salt anomalies is much more complicated. We will quantify the heat loss of this eddy in section 5.
c. Shipboard data
1) Pal LTER CTD and SADCP
Shipboard data are used to supplement the glider data and to apply mixing parameterizations. CTD profiles are collected as part of the standard sampling on the annual cruises to the WAP each austral Summer since 1993 and are currently (since 1999) collected with a dual-pumped Sea-Bird 911+ CTD system (see Martinson et al. 2008 for details). The standard grid locations are indicated as black squares in Fig. 2. The entire grid was occupied until 2008 whereas now only a subset of stations is occupied (nominally one coastal, shelf, and slope station per grid line), though this is generally complemented by various process-study CTD casts at nongrid locations that are not plotted but are utilized here.
Processed, high-resolution (5 min in time, 8 m in vertical) velocity profiles from the ARSV L. M. Gould’s hull-mounted RDI 150 kHz narrowband instrument were obtained from the Joint Archive for Shipboard ADCP (SADCP). The 150 kHz instrument provides velocity profiles good to about 300 m (depending on weather and sea state) when the vessel is on station. For mixing parameterizations we need concurrent shear and stratification profiles, so only the overlapping portion of the database is used (Januaries 2000–15; SADCP installed in mid-1999). For those applications, the hydrographic data are bin-averaged onto the same 8-m grid of the velocity data. Stratification and shear are computed as first differences on the 8-m grid as and . For bins with very weak stratification (N2 < 1 × 10−6 s−2), N2 is set to a constant value (1 × 10−6 s−2).
2) Casts on Advective Path
We draw special attention to 5 CTD casts collected in January 2014 along the Advective Path that was sampled by the glider the year prior (triangles in Fig. 2). These casts fortuitously sampled one or more eddies and indicate substantial interleaving structure, particularly in the pycnocline (Fig. 6). The layering connects cores of injected slope-type UCDW to the surrounding cooler pycnocline. These data are used to assess the importance and origin of thermohaline intrusions into the eddies.
3) S04P casts across continental slope
We suspect the eddies are generated along the continental slope. To diagnose the structure and stability of the shelf-break current upstream of Marguerite Trough we use CTD data from the WOCE S04P cruise in February 1992 (“x” symbols in Fig. 2). Though it is possible that the hydrographic structure over the continental slope has changed in the long interim period, we use these data as they provide a high spatial resolution transect across the slope, slightly finer than the 2011 reoccupation of the line. The stations are “moved” from their original location to the 200 line by tracing the isobath at the actual station location to the 200 line, a small correction (a few kilometers). Neutral density is computed using the Jackett and McDougall (1997) software and profiles of neutral density, temperature, and salinity are filtered with a third-order Butterworth filter with width 100 dbar to remove noise and filter out Charney-type instabilities. We choose to do this because surface-trapped Charney-type instabilities do not convert significant available potential energy (APE) to eddy kinetic energy (EKE) compared to the pycnocline-level instabilities we are seeking (Smith 2007). Geopotential anomaly is calculated from the smoothed temperature and salinity profiles and is linearly extrapolated to handle bottom triangles. Both neutral density and geopotential anomaly are mapped with a Gaussian weighting function to a two-dimensional grid across the slope using the WOCE global climatology vertical grid (Gouretski and Koltermann 2004) and uniform 5-km horizontal spacing.
3. Origin of mesoscale structure
a. QG model and background state
To understand the theoretical characteristics of eddies on the WAP, we begin by diagnosing the stability of the shelf-break current upstream of where UCDW intrudes onto the shelf. We follow the inviscid quasigeostrophic model of Smith (2007), assuming a local, slowly varying mean stratification N and horizontal flow U depending only on z. Assuming plane-wave solutions where is the complex amplitude of the perturbation streamfunction, ω the complex perturbation frequency, and k = (k, l) the wavevector, the linearized quasigeostrophic equations in geographic coordinates lead to the eigenvalue problem
Here, is the potential vorticity stretching operator and αx,y are the bottom slopes. We caution that the assumption of a linearly growing wave in a stable background state that is steady in time is highly idealized. Surely interactions with the evolving turbulence field of the ACC are relevant for the dynamics of the shelf-break flow. Nevertheless, we suspect that the regularity of the eddies and their common statistical metrics (Figs. 3, 4; Moffat et al. 2009; Couto et al. 2017) demonstrate that there is a preferential length scale and vertical structure for instabilities in this environment.
The background state is defined by the gridded stratification and geostrophic velocity profiles at grid location 200.150. This site is chosen as it is located midway across the slope and the geostrophic shear there agrees well with climatological SADCP shear. The current is assumed to flow parallel to isobaths, specifically at an angle of 40° north of east. The local bottom slope is calculated by fitting a 2D plane to all ETOPO1 bathymetry data (Amante and Eakins 2009) within a radius of 0.25° about the grid location.
One assumption of this model is that the background state varies slowly in the horizontal. Since there is clearly lateral shear in the shelf-break current, we need to justify excluding the potential of barotropic instability. For a generic, mixed baroclinic–barotropic instability, the contribution to the growth rate from baroclinic instability scales as and the contribution from barotropic instability scales as (Pedlosky 1987). For these data, the latter is one order of magnitude smaller than the former (σBC ~ 13 day−1, σBT ~ 2 day−1). This is corroborated by Stern et al. (2015) who found baroclinic instability to grow much faster than barotropic instability in their QG model of a similar shelf-slope configuration, though their model was fully turbulent and not limited to the linear growth stage.
b. Most unstable mode
The equation is discretized as in Smith (2007) and solved for a range of wavenumbers. For each wavenumber, the most unstable mode is that with the largest imaginary part of ω. We find the overall most unstable mode to have an inverse wavenumber |kmax|−1 = 4.4 km and growth rate 0.4 day−1 (Fig. 7). Over the wavenumber space evaluated, this is both the global and only local maximum. The mode’s vertical structure is characterized by nonzero amplitude below the permanent pycnocline in the CDW depth range with a maximum at about 600 m and a secondary maximum at about 350 m. Specifically, this depth range spans the water column presence of UCDW (T ≥ 1.7°C) and LCDW (S = Smax) within the four profiles across the slope. The instability appears to be qualitatively similar to a Phillips type instability (Phillips 1954) for the following reasons: 1) vertical structure seems to be dictated by interior sign changes in the potential vorticity gradient, 2) ~33° phase shift within the amplitude maximum indicates APE release there, and 3) horizontal scale obeys L ≈ (N/f)hpycnocline ≈ 5 km.
Overall these findings are in very good agreement with the glider observations over the Advective Path (Fig. 3). First, the vertical structure of the mode is concentrated within and below the permanent pycnocline, which is where the eddies exist. Second, the observed diameters of eddies A–C upstream and average crest-to-crest separation are ~8.3 km (averages of 2R and 2RTD) and 22.3 km, respectively, which are comparable to the theoretical diameters 2|kmax|−1 = 8.8 km and crest-to-crest separation 2π|kmax|−1 = 27.8 km for eddies spun off of an unstable Rossby wave. Temperature anomalies relative to mCDW and geostrophic currents provide some evidence for a cold cyclone between eddies B–C, though in general, only warm anticyclones tend to be found on the shelf. The eddies D–E downstream on the path are broader than A–C but have a similar separation. With the caveat that this is a linear model, note that the wavelength implies delivery of (365 days) × u|kmax|/2π = 57–114 eddies per year (for typical shelf currents of 0.05 or 0.10 m s−1), which could account for the totals observed by Moffat et al. (2009) and Martinson and McKee (2012).
c. Roles of bottom slope and current orientation
Exploring the parameter space of bottom slopes and current orientations allows us to simultaneously understand the sensitivity of the most unstable mode’s structure and growth rate to these parameters and to understand how generalizable these results are to other regions around Antarctica. Using the same geostrophic shear and stratification profiles, we repeat the above analysis for all bottom slopes between −0.15 and +0.15 at increments of 0.01 and for all current orientations between 0° (zonal) and 90° (meridional) counterclockwise from east at 10° increments and examine changes in the growth rate, wave vector, and structure of the most unstable modes. Evaluating positive and negative bottom slopes effectively allows consideration of both prograde (isopycnals slope in same sense as bathymetry; e.g., WAP) and retrograde (isopycnals slope in opposite sense as bathymetry; e.g., Ross Sea) jets.
Figure 8 shows the growth rate and inverse wavenumber of the most unstable mode over the entire parameter space. The orientation of the current has essentially no effect on the instability. This is likely because the planetary beta effect is so small at this latitude that potential vorticity gradients are dominated by the stretching term and bottom slope. Indeed the topography plays a large role in determining the strength and properties of the most unstable mode. It is found that negative bottom slopes are stabilizing while increasing positive bottom slopes are destabilizing to a point and then stabilizing. Blumsack and Gierasch (1972) demonstrate that the relevant parameter for the stability problem under QG scaling is not the bottom slope itself but rather the ratio of the bottom slope to the isopycnal slope, δ ≡ α/s (Blumsack and Gierasch 1972; Poulin et al. 2014). In the QG model of Smith (2007) that we use, the bottom slope only enters the problem as a boundary condition in the bottom layer’s potential vorticity gradient. When the bottom slope exceeds the slope of the 1028.0 kg m−3 neutral surface, the potential vorticity gradient vanishes in the lower layer which has the effect of suppressing the growth rate in accordance to the Charney–Stern criteria (Pedlosky 1987). Nevertheless, this does not change the fact that there are still sign changes in the interior potential vorticity gradient and therefore instability persists where δ > 1 or δ < 0 (Isachsen 2011). These regions of weaker instability with inverse wavenumbers of 4–5 km all have similar vertical structure as the most unstable mode obtained in the realistic scenario analyzed earlier.
The fastest growing modes for prograde currents with relatively flat bottoms (0 ≤ δ ≤ 1) are qualitatively consistent with Eady modes for the following reasons: 1) they have amplitude maxima at both ~350 m and near the bottom; 2) they have an inverse wavenumber L ≈ NH/(1.6f) ≈ 9.4 km, probably since the weakly varying N is dynamically similar enough to the uniform N of Eady’s model. These modes have a pronounced spike at ~350 m near the UCDW temperature maximum and, in regions with a flatter bottom slope than the WAP, could also be expected to contribute to exchange of CDW. The region between −2 ≤ δ < 0 represents a sort of transition between the Eady-type modes and the Phillips-type modes. They have smaller length scales (~2 km) and are bottom-boundary trapped. This trend fits the qualities described by Blumsack and Gierasch (1972) for increasingly negative slope parameter, namely decreasing length scale and boundary trapping.
4. Attenuation of mesoscale structure
Having identified a plausible origin for the UCDW eddies, we here consider the processes responsible for their decay. Box inversions of steady heat budgets on the WAP point to diapycnal mixing between overlying remnant winter mixed layer Winter Water (WW) and a constantly replenished CDW layer as the maintenance of the permanent pycnocline. For some context, using a steady advective–diffusive balance, Klinck et al. (2004) place an upper bound on the vertical (lateral) diffusivity of heat at 7.7 × 10−4 m2 s−1 (1600 m2 s−1) in the limit of no lateral (vertical) mixing. In an earlier study, Klinck (1998) used seasonal changes in water properties and a similar integrated budget to find a vertical (lateral) diffusivity of heat as 1 × 10−4 m2 s−1 (37 m2 s−1). Martinson et al. (2008) used interannual variability of WW heat content and assumed a UCDW replenishing time to suggest a vertical diffusivity of 8.5 × 10−5 m2 s−1. Our study, however, focuses on how the properties of the subpycnocline “box” are set by the mixing of advected parcels of UCDW, within which the vertical and lateral gradients are quite different from those used in mean-shelf balances.
a. Shear-driven instability
Shear-driven instability is thought to be important in maintaining the permanent pycnocline and is thought to yield larger heat fluxes on the WAP than double-diffusive instabilities (Howard et al. 2004). Regarding sources of shear, internal tides are likely not important but near-inertial waves may be. Howard et al. (2004) suggest the semidiurnal tide is weak, as is the stratification, making baroclinic conversion unlikely. Further, Beardsley et al. (2004) use rotary spectral analysis of drifter velocity to show power in the semidiurnal band is two orders of magnitude less than that in the near-inertial. In general, we might expect that vertical mixing should be elevated over seamounts or within cross-cutting canyons. At seamounts, internal waves with frequency at the critical slope may be generated whereas in canyons internal waves may become trapped and focused toward the canyon head (Gordon and Marshall 1976).
Because we do not have microstructure measurements within eddies or on the surrounding midshelf, our approach is to use SADCP and CTD data to estimate diffusivities, temperature gradients, and heat fluxes across the permanent pycnocline, which we define globally to be between 98 and 250 m. We exclude grid stations below the 000 grid line as the hydrography there is very different (very deep and cold remnant winter mixed layers). Note that winter mixing by entrainment of the thermocline during brine rejection is not considered here, but is important in the annual heat budget (Martinson and Iannuzzi 1998).
1) General shear-driven instability
The method of Pacanowski and Philander (1981, hereafter PP81) computes a diffusivity Kz(z) as a function of the Richardson number, with the idea being that when shear overcomes stratification, instability and mixing result. This method was developed for steady currents in equatorial ocean models and assumes nothing about the underlying sources of shear. Nevertheless, it has been applied to the WAP (Howard et al. 2004) and other high latitude environments (Dewey et al. 1999). For every 8-m binned CTD cast we have we pair it with all concurrent on-station SADCP profiles and compute a time-averaged Richardson number . Because this is essentially a space–time-averaged Ri based on finite-differenced data, it is best interpreted as a probabilistic measure of potential instability at space (and time) scales lower than the differencing (and averaging) scales. The PP81 parameterization is
and a vertical heat flux can then be calculated at each depth as
2) Internal wave parameterizations
Various authors have modeled the energy transfer through the steady Garrett–Munk (GM; Garrett and Munk 1975) internal wave vertical wavenumber spectrum via nonlinear interactions down to dissipation scales and turbulence production by scaling the GM shear spectrum by observed shear variance (e.g., Gregg 1989, hereafter G89). By assuming a steady turbulent kinetic energy balance and assuming a constant mixing efficiency (Γ = 0.2), one can then estimate a diffusivity coefficient from parameterized dissipation via the relation . We use a modified version of the G89 parameterization, which is
In this expression, is the 8-m SADCP shear squared, multiplied by 2 (correcting for the finite-difference filter; Gregg and Sanford 1988), and then squared again and averaged in time; is the GM shear spectrum (with local scaling parameters for the WAP shelf; see appendix) integrated up to a cutoff wavenumber β = 0.2π rad m−1, squared, and then multiplied by 2 (assuming Gaussian statistics); and K0 = 5.6 × 10−6 m2 s−1 is a nominal diffusivity for the underlying internal wave spectrum (see appendix). The functions h and j are corrections added later and are not part of the original G89 model.
First, h(Rω) is a polynomial function of Rω, the N2-normalized-shear–strain variance ratio, which is designed to adapt the model to non-GM wave fields with different aspect ratios and frequency content:
The GM spectrum has Rω = 3 (h = 1), and larger values of Rω indicate elevated importance of near-inertial waves. Most of the open ocean has Rω ~7 (Kunze et al. 2006), and values in the Southern Ocean are generally between 8 and 12 (Thompson et al. 2007; Naveira Garabato et al. 2004). Calculating shear and strain spectra for each cast as in Kunze et al. (2006), we find Rω = 9 (h = 0.42) averaged over the LTER grid, and, owing to the difficulty involved in estimating strain, we use that constant value for all casts.
The second term j(f, N) is a correction for latitude arising because the rate at which waves are Doppler shifted depends on the ratio of their horizontal and vertical wavenumbers, which depends on Coriolis frequency f (Gregg et al. 2003):
3) Results and interpretation
Fundamental results are presented in the form of profiles of N2, S2, and Kz averaged over the Shelf region (see Fig. 2 for regional boundaries) in Fig. 9. In general, S2 decreases more slowly with depth than does N2, becoming nearly white below the pycnocline and yielding diffusivities that increase with depth. To confirm that the apparently unusual behavior of S2 with depth is real, we compare the SADCP shear at station 200.140 to the shear calculated from a set of six moored current meter records at that site [corrected for finite differencing following Gregg and Sanford (1988) using local GM scaling]. The current meter shear (bold blue line) and SADCP shear (thin black line) at site 200.140 agree remarkably well, suggesting that, at least at this site, the large shear variance (and hence the large diffusivities) at depth is real (Fig. 9b).
While the SADCP shear variance represents an integral across all frequencies, we can bandpass filter the current meter observations between [f, N0] to retain only shear in the internal wave band. Doing so (dashed blue line in Fig. 9b) suggests that the internal wave band shear profile has the same shape as the total shear profile but that SADCP shear might be overestimating internal wave shear by a factor of ~1.6. Because we cannot evaluate this relation at other sites, we do not attempt to make any correction for non-internal wave shear. The histogram of in the pycnocline depth range (Fig. 9e) suggests observed shear variances are generally a factor of 2–3 greater than GM. Spectra of shear and strain in the pycnocline confirm that Rω = 9 is a good fit on the WAP shelf (Fig. 9d) and reveal an elevated importance of near-inertial waves in comparison to a standard GM spectrum.
Both the G89 and PP81 diffusivities increase with depth for reasons discussed above. Diffusivity values in the permanent pycnocline are small (<10−5 m2 s−1), so much so that the asymptotic lower limit of the PP81 method renders it inapplicable there (Fig. 9c). Although the two methods agree better deeper in the water column, all further analyses will exclusively consider the G89 parameterization. The shear (Fig. 9b) and diffusivity (Fig. 9c) profiles are presented as means with the shaded region enclosing two standard errors, where both the means and standard errors are computed on log-transformed data. Owing to the large number of stations sampled, uncertainty in the mean is small.
To get a sense of how diffusivity and vertical heat fluxes vary regionally, Fig. 10 shows composite profiles of diffusivity and heat flux for data partitioned by region alongside histograms of the pycnocline-averaged quantities. The small diffusivities that characterize the permanent pycnocline operating on the background T profile yield vertical heat fluxes across the permanent pycnocline of <1 W m−2 which is smaller than values reported by Howard et al. (2004), who used PP81. Diffusivities increase approximately log-linearly with depth to about 10−4 m2 s−1 at 300 m (Fig. 10), which is as deep as the SADCP yields useable data.
As in Fig. 9, the central values are computed on log-transformed data, however now the spread is indicated by two standard deviations in order to indicate potential values instead of uncertainty in the mean. While the central values are small, the spread at any given depth is large and spans about two orders of magnitude, suggesting that the potential for strong mixing at any depth is high and that mixing is intermittent. Pycnocline diffusivities are significantly different between the Shelf and Slope, the Shelf and Coast, and between the Northern Shelf and Southern Shelf at an α = 10% level. The variance seems to be larger in the Coast. The larger variance there is driven by larger variance in the shear. We also constructed two composites that should be more representative of the environment that most of the eddies are in. The first averages all profiles in the subregion of the Southern Shelf that is bounded by the 250 and 450 lines between stations 030 and 130 (the vicinity of Marguerite Trough) and the second averages only those profiles in that domain with a Tmax ≥ 1.55°C. The heat fluxes in those regions are qualitatively different with a much greater heat flux divergence centered about the Tmax but are statistically indistinguishable from the more comprehensive Southern Shelf in terms of cross-pycnocline metrics. However, compared to the Northern Shelf, the Southern Shelf and both Marguerite Trough composites have significantly lower pycnocline diffusivities, significantly lower shear variance ratios, and insignificantly larger pycnocline heat fluxes. The results suggest that any enhanced vertical heat flux within an eddy compared to its surroundings is due to the altered gradients of the temperature profiles and not the shear.
While the SADCP does not allow for estimates of diffusivity or heat flux below about 300 m, extrapolating the log-linear trend in Kz with depth and using observed temperature gradients suggests that below the Tmax, the vertical temperature gradient is about one order of magnitude smaller than that in the permanent pycnocline, however the diffusivity should be about two orders of magnitude larger. This combines for a heat flux about 10 times greater below the eddy, and the spatial variability in the heat flux profiles is consistent with this scaling. The stratification on the Slope and Southern Shelf is dominated by pure UCDW and a middepth temperature maximum as in the eddies. The middepth maximum in temperature results in an upward heat flux above the Tmax and a downward heat flux below, the former of which is apparent in the observed profiles and the latter of which is implied by the trend toward a ~0 W m−2 heat flux at 300 m (Fig. 10). Because the downward heat flux is about 10 times greater than the upward heat flux, the temperature maximum is mixed downward in the water column as the UCDW is advected northward and shoreward. This is consistent with the positive vertical heat flux increasing with depth on the Northern Shelf and Coast stations (Fig. 10).
b. Thermohaline intrusions
Four of the five CTD casts taken along the Advective Path reveal that within and surrounding the eddies there is substantial along-isopycnal temperature and salinity variance associated with thermohaline intrusions (Fig. 6). Those data are plotted as profiles of temperature in Fig. 11 along with a background profile, which is constructed as a running median in density space with a window size of 0.015 kg m−3. The UCDW eddies are essentially moving fronts. Joyce (1977) derived a model in which medium-scale advection of heat and salt across a density-compensated large-scale front is balanced by small-scale diffusion across thermohaline intrusions, attenuating their T–S characteristics. The medium-scale advection is taken to be in the form of alternating interleaving structures, initiated by velocity perturbations with the energy source coming from thermoclinic energy of the cross-frontal property contrasts. His model makes no distinction as to what small-scale processes conduct the mixing.
The model balance leads to an effective cross-frontal diffusivity for T of
where primes indicate small-scale variance, overbars indicate large-scale averages (scales larger than intrusions), and the vertical eddy diffusivity is taken to encompass all mixing processes. To apply Eq. (4.6) we estimate the cross-frontal temperature gradient from the glider data along the Advective Path by averaging the lateral gradients of the 5 eddies encountered (A–E; 8.5 × 10−5 °C m−1). The small-scale intrusion variance in the vertical is derived from the CTD profiles at approximately the same location: the background profile is subtracted and the variance of the derivative of the residual is calculated. We calculate Kh in two ways: once using a constant Kz typical of the midwater column (10−4 m2 s−1; Hebert et al. 1990, Pelland et al. 2013) and once using the diffusivities as parameterized from the modified G89 method. The results are 3.2 ± 1.4 and 7.9 ± 4.2 m2 s−1, respectively (two standard errors indicated). We now consider processes that may be responsible for generating the thermohaline intrusions.
1) Double-diffusive growth
Theories for the growth of thermohaline intrusions (e.g., McDougall 1985a,b) in a region of density-compensated thermohaline gradients depend on one of the components of density being unstably stratified in order for an infinitesimal disturbance to grow to finite length. Figure 11 shows that the WAP water column between the WW Tmin and the UCDW Tmax is diffusively unstable and that the layer immediately below the UCDW core can be salt finger unstable. One way to quantify how much the unstably stratified component contributes to the stability of the water column is via the density ratio . A statically stable column is diffusively unstable if 1 < Rρ < ∞, salt finger unstable if 0 < Rρ < 1, and doubly stable if Rρ < 0. As the ratio approaches 1, double-diffusive fluxes increase because the effect of the gravitationally unstable gradient is increasingly canceled by that of the stably stratified component’s gradient. Bormans (1992a) found that diffusively unstable fronts had lateral intrusion heat fluxes significantly exceeding those of doubly-stable fronts only when Rρ went below 1.54. Values in the lower pycnocline/eddy upper hemisphere occasionally breach that threshold, and there is a correspondence between the locations of largest temperature variance and the locations of density ratio nearest to 1 (Fig. 11, lower panels).
2) Internal wave advection
Another possibility is that the intrusion-scale temperature variance is caused by advection due to internal waves. If we consider the simplest advection equation and the square of its Fourier transform,
we can rearrange it to construct an inequality for the maximum contribution by internal waves to the intrusion-scale (tilde) temperature variance:
Here we have assumed that all waves are purely horizontal (inertial waves) and we have differentiated in depth by multiplying each side by β2, the vertical wavenumber squared. Shear spectra are calculated over our defined pycnocline range and temperature spectra are calculated over [180, 400] m, encompassing the eddy itself, the range of intrusions, and avoiding remnant diffusive staircases. The large-scale (overbar) temperature gradient is again the average across eddies A–E (Table 1). Figure 12 shows the two sides of Eq. (4.8) and it implies that observed temperature variance is greater than the maximum contribution of internal waves acting on the cross-eddy temperature gradient. It is worth noting, however, that if we instead use the lateral temperature gradient of only the three upstream eddies (A–C; 1.2 × 10−4 °C m−1) then the 90% error bars in Fig. 12 begin to overlap in the range of intrusion thicknesses.
3) Geostrophic turbulence
Under geostrophic turbulence, the steep roll-off of energy spectra in comparison to the more gradual roll-off of potential enstrophy and tracer variance spectra means that the density field is dominated by the large, low-mode energy containing eddies that stir along isopycnals (such as the first-mode baroclinic instability; section 3) whereas the T–S variance is dominated by small-scale filamentary features that must be density compensated (Smith and Ferrari 2009). The cascade of tracer variance is halted at high vertical wavenumbers by vertical mixing, which counters the variance production by mesoscale stirring. If the along-isopycnal passive tracer variance surrounding the eddies was generated by their own stirring, we would expect their sizes to match the eddy length scale associated with tracer variance generated by stirring across a mean gradient in accordance with mixing length theory (Tennekes and Lumley 1972):
where C is some passive tracer, x is the cross-slope coordinate, and overlines indicate an average in time or space.
To unite density-compensated variations into a single variable we consider spice (Flament 2002) which effectively serves as a quantification of distance in T–S space orthogonal to isopycnals. We compute spice profiles C(z) from the glider data along the Advective Path in addition to high-passed spice profiles CHP = C − CLP, which are the total profiles minus a filtered profile (triangular filter with width 25 m), thus retaining only the interleaving layers. Inspection of composites of across eddies A–E, where the over line indicates averaging in the vertical between [σshal, σdeep], reveals that thermohaline variance is indeed largest at the eddy edges (Fig. 13). To visualize how spice variance is distributed more generally along the Advective Path, Fig. 14a shows a section of C on isopycnals and Fig. 14b collapses along-isopycnal variations into profiles of and , where now primes are relative to the along-isopycnal mean and the products of anomalies are averaged on each isopycnal. Unsurprisingly, exhibits largest values in between [σshal, σdeep], being overwhelmed by the presence or absence of warm core eddies. On the other hand, is less affected by the presence of eddies, making it much smaller within [σshal, σdeep] but of comparable magnitude above.
Estimation of the mixing length requires knowledge of the mean spice gradient that is presumably stirred. Figure 14c shows time-averaged spice along the 300 line (location of Marguerite Trough) as computed from the historic Pal LTER CTD data. There is a very strong gradient at the continental slope (180 station; gray line in Fig. 14d), however, unlike other grid lines, on the 300 line there is also a secondary region of large cross-slope spice gradient within [σshal, σdeep] near the shoreward 130 station (blue line in Fig. 14d). We suspect that the gradient there is associated with the quasi-permanent diversion of the shelf-break current onto the shelf upon interacting with Marguerite Trough (Fig. 2).
Importantly, it is not known whether the interleaving layers are generated locally or remotely, so it is not obvious which is the relevant gradient to be stirred. Likewise, it is not immediately clear whether the relevant variance is or (i.e., should the high spice cores of the eddies be removed?). A conservative estimate is to assume that the variance is strictly that of the interleaving layers and that the gradient that is stirred is the local gradient within Marguerite Trough. That would suggest = 5.6 km averaged between [σshal, σdeep] (with range 3.4–7.5 km). An alternate perspective is to treat the warm eddy cores themselves as filamentary structures associated with stirring of the larger cross-slope gradient. In that case, = 11.1 km (with range 3.4–17.2 km). Either way, the estimates are of the same order of magnitude as R and |kmax|−1. This suggests that the observed eddies are of sufficient size to generate the observed spice variance by stirring the mean spice gradient.
It is not obvious what causes the interleaving. Internal waves cannot provide enough temperature variance unless the lateral temperature gradient is on the larger end of those observed and the internal wave field is strongly biased toward inertial waves (which it somewhat is with Rω = 9). On the other hand, double diffusion, while qualitatively consistent with the vertical structure of the observations, tends to have slow growth rates. A more likely candidate appears to be stirring by the eddies of the cross-shelf temperature gradient within Marguerite Trough. Spice variance is largest at the eddy edges and its magnitude is consistent with the eddy sizes in accordance with mixing length theory. Thermohaline intrusions appear to be important for the decay of the eddies, particularly in their upper hemispheres, as cold neighboring water is brought inwards and warm slope water is ejected.
c. Frictional spindown
What appears to be cooling of the eddies may instead be redistribution of heat via a change in the buoyancy distribution of the eddy. Under the assumptions of a strictly vertical exchange of both mass and momentum, spindown of a geostrophically balanced along-front (or azimuthal) flow in the ocean interior leads to a secondary radial circulation that redistributes buoyancy so as to flatten isopycnals (Garrett 1982; Bormans 1992b). The potential vorticity equation takes the form of a diffusion equation in which the local time evolution of buoyancy B is balanced by a viscous spreading and a diapycnal diffusion of mass:
Our data are incapable of being used to invert such a model, but we can make a few comments. The geostophic shear we calculate is evaluated at approximately the depth of the eddy core and is relative to the bottom, a depth range where we have previously shown there is strong potential for diapycnal mixing. It is possible that the turbulent mixing itself flattens isopycnals without the need for friction to reduce the azimuthal geostrophic flow. We do not know the Prandtl number within these eddies, but it is generally thought to be of order 1. Assuming that it is 7 (canonical seawater value), an upper bound for a middepth value of the effective viscous spreading coefficient is of order 1 m2 s−1, which is comparable to the lateral diffusivity estimated in the previous section. There is insufficient evidence to confirm or refute any role of spindown in the eddies’ decay.
Here we use a simple diffusion model informed by the results of the previous section to simulate the total heat loss of the eddy we tracked over crossings 1–5. We would like to emphasize that we do not expect to reproduce the small-scale structure of the eddy’s cooling but instead aim to capture the bulk heat loss. The geometry of the first crossing gives us confidence that we passed directly through the eddy’s center and through its entire diameter, therefore we set the initial condition as the hydrographic structure of the eddy as sampled on crossing 1. We identify the eddy’s geographic center (x1, y1) and radius (R1 = 3.4 km) as the mean and standard deviation, respectively, of a 2D Gaussian fit to the field of heat content per unit area integrated between σshal and σdeep. We then assign each profile a radial distance from that center and grid the T and S on each isopycnal to construct a 3D field for the eddy. Outside of the eddy the domain is populated by the mCDW profile.
a. Diffusion model
We treat temperature and salinity as passive tracers and allow the eddy to diffuse in a domain x = [−10 km, +10 km], y = [−10 km, +10 km], and z = [0 m, 480 m] with 250-m grid spacing in the horizontal and 2.5-m grid spacing in the vertical. The equation for temperature is
and it is integrated for 4.5 days with a time step of 150 s using a finite-difference control volume approach. The diffusivity values are constant in time and are based on those parameterized in the previous section. The southern shelf composite Kz profile (Fig. 10) increases log linearly with depth from a value of 4.23 × 10−6 m2 s−1 at 120 m, however we found better results are achieved if we set Ktop to be 5 times that value, which is within the reported 95% confidence interval.
b. Fit to data
We use a diffusion-only equation because the eddy’s geographic coordinates [x(t), y(t)] are additional unknowns. Therefore we assume that the eddy is advected by a depth-invariant current and seek an optimal [x(τi), y(τi)] at each crossing i. We do this by performing a grid search by moving the modeled eddy to each location in a grid encompassing the glider track and then subsampling the modeled field along the glider track and calculating the depth-integrated heat content along the track. The position k that minimizes the root-mean-square (RMS) error between and defines the eddy center’s location at time τi. On each crossing, there are a few profiles where thermohaline intrusions significantly attenuate that profile’s heat per unit area and the presumed symmetry of the eddy and therefore they are excluded in the RMS fit (but are included when evaluating total along-track heat content).
To place error bars on our modeled heat content, we perturb the best fit eddy location by 1.2 km radially, again subsample the modeled field along the glider track, and then compute total along-track (x) and vertically (z) integrated heat content:
The spread in provides a measure of uncertainty in the heat content of the best-fit eddy.
c. Results and interpretation
The best fit locations and modeled heat content are shown in Fig. 15. The predicted eddy track flows counter clockwise about the seamount east of Marguerite Trough (near 300.060), consistent with the mean flow, and the glider’s excursion to the south confirms that this eddy did not go into Marguerite Bay. We are confident that the glider was indeed sampling the same eddy. Simulated particle trajectories (not shown) originating at [x(τ1), y(τ1)] and using all possible 4-day sequences of the time-varying, depth-averaged current fields from Dinniman et al. (2011) confirm that parcels of water can follow the best-fit trajectory and can traverse the proposed distance traveled by the eddy over the total observation time.
In general, there is good agreement between the time series of observed and modeled along-track heat content (Fig. 15b). Particularly, the agreement is good near the edge of the eddy which contributes more to the integral of total heat content. There is a wide variety of intrusions, filimentation, and other submesoscale variability along the track. For example, the fourth crossing appears to skirt the eddy edge but contains a profile of water ejected from the eddy core. The second and third crossings appear to show the eddy core “split” by the cool layer that began penetrating the eddy during the first crossing. Obviously the model does not simulate these processes though it does seem capable of parameterizing their consequences. The total integrated heat contents fall within the range of perturbed model fits (Fig. 15c).
The results can be used to estimate the rate of heat loss of the eddy between crossings 1 and 5. The total heat content of the eddy during crossing 1 is easily obtained by integrating the model initial condition (itself a fit to the data) vertically between σshal and σdeep, and then laterally out to R1. At crossing 5, we obtain Qmax,chord and Rchord by fitting Eq. (2.1) to the data but to calculate the total heat content we need a cross section through the eddy center. Simple trigonometry shows that and , where D is the distance between the true eddy center and the chord-length center. Using D = 2.5 km (Fig. 15a), the total heat content is then obtained from Eq. (2.2) by using Qmax,real and Rreal and integrating out to R1. The initial and final heat contents are Qtot(τ1) = 7.9 × 1015 J and Qtot(τ5) = 5.5 × 1015 J, implying a cooling rate ∂Q/∂t = 7.0 × 109 J s−1. It is worth noting that the modeled heat content on crossing 5 is 5.6 × 1015 J, only 2% greater than observed.
6. Discussions and conclusions
The first goal of this study was to understand the origin of UCDW eddies on the shelf. We conducted a linear stability analysis of the shelf-break current to show that its most unstable mode has a CDW-level amplitude maximum and an inverse wavenumber of 4.4 km. Glider observations confirmed that UCDW eddies on the shelf have diameters and crest-to-crest separation consistent with that metric and their anomalies are confined to the subpycnocline layer. The strong shear and weak stratification allow for sign changes in the interior potential vorticity gradient which cause instabilities smaller than the Rossby radius to persist even when Eady modes are suppressed by the bottom slope. This, in combination with the insensitivity to the planetary beta effect, suggests to us that similar UCDW eddies should be common around the Antarctic margins. However, we must be cautious with any broad extensions. While we have shown that either eastward or westward currents can support eddies, westward currents are often associated with the Antarctic Slope Front and very different stratification. The demonstrated sensitivity of eddy fluxes of CDW to other parameters that were not considered here, such as depth of the continental shelf or prominence of the Antarctic Slope Front (Stewart and Thompson 2015), restrict any generalizations and suggest eddy fluxes are still likely highly localized (Stewart et al. 2018) with one region of importance being the WAP.
The second goal of this study was to identify the major mixing processes that disperse eddy heat into surrounding water. Shipboard CTD and SADCP measurements were used to show cross-pycnocline diffusive fluxes should on average be <1 W m−2 while fluxes below the eddy Tmax should be an order of magnitude larger. CTD and glider profiles reveal that interleaving layers and thermohaline intrusions are ubiquitous, particularly in the lower pycnocline/eddy upper hemispheres where they link warm slope-type UCDW to the cooler waters within the adjacent water column’s broad pycnocline. This water is most quickly eroded over the Advective Path while the deeper, core Tmax at ~250 m persists longer. The erosion of the core heat content per unit area is more apparent than that of the temperature maximum, implying erosion at the eddy boundaries and/or a redistribution of heat within the eddy interior. Consistent with this, high-pass-filtered spice variance is largest at the eddy edges.
The origin of the thermohaline intrusions is ambiguous. The magnitude of their thermal variance is larger than would be expected if due solely to internal wave shear, although this conclusion is highly dependent on the magnitude of the lateral temperature gradient across the eddy. In addition, thermal variance is largest above and just below the eddy core, both regions susceptible to double-diffusive instabilities. Instead, given the observed along-isopycnal spice variance, the mean spice gradient, and the size of the first-mode instability, a more plausible explanation for the origin of the interleaving layers may be stirring by the eddies themselves. A recent modeling study (Stewart et al. 2018) suggests that eddies may transfer heat across the Antarctic slope primarily by along-isopycnal stirring as opposed to advection by their overturning streamfunction. This also supports stirring as an explanation for the thermohaline variance.
The flexibility afforded by Slocum gliders allows for real-time, data-adaptive sampling. One of the novel aspects of this study is the collection of cross sections through an eddy as it traversed the shelf, providing a first estimate of the eddies’ rate of cooling, which is consistent with that of a 3D diffusion model applied to the initial crossing. A lower limit on the time required to eliminate all heat relative to the mCDW profile within one radius from the eddy center is given by Qtot(τ1)/(∂Q/∂t) = 13.1 days which assumes that the gradients are constant in time. Because lateral mixing dominates, an alternate cooling time is given by the solution to the two-dimensional diffusion equation with homogeneous lateral boundary conditions at infinity and a Gaussian initial condition of radius R. That solution implies that the heat content per unit area at the eddy center decreases as R2/(R2 + 4Kht). For the tracked eddy a 50% decrease is achieved after 10.3 days. For a typical eddy (R = |kmax|−1), a 50% decrease takes 17.5 days.
We can synthesize these results with a heat budget integrated over a homogeneous eddy whose dimensions and gradients do not change in time:
In this formulation, ε collects the mismatch between the left and right hand sides. It is meant to represent the effects of noise in the data but also captures the effects of processes that restructure the eddy. Because we only know ∂Q/∂t for the tracked eddy (7.0 × 109 J s−1), we balance the budget with that eddy’s gradients and geometry averaged across crossings 1 and 5. For the tracked eddy, R = 3.4 km, H = 188 m, 〈∂T/∂z|top〉 = 0.019° C m−1, 〈∂T/∂z|bot〉 = −0.0014° C m−1, 〈∂T/∂r|r=R〉 = 1.1 × 10−4 °C m−1, and the diffusivities are those used in the model, evaluated at the average depths of σshal and σdeep. The balance implies the following conceptual model of eddy cooling. About 2% of initial heat loss occurs diapycnally through the permanent pycnocline; about 3% occurs diapycnally through the bottom and 95% occurs laterally. The lateral heat loss, which occurs in part through thermohaline intrusions, directly warms the surrounding pycnocline base and middepth waters. The diapycnal mixing through the bottom also works to flatten isopycnals at and below the Tmax which reduces the azimuthal geostrophic current. The mismatch ε/(∂Q/∂t) = 10% is small, suggesting redistribution processes such as viscous spreading are of secondary importance compared to lateral mixing. To place this number in context, if the lateral diffusivity was increased from 3.2 to 3.6 m2 s−1 (which is within one standard error) then ε goes to zero. The budget is shown schematically in Fig. 16 with fluxes estimated for a more typical eddy (gradients and geometry from Table 1).
Importantly, a preference for the eddies to mix laterally and downward suggests that the eddies are good at redistributing heat rather than immediately venting it to the atmosphere. This sheltering has implications for the ability of intra and subpycnocline heat to persist and make its way shoreward toward marine-terminating glaciers even as the eddies themselves cease to remain coherent structures.
McKee and Martinson were supported by NSF Grants OPP-08-23101 and PLR-1440435; Schofield was supported by NSF Grant ANT-08-23101 and NASA Grant NNX14AL86G. Richard Iannuzzi processed the CTD and glider data, made Fig. 1, and provided valuable comments. We thank David Aragon and the Rutgers COOL group for deploying the glider and providing logistical support. We also thank Pat Caldwell at the JASADCP for culling and providing the SADCP data and Mike Dinniman for providing model current fields. We also acknowledge the captains and crews of the L.M. Gould and the N.B. Palmer for their efforts in acquiring many years of LTER data. Insightful comments from the editor and two anonymous reviewers greatly improved this manuscript. This is LDEO contribution 8310 and Palmer LTER contribution 620.
The Garrett–Munk Spectrum on the WAP Continental Shelf
The Garrett–Munk (GM; Garrett and Munk 1975) spectrum provides an empirical description of the internal wave field in a flat-bottom open ocean away from boundaries and is frequently used to provide a statistical measure of internal wave variability in mixing parameterizations of nonlinear wave–wave interactions (e.g., G89). The GM spectrum ignores the presence of upper and lower boundaries and vertical structure in the wave field is manifested through a series of modal structures. Two fundamental parameters that scale the spectrum are a nondimensional energy E (canonically 6.3 × 10−5) and a stratification depth b (canonically 1300 m). The GM shear spectrum is
where β is the vertical wavenumber (rad m−1), is a reference wavenumber, N0 = 0.0052 rad s−1 is the canonical reference stratification, and j* = 3 is the canonical mode number. A roll-off of −1 is imposed at wavenumbers greater than βc = 0.2π rad m−1.
On the WAP continental shelf with typical bottom depths of ~400 m, the nondimensional energy E and the stratification scale b are likely quite different from canonical values (Levine 2002). Here we maintain the GM functional form (thus neglecting vertical boundaries) but adjust E and b in accordance with the environment of the WAP. One way to do so (section 4c of Levine 2002) is to first estimate the product N0b from the observed stratification N(z) and then to estimate the energy density in the internal wave band. We estimate N(z) by first computing an average density profile from all summertime (January and February) Shelf CTD casts. For all casts, the depths of a set of evenly spaced isopycnals are calculated via linear interpolation and each isopycnal is assigned its mean depth. The resulting potential density profile is then interpolated to an 8-m grid and finite-differenced in order to obtain a profile of mean N(z). The product N0b, which scales the wavenumber bandwidth, is defined as
so that b = 288 m. The profile of mean N squared is shown in Fig. 9a. The lower integration bound is chosen as an average depth of all shelf stations and the upper integration bound is chosen as a typical seasonal mixed layer depth since overturns in the averaged density profile become apparent shallower than that [because the depth of the seasonal mixed layer varies regionally and interannually, the mean N(z) profile is quasi-exponential with no apparent seasonal mixed layer]. Parameters N0 and j* are left at their canonical values.
To estimate the nondimensional energy density, we first estimate the kinetic energy density KE via baroclinic energy spectra from the moored current meter data of Martinson and McKee (2012). Their current meter data generally sample once per hour, however in 2007 they have one year-long record at site 300.100 that sampled once every 20 min with two current meters spanning the permanent pycnocline (one each nominally at Tmin and Tmax). We define the barotropic current as the average of the two records and subtract that out. Then, the spectrum of the residual baroclinic current at each depth is calculated via squaring the Fourier transform of the complex time series, is averaged in the vertical, and is integrated over [f, N0]. While the Nyquist frequency for this sampling is 0.5N0, the −2 roll-off in the internal wave band reveals that most of the contribution to the integral comes from the lower portion of the band. The depth average of the two integrated spectra times 1/2 yields the kinetic energy density (KE) and the nondimensional energy density is given by
We now consider how the modified GM spectrum translates to parameterized diffusivities. Dissipation in the G89 model is given by
where the constant ε0 is derived from wave–wave interactions in the GM energy flux spectrum. The assumption of a steady turbulent kinetic energy balance and a constant mixing ratio yields the relation for the diffusivity
as presented in the text [but here without the corrections h and j; Eqs. (4.4) and (4.5)]. Under the canonical GM spectrum, K0 = 5.0 × 10−6 m2 s−1 whereas for the local parameters derived above K0 = 5.6 × 10−6 m2 s−1, which is a small difference. In addition, compared to shear variances using canonical GM parameters, integrated shear variances using the corrected GM parameters yield diffusivities that are generally up to a factor of 2 greater in the upper water column but almost equal near the Tmax. Given that the G89 method is demonstrated to be accurate within a factor of 2, all of this suggests that the improvements afforded by the corrected GM spectrum are small.