Ocean ventilation is the process by which climatically important tracers such as heat and carbon are exchanged between the atmosphere and ocean interior. In this paper a series of numerical simulations are used to study the interaction of submesoscales and a topographically steered jet in driving rapid ventilation. The ventilation is found to increase both as a function of wind stress and model resolution, with a submesoscale-resolving 1/120° model exhibiting the largest ventilation rate. The jet in this simulation is found to be persistently unstable to submesoscale instabilities, which are known to feature intense vertical circulations. The vertical tracer transport is found to scale as a function of the eddy kinetic energy and mean isopycnal slope, whose behaviors change as a function of the wind stress and due to the emergence of a strong potential vorticity gradient due to the lateral shear of the jet. These results highlight the importance of jet–submesoscale interaction as a bridge between the atmosphere and the ocean interior.
The oceanic overturning circulation plays a governing role in Earth’s climate because of the enormous capacity of the ocean to store heat and CO2. The Southern Ocean, which surrounds Antarctica, plays a disproportionate role in this overturning circulation because it is one of the main areas where deep waters rise to the surface to exchange heat and CO2 with the atmosphere (Le Quéré et al. 2007; Marshall and Speer 2012). This exchange process, known as Southern Ocean ventilation, therefore controls the uptake and storage of heat and CO2 into the global ocean. Changes in Southern Ocean ventilation have been closely linked to abrupt climate transitions in Earth’s past, but the driving mechanisms for this are not understood. As a consequence we do not know the potential for abrupt climate change under anthropogenic forcing in a future climate.
The rate-limiting step in the exchange of tracers between the atmosphere and the deep ocean is the exchange of tracers across the base of the mixed layer, which is strongly modulated by the presence of submesoscale eddies. Submesoscale eddies are ubiquitous features in the oceanic boundary layers, and exert a strong control over the boundary layer shear and stratification at fast time scales. The vertical circulations associated with submesoscale instabilities and turbulence are well appreciated for their role in modulating the biogeochemistry and ecology near the ocean surface (Mahadevan 2016; Lévy et al. 2018), and for their dynamical role in establishing a “bridge” between the atmosphere and the oceanic interior (McWilliams 2016; Balwada et al. 2018). Tracers and gases crossing this bridge are able to be subducted or entrained rapidly, making submesoscale effects a key contributor to global ocean–atmosphere exchange despite their transient nature and intermediate [ km] spatial scales (Siegelman et al. 2020).
The enormous variety in the behavior of submesoscale processes has generally necessitated studying each of them in isolation [e.g., idealized studies such as Boccaletti et al. (2007), Bachman and Taylor (2014), and Jiao and Dewar (2015), which focus in individual submesoscale instability types]. However, the groundwork laid by these early studies has allowed progress toward more nuanced topics, such the interactions between these processes (Stamper and Taylor 2017; Stamper et al. 2018) and with atmospheric forcing (Bachman et al. 2017a; Callies and Ferrari 2018). One such topic which has recently received attention concerns the effects of submesoscale instabilities and waves in the presence of a lateral potential vorticity (PV) gradient (Stamper et al. 2018; Taylor et al. 2018). In this context the PV gradient is generally conceptualized as being due to a laterally sheared jet associated with a mesoscale density front in thermal wind balance. Taylor et al. (2018) show that submesoscale instabilities are capable of growing and rapidly subducting surface water along the front even though their nonlinear evolution is arrested by the presence of the jet. It is suggested that these types of interactions between submesoscales and jets, despite their sporadic occurrence and transience, may contribute a nonnegligible fraction of total subduction when accumulated globally.
It remains an open question whether submesoscales are required for amplified subduction along jets, or if the jets themselves can induce strong subduction by permitting deep isopycnal layers to outcrop. The latter possibility was raised by Klocker (2018, hereafter K18), who showed that when the winds are sufficiently strong the flow can undergo a transition from a state dominated by turbulence to a state dominated by jet dynamics. In the turbulence-dominated regime baroclinic eddies were capable of restratifying the surface boundary layer to prevent outcropping of deep isopycnal layers, essentially keeping a “window” from the surface to the ocean interior shut. In the jet-dominated state the jet itself establishes a potential vorticity gradient that inhibits the formation of the baroclinic eddies. In turn, this allows the mean isopycnal slope within the jet to steepen and generates a larger Ekman buoyancy flux that destratifies in the surface boundary layer. This effect is so strong that it can create an unstably stratified mixed layer that connects with the isopycnal layers beneath the pycnocline, creating the conditions for greatly enhanced ventilation rates.
Indeed, the simulations conducted in K18 showed a massive, nonlinear increase in the ventilation rate as the wind stress increased. Because these models were run at 1/20° resolution it was assumed that the grid scale was too coarse to resolve submesoscale effects, and that the increase in ventilation was driven primarily through the outcropping of the deep isopycnals. However, it is possible that the resolution of this model was sufficient to marginally resolve larger submesoscale dynamics, so it remained unclear to what extent the jet and the submesoscales conspired to produce the ventilation increase.
In this paper we aim to answer the question of why such amplified tracer ventilation was observed as the wind stress increased in K18, and to assess which mechanism—submesoscale dynamics or the jet—was chiefly responsible. To avoid the ambiguity about whether submesoscale dynamics were present in the K18 simulations, here we conduct a much higher-resolution, submesoscale-resolving simulation at 1/120° using the same atmospheric forcing but in a smaller subdomain of the K18 Kerguelen Plateau simulation. With this new model we intend to provide context for the results of the K18 simulations and explain with greater nuance the interaction between submesoscales, jets, and atmospheric forcing. Perspectives on the role of these processes and their parameterization in climate-scale ocean models are discussed in section 5.
2. Comparison of simulations at different model resolutions
It is well known that the dynamics resolved in numerical ocean models are highly dependent on model resolution, and thus so are important properties of the flow such as stratification, tracer transport, and vertical circulations (Kirtman et al. 2012; Newsom et al. 2016). Several effects of increasing resolution are documented in the literature pertaining to submesoscale eddies, such as the flattening of kinetic energy spectra (Capet et al. 2008; Barkan et al. 2017), increase in RMS values of vorticity, horizontal convergence and strain rate (Rocha et al. 2016; Choi et al. 2017), and change in boundary layer depth (Marchesiello et al. 2011; Bachman et al. 2017b). As grid resolution is increased the resolved dynamics may undergo abrupt regime changes—for example, a series of Southern Ocean simulations conducted by Bachman et al. (2017b) at different resolutions showed sharp increases in RMS vertical velocity upon crossing the 1/96° (submesoscale permitting) and 1/192° (internal gravity wave permitting) thresholds. Indeed, the Kerguelen Plateau simulations conducted in this analysis (see Fig. 1 and appendix A) hint at similar behavior as the resolution is increased (Fig. 2, left column), where the transition from 1/4° (KERG4) to 1/120° (KERG120) resolution sees the emergence of small vortices, filaments, and a conspicuous quasi-zonal jet skirting the northern edge of the plateau (hereafter referred to as “the jet”). This noteworthy change in behavior is why the resolution dependence of the tracer ventilation will be our starting point for investigating the results seen in K18.
As in K18, the wind and buoyancy forcing in these simulations are constant in time and only vary in the meridional direction (see K18, his Fig. 1b), and are calculated as the annually and zonally averaged values from the year 2005 from the Southern Ocean State Estimate (Mazloff et al. 2010). The wind direction is purely zonal, and the meridional wind stress profile obtained via the zonal- and time-averaging is referred to as τ. In K18 additional experiments were run by varying the strength of the wind stress from 0.5τ to 2.0τ, and much of the analysis concerned the stronger wind stress experiments due to the emergence of the jet. Here we also focus the most attention on the experiments with 2.0τ wind stress, but we will analyze the weaker wind cases later in the manuscript. We will be explicit about the wind strength shown in the figures either through the figure captions or labels.
K18 notes that the jet and the underlying oceanic density structure undergo a striking transition as the wind stress is increased, marked by an increase in the time-mean isopycnal slope, destratification of the surface mixed layer, and a direct outcropping of deep isopycnal layers (see K18, his Fig. 5), all of which accompanied an abrupt and significant increase in ventilation. The behavior we observe at 2.0τ wind stress is generally consistent with K18 across different model resolutions (Fig. 2). As the resolution is increased the jet migrates slightly southward toward Kerguelen Island, but the eastward-flowing jet core remains distinct from the surrounding fluid in all cases (center column). The jet itself becomes slightly tighter at higher resolution and begins to develop more complicated structure, including the appearance of multiple, weaker retrograde jets in the KERG20 and KERG120 cases. The buoyancy frequency (right column) also behaves fairly similarly between the KERG4, KERG10, and KERG20 cases. In each case there is a direct outcropping of isopycnal layers that extend beneath 500-m depth within the jet (white lines), though the packing of layers tends to be tighter and their extent deeper at higher resolution. The stratification in the fluid surrounding the jet is marked by a shallow mixed layer within the top 100 m overlying a pronounced pycnocline. Within the jet itself, however, the stratification is weak at all depths and the pycnocline is absent, and persistent convectively unstable conditions even appear within the top 200 m (red patches). It is noteworthy that this does not occur in the KERG120 case—generally, the pycnocline is shallower and weaker in all regions of the simulation, and no convectively unstable regions appear even though the isopycnal layers remain nearly vertical within the jet. KERG120 is also the only case where mixed layer depth (green line) does not experience a pronounced deepening along the northern flank of the jet. These differences between KERG120 and the lower-resolution cases indicate some emergent dynamics at very high resolution that exert a strong control over the ocean density structure, and warrants an investigation into their effects on tracer ventilation.
Figure 3 shows time series of the change in mean tracer concentration between 500- and 1500-m depth, 49°–39°S, and 61°–84°E. These bounds were chosen to match the depth range used in Fig. 9 of K18, but also were restricted to the smaller domain size used in KERG120 while avoiding the sponge layer on the lateral domain boundaries. In this figure it is clear that the model resolution has a significant impact on the amount of tracer that can be subducted into the ocean interior, where the steady concentration increases by over an order of magnitude from 1/4° resolution (7.29 × 102 m−3) to 1/10° (1.96 × 104 m−3), by nearly a factor of 4 from 1/10° to 1/20° (7.52 × 104 m−3), and by another 40% from 1/20° to 1/120° (1.04 × 105 m−3). In conjunction with Fig. 2 and the results of K18, this suggests that the outcropping of isopycnals and weak mixed layer stratification are necessary but not sufficient conditions for tracer ventilation. That is, the dependence of tracer concentration on model resolution implies that intermediate-scale processes, herein interpreted as submesoscale dynamics (5), play a central role in driving tracer subduction at the jet.
In the following sections we will examine the conditions that lead to submesoscale instability in the jet and how the jet itself modulates submesoscale turbulence and energy. We will then return to the perturbation experiments from K18, wherein the wind stress is varied in several models at 1/20° resolution, to give a more nuanced interpretation of the role of submesoscale and jet dynamics in driving ocean ventilation.
3. Submesoscale instabilities in the KERG120 simulation
The ocean surface boundary layer (SBL) lies at the contact point between the ocean and atmosphere, and thus the stress and buoyancy forcing from the atmosphere exert a heavy influence on the processes occurring therein. Boundary layer turbulence is the principle by-product of the atmospheric forcing, serving to transport momentum, heat, and tracers through a well-mixed layer which is often (somewhat imprecisely) considered as analogous to the boundary layer. Here we remind the reader that, in the context of submesoscale turbulence, it is important to keep the concepts of the mixed layer and boundary layer separate. That is, the notion of a boundary layer arises because the fluid stress must be continuous at the air–sea interface. Because the ocean and atmosphere move relative to each other, the stress matching results in a number of frictional processes that affect the shear and stratification on both sides of the interface (Thomas and Ferrari 2008). The resultant surface waves and turbulence tend to keep the layer weakly stratified (hence “mixed” layer), but dynamically this is not constrained to be the case (Callies and Ferrari 2018; Young 1994; Bachman and Taylor 2016).
Both the wind and surface heat loss are well-known drivers of turbulence in the SBL (Shay and Gregg 1984). The frictional force associated with the wind drives Ekman flow that is orthogonal to the wind direction and is given by , where is the vertical unit vector and f is the local Coriolis parameter. Ekman flow that is directed across surfaces of constant buoyancy, as is the case when the wind is directed along geostrophically balanced fronts (i.e., parallel with the geostrophic shear), advects lighter water over dense in the case of “up-front” winds (wind directed against the shear) or denser water over light for “down-front” winds (wind blowing with the shear). In the latter case a variety of submesoscale instabilities are excited (Thomas et al. 2013). In all, the dynamics that are driven by the interplay between the wind and oceanic density field are a function of the wind direction, wind strength, and oceanic frontal strength, all of which rapidly evolve in time and space.
The bathymetry of the Kerguelen Plateau constitutes a barrier to the flow of the Antarctic Circumpolar Current, and is responsible for a significant standing meander that skirts the northern edge of the plateau (Damerell et al. 2013). The eastward-flowing jet associated with these fronts also follows the bathymetry of the plateau (Park et al. 1993; Sparrow et al. 1996), as seen in Fig. 2, flowing in a predominantly zonal direction as it passes by its northernmost edge. The westerly winds in the KERG simulations are thus roughly aligned with the jet and are directed down-front, resulting in Ekman advection of cold, dense fluid from the southern side of the front to the warmer, lighter north side. The time-mean Ekman buoyancy flux (EBF), given by
is a way of measuring this transport, where |∇hb| is the lateral buoyancy gradient and ρ0 is a reference density (Fig. 4a). Positive values reflect the destabilizing nature of this transport, indicating a tendency to generate convective turbulence and destratify the boundary layer.
The dynamical environment set up by the wind, jet, and density front thus precondition the flow to a variety of dynamical instabilities that occur due to weak stratification or large shear. Several types of submesoscale instability have been identified which are capable of significant vertical transport at small time scales. Such instabilities are often described based on the condition fq < 0, where f is the local Coriolis parameter and q is the Ertel potential vorticity (PV). Specific instabilities can be identified by writing the PV as
where N2 = ∂b/∂z and ω is the absolute vorticity vector split into its vertical (ωυ = f + ζ) and horizontal (ωh) components. Assuming a flow in geostrophic and hydrostatic balance, one may leverage the thermal wind relation,
where Ro = ζ/f is the vorticity Rossby number and Ri = N2f2/|∇hb|2 is the balanced Richardson number. The expression (4) reveals that the condition fq < 0 can be met in multiple ways: N2 can be negative, indicating convective instability; sufficiently strong anticyclonic vorticity, associated with negative values of Ro and indicating centrifugal instability (CI); and strong baroclinicity combined with weak but stable stratification, reflected in small (positive) values of Ri and associated with symmetric instability (SI). We also will consider a fourth type of submesoscale instability, ageostrophic anticyclonic instability (AAI; McWilliams et al. 2004; Wang et al. 2014), which is induced when the horizontal strain rate, , is strong enough that . We acknowledge that realistic flows rarely satisfy the formal requirements by which these instabilities are derived (linearized, geostrophic, hydrostatic, etc.) and that instabilities are often mixtures of the above types (e.g., Stamper and Taylor 2017; Grisouard 2018); nonetheless, the linear framework often proves to be a useful means of examining submesoscale dynamics in models and observations (e.g., Adams et al. 2017). Note that this analysis excludes mixed layer baroclinic instability (e.g., Boccaletti et al. 2007; Fox-Kemper et al. 2008), which can occur even when fq > 0.
The time-mean Ri is shown in Fig. 4b and affirms that areas of strongly positive EBF are associated with weakened stratification and lower values of Ri. The jet region features particularly small values of Ri in comparison with the rest of the domain, though the time-mean values are well above what would be required for fq < 0 [for a stably stratified flow with Ro = 0 the flow would be unstable for Ri < 1, whereas the jet region features time-mean ]. This occurs because the density fronts and filaments at which Ri ~ 1 are in reality much smaller in scale and highly transient, as can be seen in the instantaneous snapshot in Fig. 4c. The domain is filled with such features, which also tend to occur at the fringes of mesoscale eddies and filaments (Brannigan et al. 2017). The low time-mean Ri in the jet thus indicates that these features are much more likely to occur there and would tend to be highly concentrated.
Figure 4 thus suggests that, at any given time, individual density fronts within the jet would be likely to undergo submesoscale instability, and that collectively these instabilities would occur quite often. To investigate this hypothesis, the model output was used to measure the frequency at which the flow is unstable to the instabilities mentioned previously. For each output snapshot the values of Ri, Ro, ωυ, and S were diagnosed at each grid point and the conditions for instability were checked (see appendix B). Figure 5 shows maps resulting from this procedure and confirms that the jet region is persistently unstable to all types of instability, albeit at different frequencies. CI occurs most frequently (Fig. 5b), at around twice the frequency of AAI (Fig. 5c) and 5 times the frequency of SI (Fig. 5a). Altogether, the flow within the jet is unstable to some form of these submesoscale instabilities around 25% of the time.
Submesoscale turbulence is also well known for its role in the oceanic energy budget, which is to act as a conduit to transfer energy from the large-scale mean circulation down to dissipative scales (McWilliams 2016). Mathematically this transfer takes place through energy conversion terms in the averaged kinetic and potential energy budgets, namely, through the geostrophic shear production (SI), lateral shear production (CI), vertical buoyancy flux (from convective and mixed layer baroclinic instabilities, e.g., Boccaletti et al. 2007), and combinations of the above (AAI). Diagnoses of these energy conversion terms from the model output gives an indication of the types of instability that feed the turbulence, and are shown in Fig. 6. The notation is such that overbars indicate a time average over 1 year of model output, and primes indicate a deviation from this average derived from instantaneous snapshots (see also appendix A).
The lateral shear production in Fig. 6b is the dominant energy source for turbulent kinetic energy, which is consistent with CI being the most common form of submesoscale instability in Fig. 5. The vertical buoyancy flux (Fig. 6d) is also significant, and is most commonly associated with standard, mixed layer baroclinic turbulence, though buoyancy production can also be generated by SI, CI, and mixed instabilities (Grisouard 2018). The geostrophic and ageostrophic shear production (Figs. 6a,c) are of minor significance in the jet.
In summary, the wind stress and jet conspire to create a flow that is frequently unstable to several forms of submesoscale instability. The density front underlying the jet creates a large potential energy reservoir that these instabilities can draw upon, as does the lateral shear associated with the jet itself. As noted in K18, the shear in the stronger wind stress cases also creates a potential vorticity gradient that partially suppresses mixed layer turbulence. A combination of this effect and persistent, destratifying EBF prevent baroclinic turbulence from tilting the density front over and cutting off ventilation to the ocean interior. Thus, in the jet there exist both a “window” for ocean ventilation (weak stratification and a direct connection to deep density layers) and a powerful mechanism for vertical transport (submesoscale instabilities), and we now proceed to analyze how these mechanisms modulate tracer transport in the KERG20 and KERG120 simulations.
4. Eddy energy and isopycnal slopes modulate ventilation rates
The proliferation of submesoscale instabilities within the jet lead to significant variability at fast time scales. K18 notes that this variability, which was measured as the eddy kinetic energy in his Fig. 5c, is tempered somewhat in the stronger wind stress simulations by the emergence of a potential vorticity gradient associated with the jet’s lateral shear. However, these simulations also exhibited progressively greater tracer ventilation, creating a somewhat paradoxical scenario where weaker variability corresponded to stronger transport. We now return to the KERG20 simulations from K18 to revisit those results from the joint perspective of submesoscale and jet dynamics.
The roles of these processes in driving subduction can be understood by analyzing the variability and stratification within the jet. Figure 7 shows vertical structures of several key variables, which are averaged in time and meridionally from 45.5° to 44.5°S at 70°E (along the green line in Fig. 1), corresponding to the width of the jet. Figure 7a shows a similar result to that observed in K18 (his Fig. 7c), where the eddy kinetic energy,
steadily decreases as the wind stress increases in the KERG20 simulations. The EKE from the KERG120 simulation (black line) is also shown for comparison, and, naturally, is significantly larger due to its greater resolution of submesoscale turbulence. Here we note that, in general, the increase or decrease of EKE with increasing model resolution depends on the types of dynamics that become resolved on the high-resolution grid and the nature of the averaging operator. Since here we are considering the dynamics of the jet, which is relatively narrow and largely free of mesoscale turbulence and large eddies (Fig. 2j), the emergence of submesoscale features is the main difference between these two models and their resultant EKE fields.
The behavior of the EKE in Fig. 7a is at odds with the notion that submesoscale turbulence causes greater ventilation at higher wind stresses, meaning that other aspects of the flow must be taken into account as well. To this end, Figs. 7b and 7c show the time-averaged vertical and lateral buoyancy gradients, or and , respectively. As stated previously, interacts with the wind stress to produce a destratifying EBF which exerts a strong control on , so it is useful to consider first. Figure 7c shows that stronger wind stresses tend to increase the magnitude of the lateral buoyancy gradient, though the values differ between all of the KERG20 simulations only by around 20%. As seen with the EKE, tends to be stronger in KERG120 due to tighter fronts and filaments being permitted at higher resolution, an effect which has been noted previously (Bachman et al. 2017b; Fox-Kemper et al. 2011). The increased EBF at higher wind stresses in turn produce lower stratification in Fig. 7b. It is noteworthy, however, that the behavior in KERG120 does not follow this pattern—the stratification in this simulation is significantly stronger than that of the 2τ KERG20 simulation, despite having far larger .
The combined effects of the buoyancy gradients are encapsulated by the time-mean isopycnal slope, , which is shown in Fig. 7d. As before, there is a clear trend in the KERG20 simulations, with stronger slopes corresponding to stronger wind stresses. The slope in KERG120 is much weaker than in the 2τ KERG20 simulation due to its stronger vertical stratification.
The issue now is how to synthesize the information from each of these panels to explain the increased ventilation at higher wind stresses in KERG20, as well as the even stronger ventilation yet in KERG120. Figure 7e indicates that the square of the eddy vertical velocity follows this pattern, with greatly increased velocities as the wind stress increases and an even larger jump in KERG120 (note the logarithmic scale). We approach this problem by seeking a scaling for this vertical velocity, starting with the continuity equation for the eddy velocities,
This equation suggests a scaling w′ ~ u′(H/L), where H and L are vertical and horizontal length scales, respectively. Here we follow the reasoning of Fox-Kemper et al. (2008) and assume that parcel displacements occur along paths that are proportional to the mean isopycnal slope, such that ratio of these length scales is H/L ~ |S|. We note that this scaling leaves room for a nondimensional prefactor, and because the dynamics considered here are not restricted to baroclinic mixed layer eddies or their buoyancy fluxes, we do not necessarily expect this prefactor to match those discussed in Fox-Kemper et al. (2008).
Altogether the scalings above lead to an expectation that
for a prefactor C to be determined. We empirically determined the optimal value for C to be 0.06, which optimizes the fit between the diagnoses of and its scaling shown in Fig. 7f. In this panel each quantity is averaged over the mixed layer depth of the respective simulation, which is indicated by the colored hash marks on the left side of Fig. 7a. The scaling (+ symbols) shows good agreement with the diagnosis (circles) across all values of the wind stress, including the KERG120 simulation. For reference we also include the mean tracer concentrations within the KERG120 domain for each simulation (squares), whose ordinate is shown on the right side of the panel and which also trends similarly to the vertical velocity.
The cumulative results from the suite of simulations presented in this manuscript allow us to make the following determinations about the dynamics driving tracer ventilation. Direct outcropping of deep isopycnals can occur even in relatively coarse-resolution models (Fig. 2), but the very high sensitivity of the interior tracer concentration to model resolution (Fig. 3) suggests that submesoscale processes play a key role. Reexamination of the KERG20 simulations confirms that the results in K18 were, in fact, a subtle combination of submesoscale and jet dynamics. The transition of the jet from a baroclinic to a barotropic regime in K18 was accompanied by the emergence of a potential vorticity gradient that suppressed EKE, but also coincided with a large increase in EBF that destratified the upper ocean and steepened the isopycnal slope. The steepening more than compensated for the reduced EKE, and led to greatly increased ventilation rates in the strongest wind stress cases. These interpretations were made possible by the additional KERG120 simulation introduced here, which had reduced isopycnal slopes but massively increased EKE due to the added resolution. This simulation had greater ventilation than all of the KERG20 cases, implicating submesoscale processes as the drivers of the tracer subduction and the isopycnal slope as simply modulating the subduction rate.
The mechanisms of Southern Ocean ventilation remain an important topic in the study of Earth’s climate system, and have spurred significant interest from ocean modelers and observationists alike. The simulations from K18 provided an in-depth analysis of one such mechanism—eddy suppression by a barotropic jet—that highlighted jet dynamics as a possible new frontier in the quest to understand ventilation. Those simulations were advantaged by their sufficiently high resolution to resolve the eddy momentum fluxes that sharpened and maintained the jet, but were disadvantaged by their insufficient resolution to unambiguously resolve submesoscale processes. Given the spate of recent research on intense subduction driven by submesoscales, this paper was conceived as a revisitation of K18 to understand whether submesoscales played a role in those results.
In this manuscript submesoscale instabilities were found to have played a principal role in the ventilation observed by K18. This conclusion was made possible by a suite of new simulations spanning resolutions from 1/4° to 1/120°, which showed that ventilation was extremely sensitive to the model resolution despite the simulations having similar (weak) stratification profiles. The 1/120° simulation, labeled KERG120, featured a jet that was found to be persistently unstable to a variety of submesoscale instabilities despite its strong lateral shear and potential vorticity gradient. These dynamics recall recent work by Taylor et al. (2018), who found that such instabilities can lead to elevated, sustained subduction along barotropic jets even though their nonlinear evolution is arrested.
The 1/20° “KERG20” simulations of K18 were then reexamined for a relationship between the eddy energy and ventilation rate. The key finding in this reanalysis was that the ventilation rate trends with the eddy vertical velocity, which was found via kinematic arguments to scale with the product of eddy kinetic energy and isopycnal slope. The jet and submesoscales were thus determined to play dual (but not necessarily complimentary) roles in the tracer ventilation, i.e., the barotropization of the jet leads to tighter lateral density gradients and steeper slopes, but weakened submesoscale energy due to its potential vorticity gradient. In the KERG20 simulations the former overwhelmed the latter and led to massively increased ventilation. It was only after factoring in KERG120, with its shallower slopes but far greater eddy energy, that this nuance became apparent.
This study adds to the ever-growing body of research highlighting the importance of submesoscale processes as a bridge between the atmosphere and oceanic interior. It also highlights the need for skillful submesoscale eddy parameterizations, particularly for climate-scale ocean models whose grids are far too coarse to resolve even the largest submesoscale dynamics. With regard to the topic of parameterizations, the simulations presented here show that it is possible to resolve a rich field of fronts, filaments, and jets in the submesoscale gray zone (i.e., the 1/10°–1/20° range that is presently considered state-of-the-art for global ocean modeling) while still missing the intense vertical circulations that dominate tracer subduction. Our results also hint that submesocale-resolving “truth” simulations, which are employed in the classical approach of diagnosing turbulent statistics for parameterization development, may need to be well under subkilometer resolution to fully capture all the relevant vertical circulations. Naturally, these issues are beyond the scope of this paper, but reaffirm the benefit of multiscale modeling approaches such as are employed here.
We thank Laure Zanna for useful discussions that helped to inspire this manuscript. This material is based upon work supported by the National Center for Atmospheric Research (NCAR), which is a major facility sponsored by the National Science Foundation under Cooperative Agreement 1852977. We would like to acknowledge high-performance computing support from Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation. This research was supported under Australian Research Council’s Special Research Initiative for Antarctic Gateway Partnership (Project ID SR140300001), and was undertaken with the assistance of resources from the National Computational Infrastructure, which is supported by the Australian Government.
The Massachusetts Institute of Technology general circulation model (Marshall et al. 1997) was used to conduct a series of ocean simulations of the Kerguelen Plateau region at different resolutions. The advection schemes, equation of state, eddy parameterizations, and model domain for the 1/4°, 1/10°, and 1/20° simulations are the same as was used in K18. The 1/120° simulation was conducted in a subdomain of the K18 models, and extended from 60° to 85°E in the zonal direction and from 50° to 38°S in the meridional direction. This domain was sufficient to resolve the Kerguelen Islands and the primary topographic jet that was the focus of the study by K18. The resolution of this simulation was 1/120°, or about 590 m per grid cell in the zonal direction at the southern boundary and 720 m at the northern boundary. The vertical grid consisted of 150 layers of varying thickness, ranging from 10 m at the surface to 50 m at depth. The 1/4° and 1/10° simulations were run for 15 years, the 1/20° for 10 years, and the 1/120° for 5 years. The horizontal eddy viscosity for each simulation was scaled according to the grid resolution.
The model setup for the 1/20° simulations is as described in K18. The model setup for the 1/4° and 1/10° simulations is identical with that for the 1/20° simulations apart from the change in horizontal resolution. No parameterization of mesoscale eddies was used in any of these simulations. For the 1/120° simulation the model bathymetry was interpolated from the global GEBCO bathymetry product (www.gebco.net), which has a resolution of 15 arc s. The maximum depth of the domain was set to be 5000 m. Open boundary conditions were used to force the model at the lateral boundaries, where the velocity, temperature, and salinity fields are forced using daily output from the larger K18 domain. A 1/2°-wide sponge layer was used to relax the model to the boundary conditions, with a 1-day relaxation time scale at the inner edge of the sponge and a 4-h time scale at the outer edge of the sponge. The wind and buoyancy forcing was the same as was used by K18, and was applied as meridionally varying but constant-in-time forcing over the model domain.
The 1/120° simulation was initialized by interpolating output from year 20 of the K18 simulation to the higher-resolution model grid used here. The model was spun up for 1 year, which we deemed to be sufficiently long for the dynamics to equilibrate given the fast [ day] time scales for the submesoscale field to evolve. After the 1-yr spinup, a passive tracer was introduced at the ocean surface over the next 5 years of model integration. The tracer concentration was relaxed to a value of 1 at the surface everywhere except for within the sponge layers, where it is instead relaxed to zero at all depths. A relaxation time scale of 2 h was used, which allowed for a near-instantaneous but numerically stable tracer forcing to the surface (one) and sponge-layer (zero) values. Since the time-mean flow in the domain strongly tends to travel in a northeasterly direction, forcing the tracer in this way allows us to examine the subduction of the tracer due to submesoscale and frontal processes within the domain without the tracer concentration piling up in the interior.
Instantaneous snapshots of the tracer concentration, velocities, potential temperature, and salinity were written out every 6 h for 1 year of simulated time. The 6-h frequency was chosen due to the fast growth time scales for the submesoscale instabilities, as well as the speed (i.e., short advective time scale) of the jet; even 1-day snapshots were found to be too long to capture much of the fast dynamics, and time-averaging was found to smooth the output to an unsatisfactory degree. In all, this amounted to 1440 snapshots of each field, which were sufficient to infer both the time-mean behavior of the jet and dynamical instabilities therein. The notation throughout the article is such that overbars indicate a time average over 1 year of model output, and primes indicate a deviation from this average derived from instantaneous snapshots.
We note that the definition of the time-mean isopycnal slope, , involves averaging the magnitude of the lateral buoyancy gradient vector, , which is a nonlinear term. This average is taken over motions acting on both short and long time scales, which are coupled inside the vector magnitude operator. It is thus a pertinent question whether averaging over a large number of instantaneous fields, in which small-scale fronts and flow structures proliferate (seen, for example, in Fig. 4c), gives a different result than taking the magnitude of the time-averaged buoyancy (Fig. 4b). That is, does |S| represent the average slope of the jet front, or is it biased more to represent the slope of smaller individual fronts within the jet? In the analysis we experimented with both approaches and found that both give quite similar results for |S| (not shown). This is due in part to the large number of realizations used in the time-averaging, where 1440 snapshots effectively smoothed out any effects of individual, strong flow structures. We also note that space and time are highly correlated in these dynamical systems, where long time scales generally correspond to large spatial scales, and vice versa. Thus, a long time average behaves akin to a spatial average over large scales, and we assert that the large number of snapshots used in the averaging means that |S| measures the slope of the large-scale jet front.
Calculating the Frequency of Submesoscale Instability
The procedure used to diagnose the results in Fig. 5 is as follows. For each output snapshot the values of Ri, Ro, ωυ, and S were diagnosed at each grid point. The conditions for instability were then checked in the following order. If Ri < 1 the flow was deemed to be SI-unstable and the search was ended. If Ri > 1 but Ro < −1 the flow was deemed CI-unstable and the search was ended. If both Ri > 1 and Ro > −1 then the grid point was checked for whether . However, because AAI can occur under less stringent conditions than the other instability types (Wang et al. 2014), much of the domain was AAI-unstable at any given time, though often with negligibly small growth rates (e.g., McWilliams et al. 2004). Therefore, even if at a given point, it was deemed AAI-unstable only if S > 10−4 s−1, which we treated as a proxy condition for AAI growth rates similar to those of SI and CI (e.g., Dewar et al. 2015; Stamper and Taylor 2017). Each grid point was thus potentially unstable to more than one of these instability types, but the algorithm prevented any point from being flagged for more than one.
After each grid point was flagged for a particular instability type (or none), the grid was coarse-grained into boxes consisting of 10 × 10 grid cells (approximately 6 km × 6 km). Each coarse grid box was then flagged if it contained an “unstable point” of a particular type. This procedure was repeated for all 1440 snapshots in the output ensemble. Finally, for each coarse grid box we counted how many times it was flagged as unstable and divided this number by 1440, treating the result as an indicator of how often submesoscale instability occurred. Using only the condition AAI is the most common type of instability, with much of the domain unstable over 90% of the time (not shown). The additional condition S > 10−4 s−1 (i.e., restricting so that the AAI growth rates are fast enough to be “submesoscale relevant”) reduces the frequency substantially and essentially isolates the unstable locations inside the jet (Fig. 5c), similar to the other instability types.