The ocean’s mechanical energy budget can be used to quantify the magnitude of reservoirs of energy (kinetic and potential) and the rate of conversion between these reservoirs. These conversions include the generation of eddy kinetic energy via the action of baroclinic and barotropic instability. In this study, we estimate these conversion rates using a thickness-weighted energetics framework along the lines of Bleck (1985); Aiki et al. (2016) and Barthel et al. (2017). The thickness-weighted framework is formulated in density coordinates rather than the more common approach in Eulerian coordinates (e.g. Lorenz, 1955; Chen et al., 2014). This approach is justified by upwelling in the Southern Ocean being an approximately adiabatic process – diapycnal mixing in the ACC has a much lower magnitude than isopycnal mixing (Naveira Garabato et al., 2016). The advantage of the thickness-weighted energetics framework is that it explicitly includes the energy conversion related to interfacial form stress, associated with baroclinic instability (Aiki et al., 2016; Barthel et al., 2017).

Barthel et al. (2017) expanded the two-layer energetics framework of Aiki et al. (2016), based on the thickness-weighted energetics theory of Bleck (1985), to include a free surface and a friction parameter. We extend the two-layer framework to three layers, which can be directly extended to n levels for use with global ocean model output. The three layers of our energy framework, with densities ρ₁< ρ₂ < ρ₃, are each associated with a kinetic energy reservoir,

with h_i and u_i = (u_i , v_i ) the thickness and horizontal velocity respectively of each layer i (Figure 1). Vertical velocities are neglected via a standard scale analysis, since horizontal scales in the ocean are much larger than the vertical (Cushman-Roisin and Beckers, 2011). ρ₀ is a reference density, 1035 kg m^-3. The available potential energy (APE) can be defined for each interface from a reference height (here, defined as the time-mean interface depth, though the choice does not affect the final result):

Here, the APE₀ is the available potential energy of the ocean’s free surface, defined by deviation η₀, and APE_j (j = 1, 2) is the available potential energy associated with each of the two internal interfaces between density layers of the ocean, with interface deviation η_j (Figure 1). The densities of the three layers are used to define reduced gravities g′₁= (ρ₂ – ρ₁)g/ρ₀ and g′₂ = (ρ₃ – ρ₂ ) g/ρ₀.

Following Aiki et al. (2016) and Barthel et al. (2017), the Navier Stokes equation for a rotating, Boussinesq fluid is used to describe the time dependence of the kinetic energy reservoirs, and is simplified by the layer-thickness continuity equation, under the assumption that the fluid is adiabatic. Time dependence of the APE reservoirs can also be calculated by differentiating η_j with respect to time, and simplifying with the Montgomery potential anomaly Φ_i (Cushman-Roisin and Beckers, 2011),

An equation for the advection of the Montgomery potential flux is also formulated (see Aiki et al. (2016) for details).

Having defined energy reservoirs and their time-dependence, Reynolds averaging is introduced to decompose quantities into mean and eddy quantities. The layer thicknesses and interface heights are decomposed as h i ( x , t ) = h i ¯ ( x ) + h i ′ ( x , t ) , with overlined quantities the time-mean and primed quantities the deviation such that h i ′ ¯ = 0 .

The time-mean available potential energy reservoirs are thus split into mean (MAPE) and eddy (EAPE) components:

For the horizontal velocities, we could use the same (Ᾱ, A′) Reynolds decomposition. However, following Bleck (1985), we instead choose to decompose the layer horizontal velocities into a thickness weighted mean velocity and deviation (Â, A"), defined as u_i = û_i + u_i ", where the thickness weighted mean is expressed as

The use of this thickness-weighted mean is a key difference between the thickness-weighted energetics framework and that of Lorenz (1955). The kinetic energy is thus decomposed into mean (MKE) and eddy (EKE) components:

As in Aiki et al. (2016), the equations describing the time-dependence of these reservoirs are low-pass filtered in time, yielding equations describing the conversions into and out of the mean energy reservoirs. These mean equations are subtracted from the time-average of the original time-dependent equations, yielding equations describing conversions from the eddy energy reservoirs. The equation for the eddy kinetic energy in each layer is

Here, F_τ,i parameterises frictional forces (both interior viscosity and bottom drag) in each layer. The three-layer eddy kinetic energy budget, Eqn. (6), is identical to the two-layer framework of Barthel et al. (2017), but there are minor changes in the APE equations due to the addition of a third APE reservoir. Some terms appear in the equations for multiple energy reservoirs with opposite sign, indicating a transfer of energy between those reservoirs. An energy diagram can thus be constructed, in which conversions of energy act to link reservoirs (Figure 2). The mean and eddy Montgomery potential flux equations are also used to link reservoirs, where terms are not directly found in the equations for two different reservoirs. Two conversions linking mean and eddy reservoirs are identified:

These two conversion terms, the interfacial form stress conversion and Reynolds stress conversion, can be taken as indicators of the action of baroclinic and barotopic instabilities respectively. We calculate these terms, as well as the eddy kinetic energy generation terms in Eqn. (6), using output from a high-resolution ocean model to investigate the spatial distributions of instability processes involved in Southern Ocean upwelling.

The model used in this study is a high-resolution global ocean and sea-ice model, ACCESS-OM2-01 (Kiss et al., 2020), which has 75 vertical levels and 0.1° spatial resolution. The model is hence eddy-resolving at Southern Ocean latitudes (Hallberg, 2013). ACCESS-OM2-01 uses ocean numerics and configurations based on the Modular Ocean Model version 5 (MOM5) (Griffies, 2012). The model is forced at the ocean surface by a repeat year forcing (May 1990-April 1991, a neutral year with regard to interannual climate modes (Stewart et al., 2020) of the JRA55-do atmospheric reanalysis (Tsujino et al., 2018). Using a repeat year forcing ensures the mean state of the ocean is relatively constant, for ease of calculating eddy-related quantities, and to ensure that interannual variability is due to oceanic processes. The model has been spun up for 270 years to decrease initial model drift, and the following 10 years of model output are analysed.

ACCESS-OM2-01 has an overturning circulation within the range of observation-derived estimates (Figure 3). The circulation can be analysed using the meridional overturning streamfunction, defined in density-latitude coordinates as

where the ≡ symbol describes the equivalence of the depth coordinate (left) and density coordinate (right) calculations. This expression is the cumulative integral from the bottom of the ocean (σ_B) of the total northward transport in density layers v h ¯ ( x , y , σ ) , integrated over all longitudes and time-averaged. Figure 3A shows the meridional overturning streamfunction calculated in σ₁ (potential density referenced to 1000 dbar) coordinates, mapped to the zonal-mean isopycnal depths. Flows can be interpreted to follow streamfunction contours (arrows). Note that this mapping is imperfect, because isopycnal depths vary with longitude so water exists at depths greater than the zonal mean, and is only provided for visualisation purposes.

Figure 3B shows the modelled volume transport at 32 °S in density classes defined by Talley et al. (2003), compared with three observation-derived inverse model estimates (Sloyan and Rintoul, 2001; Talley et al., 2003; Talley, 2008) between 30-40 °S. Additional resolution is provided for the Circumpolar Deep Water (southward moving, upwelling water). The modelled overturning circulation generally lies within the range of these sparse observational studies, but the model has too little Antarctic Bottom Water transport at 32°S compared with observations. Actual production of bottom water near Antarctica in ACCESS-OM2-01 compares well with observations (Moorman et al., 2020), an uncommon feat for ocean models, which often do not represent the process of bottom water formation near Antarctica accurately (Heuzé, 2021). The low bottom water transport at 32 °S is hence likely to be due to enhanced mixing in the model, owing to poor resolution of downslope flows and spurious numerical mixing. The modelled bottom water is thus over-mixed with lighter densities as it travels northwards, decreasing transport in the densest class.

We have only compared the model to observations at 32°S, where there is a range of observational studies, hence we cannot comment on the accuracy of the modelled overturning at the more southern Antarctic Circumpolar Current latitudes. However, since the modelled upwelling lies within the large range of observational estimates of the overturning circulation at 32°S, the modelled overturning circulation is likely to be sufficiently realistic for the purpose of understanding the physics involved in its upwelling arm. Additionally, Kiss et al. (2020) show good agreement of other diagnostics of ACCESS-OM2-01, such as temperature and salinity transects, Antarctic Circumpolar Current transport and sea-ice coverage, with observational datasets. Solodoch et al. (2022) provide further validation of the model in the Southern Ocean. The eddy kinetic energy distribution of ACCESS-OM2-01 is also similar in magnitude and spatial distribution to satellite observations (see Supplementary Figure S1).

The thickness-weighted energetics framework requires calculations in density coordinates, but ACCESS-OM2-01 output uses depth coordinates, and thus requires a vertical coordinate transformation. We choose 31 σ₁ density levels so that layers are 100-200 m thick (Figure 4A). We choose 1000 dbar (~1000 m) as the reference pressure for potential density, because the upwelling density range is approximately centred around 1000 m (Figure 3A). We recognise that neutral surfaces would be more accurate than potential density (McDougall, 1987; Jackett and McDougall, 1997), but they are computationally expensive.

The ACCESS-OM2-01 model output is transformed into the 31 density levels following the binning method of Lee et al. (2007). The binning process is performed using daily averaged output. The high temporal resolution is important because correlations in fluctuations in layer thickness and velocity are used in energy and transport calculations. Sub-daily timescales are unlikely to be relevant for mesoscale dynamics in the Southern Ocean.

The thickness-weighted energy conversion framework requires a Reynolds decomposition into time-mean and eddy quantities, as does the calculation of eddy transport as does the calculation of eddy transport. Eddy quantities are calculated as the deviations from a 10 year time-mean. Aiki and Richards (2008) find calculations of energy conversions to be robust to the choice of time-averaging. Seasonal variability in eddy kinetic energy is predominantly in the upper 100 m of the ocean (Rieck et al., 2015), far above upwelling processes, so we do not account for seasonality in the time-mean. Eddies in the Southern Ocean have lifetimes of around 6 months (Chelton et al., 2011), suggesting a 10 year time period would be sufficient to calculate accurate time-means that are not biased by a single eddy. This definition means the fluctuation term includes all variations, including seasonality, interannual variability and waves as well as turbulence.

We define upwelling to be the eddy component of the volume transport across time-mean sea surface height (SSH) contours, which are indicative of Antarctic Circumpolar Current streamlines. The mean component of the volume transport in the interior density layers does not contribute to the net southward flow (shown later in the Results section), and therefore is not included in the analysis. Computing the eddy transport across SSH contours, rather than across lines of constant latitude, removes the majority of the ACC flow from the calculated upwelling transport. Removing the mean transport also minimises the contribution of ACC standing meanders to the upwelling transport across SSH contours. The mean ACC flow at depth does not exactly follow SSH contours,but we find that other choices of contours, such as Montgomery potential or isopycnal depth contours, yield similar results (not shown). The cross-contour volume transport across a contour C of scalar field S through an area dA tangent to the contour is

For a density layer in the ocean with layer thickness h, and nearly horizontal flows, this reduces to

where we have taken dA to extend throughout the layer, and the vector n s = ∇ S / | ∇ S | .

The quantity t(ρ,x , C ) has units of area per unit time, and is properly a linear density of volume transport, but for simplicity we refer to it also as the cross-contour volume transport. When integrated over all points along a contour C , we obtain the circumpolar transport T(ρ, C ). The three-dimensional cross-contour volume transport t can be decomposed via a Reynolds decomposition into mean and eddy components, of which the eddy component, t_eddy, is our measure of upwelling

The eddy cross-contour volume transport is calculated for 29 evenly spaced SSH contours between -0.1 and -1.5 m (Figure 4B), which cover the latitudinal extent of the Antarctic Circumpolar Current. The cross-contour transport is computed across the faces of the discrete model tracer grid cells to guarantee exact closed circumpolar contours. The eddy transport is smoothed with a 10 along-contour grid point running average to reduce noise (Figure 5A). This can be thought of as reducing the influence of the rotational component of eddy transport, which obscures the physically relevant divergent component (Lau and Wallace, 1979). The decomposition is imperfect (Griesel et al., 2009) as only a circumpolar average can truly remove all of the rotational eddy transport, and thus we do not claim to have exactly calculated the divergent eddy transport (note that no such divergent-rotational decomposition is perfect, as solutions are not unique in the presence of boundary conditions (Fox-Kemper et al., 2003). However, we expect that including part of the rotational eddy transport will only increase noise, and that hence the significance of any results would likely increase if the rotational component had been entirely removed. The choice of 10 grid cell smoothing is not unique and is chosen to be of a similar length scale to observed eddy sizes (Klocker and Abernathey, 2014), but our results are not sensitive to this choice.

Note that our definition of eddy transport (Eqn. 10) differs from the more traditional definitions of overturning circulation eddy-mean decompositions (e.g. Wolfe and Cessi, 2009). In the traditional method, the eddy component is taken as the difference between the residual transport in density coordinates and the time-mean Eulerian transport in depth coordinates mapped to density coordinates. The time-mean Eulerian transport calculation therefore includes velocity contributions from layers outside the isopycnal of interest. In contrast, our isopycnal mean method uses the time mean velocity u ¯ within a single density layer, with velocities averaged across a range of depths following vertical movement of the isopycnal. The streamfunction associated with mean transports calculated via the isopycnal mean method will not be closed due to temporal variation in the range of isopycnals present, and hence the isopycnal mean method is not ideal for streamfunction calculations, but we argue that this method is more physically justified for the purposes of analysing transport in the upwelling density range of the ocean. We present the isopycnal mean method in Figure 5A and the traditional method in Figure 5B, the latter of which contains more noise and less prominent upwelling hotspots, even after the same smoothing filter was applied. We thus use only density coordinates in all calculations.

We analyse upwelling dynamics over a density range that is representative of the zonally averaged upwelling (bounded by the orange contours in Figure 4A). We refer to this density range as the ‘upwelling arm’. The density range is chosen to select densities of water that, in the circumpolar average, have net eddy transport orientated southward and do not intersect the surface or ocean bottom during the 10 year simulation period. These criteria ensure that northward flowing bottom water and surface layer water are not included, as different dynamics are present in these boundary regions of the ocean. Specifically, flows are more likely to include diabatic components at boundaries, which contradicts the adiabatic assumption underpinning the thickness-weighted energetics framework. Outcropping density layers also need not be considered. Since isopycnals slope in the Southern Ocean, the density range satisfying these criteria varies with latitude and SSH contour. We thus choose a density range such that all contours with SSH less than or equal to -0.6 m satisfy the two criteria, resulting in the upwelling arm density range 1032.2 < σ₁ < 1032.5 kg m^-3.

The meridional overturning streamfunction in σ₁-latitude coordinates is presented in Figure 6A. Streamfunction contours are approximately horizontal in density-latitude space except near the surface (lightest colours), indicating that flow occurs nearly a diabatically along isopycnals. Thus, assumptions that upwelling occurs along isopycnals are valid except near the surface. The decomposition of the residual overturning circulation into mean and eddy components is shown in Figures 6B, C respectively, using the same Reynolds decomposition as for eddy transport which computes the mean overturning from the mean velocity on isopycnals, rather than the usual depth-space Eulerian mean calculation (see Methods). Southward eddy transport (Figure 6C) is intensified in the upwelling arm (between the horizontal orange lines) and between 500-2000 m depth (dotted and dashed grey lines).

The analysis of zonally integrated eddy transport is useful for quantifying the total upwelling. However, upwelling is known to vary with longitude (Tamsitt et al., 2017). We first consider the transport along a single contour – the -1.2 m SSH contour, highlighted in Figure 4B – to illustrate the spatial distributions of eddy and mean transport contributions. The horizontal spatial distribution of eddy upwelling across the -1.2 m SSH contour (Figure 7A) shows that southward eddy transport is non-uniform, and enhanced in regions that are colocated with or are downstream of topographic features. The along-contour cumulative integrals of eddy, mean and total transport across the -1.2 m SSH contour are shown in Figure 7B. There are five regions of high southward eddy transport (black) that occur near topography, defined by large changes in the cumulative integral. We refer to these regions as topographic hotspots of upwelling. The residual (blue) transport across the -1.2 m SSH contour is also southward, consistent with the zonally averaged the residual overturning streamfunction (Figure 6), whilst the mean (red) transport is noisy and small in the along-contour integral. The small-scale variation in mean and total transport, due to the mean Antarctic Circumpolar Current crossing the SSH contours, is not present in the eddy transport. The smooth eddy transport suggests that the Reynolds decomposition method isolates the eddy transport signal from the mean Antarctic Circumpolar Current, providing confidence that the cross-contour eddy transport method accurately estimates upwelling. There is a total of 19 Sv of southward eddy transport across the -1.2 m SSH contour, integrated vertically over the upwelling arm (Figure 7B). The eddy transport in the hotspot regions (grey bars) dominates the transport; the eddy transport in these regions contributes 76% to the total eddy transport, but occurs in only 28% of the circumpolar longitude range for this contour. A similar dominance of hotspot regions occurs across the majority of Southern Ocean SSH contours.

Upwelling is spatially variable, so we analyse the spatial distribution of upwelling across SSH contours, densities and depths at each hotspot. Figure 8A shows the eddy transport integrated along circumpolar contours for each density bin, T_eddy(ρ, C ). The SSH axis ( C ) is a proxy for latitude. The upper black line is the along-contour average surface density, and the lower line indicates where T_eddy > 0 (northward-flowing bottom water not associated with upwelling). Integrating between these black lines gives upwelling as a function of SSH contour, expressed as ∑ ρ_top <ρ<ρ_bottom T(ρ,  C _hotspot). Figure 8C shows this distribution of upwelling across SSH contours at each hotspot, C _hotspot. Hotspots are defined by the following longitude ranges: Macquarie Ridge 130 °E-155 °W, Pacific-Antarctic Ridge 155 °W-100 °W, Drake Passage 80 °W-5 °E, Southwest Indian Ridge 5 °E-55 °E, and Kerguelen Plateau 55 °E-130 °E. Upwelling occurs over the whole range of sea surface contours, indicating upwelling occurs at all points across the Antarctic Circumpolar Current. However, the SSH contours where upwelling is maximised depend on the hotspot, with the Southwest Indian Ridge, Kerguelen Plateau and Pacific Antarctic Ridge having more upwelling over southern contours and Drake Passage and Macquarie Ridge in central contours. Furthermore, Drake Passage is the hotspot with the most upwelling transport (Figure 8C), which is attributed to its substantial topography and longitudinal extent (Figure 7A), while the other hotspots have similar contributions to the circumpolar upwelling transport.

Figure 8B shows the eddy transport summed along contours at each hotspot T(ρ, C _hotspot) between the black lines of Figure 8A, averaged over all the SSH contours C _hotspot and normalised by the average layer thickness h, and hence shows upwelling in density space. The density axis of Figure 8B is then mapped back into depth space (Figure 8D). The density ranges of eddy upwelling hotspots (Figure 8B), are similar with small variations. Kerguelen Plateau in particular has its upwelling skewed to lighter densities and shallower depths (Figure 8D) than the other hotspots. The topography of Kerguelen Plateau is shallower than other hotspots, with large parts of the plateau shallower than 1000 m (Figure 7A). The interaction of the flow with topography at different depths is thus likely to be relevant to the density of upwelling water. Figure 8 shows that upwelling occurs at densities outside of the upwelling arm density range. However, the eddy transport in the upwelling arm density range accounts for ~75% of the total southward eddy transport, indicating the upwelling density arm is a sufficient representation of upwelling water.

As motivated in the Energy Conversion Framework section, we use the interfacial form stress conversion calculated from a layer-thickness energetics framework to indicate the spatial distribution of baroclinic instability in the high-resolution model. The interfacial form stress conversion (blue arrows in Figure 2) was calculated in the 31 levels of the transformed ACCESS-OM2-01 model output. From 35-70°S, the interfacial form stress conversion results in 0.29 TW converted from mean to eddy energy, consistent with the action of baroclinic instability. The spatial distribution of eddy transport is compared with contours of the interfacial form stress conversion, for the five topographic hotspots of upwelling (Figure 9), and shows that the energy conversion is enhanced in the same locations as upwelling. Since eddy transport is calculated over SSH contours, the eddy transport distributions have been interpolated between SSH contours onto the model grid to better compare with the interfacial form stress conversion, which is computed in every horizontal grid cell. Grid cells lying between SSH contours are assigned a transport weighted by the distance between the two closest points on adjacent contours. The resultant spatial distribution of eddy transport is not uniquely defined at grid points between contours since other interpolation methods exist, but provides a visual guide. Southward eddy transport regions (blue shading) have a striking spatial correspondence with positive interfacial form stress conversion values (purple contours), and conversely northward eddy transport (red colours) tends to be located where the interfacial form stress conversion is negative (orange/yellow contours). The spatial correspondence is best highlighted at the Southwest Indian Ridge hotspot (Figure 9A). There are, however, regions where the correspondence is less strong (e.g. the north-western part of the Pacific-Antarctic Ridge, Figure 9D), which may be explained by shortcomings in the assumption of adiabatic and geostrophic mean flows (see Discussion).

To better quantify the visible spatial correspondence between eddy transport and the interfacial form stress conversion, we compare the distributions of southward eddy transport and the interfacial form stress conversion along the -1.2 m SSH contour (Figure 10A). A lag-correlation analysis (Figure 10B) shows the Pearson regression correlation coefficient for each spatial lag between the two fields, where a positive lag distance indicates that the transport is upstream of the energy term. Regions of significance greater than 95% confidence are highlighted, taking into account internal correlations in the spatial distributions (Santer et al., 2000). There is a significant correlation (Pearson regresssion correlation coefficient r = 0.80) between the interfacial form stress conversion and eddy transport with zero lag (Figure 10B) for the -1.2 m SSH contour. Similar results (0.63 ≤ r ≤ 0.80) hold for all other contours with SSH less than -0.6 m. In comparison, the Reynolds stress conversion associated with barotropic instability (orange arrows in energy budget, Figure 2) is enhanced in the same regions as eddy transport, but is not correlated (Figure 10C), with maximum correlation r = 0.26 for the -1.2 m SSH contour (Figure 10D) and r ≤ 0.53 for other contours. The high spatial correspondence between the interfacial form stress conversion and eddy transport indicates that the two fields are dynamically linked, and thus that upwelling is driven by baroclinic instability. This result is consistent with the idealised study of Barthel et al. (2022). Quantifying this baroclinic conversion could thus be an avenue for improving parameterisations of eddy transport in coarse resolution models.

Past idealised and observational studies have suggested that, like upwelling, eddy kinetic energy is enhanced downstream of topography (Chelton et al., 1990), but that there is a spatial offset between eddy kinetic energy and upwelling (Barthel et al., 2022), eddy heat fluxes (Abernathey and Cessi, 2014; Foppert et al., 2017), and baroclinic activity (Bischoff and Thompson, 2014; Chapman et al., 2015). We investigate this phenomenon using ACCESS-OM2-01. The modelled eddy kinetic energy distribution (Figure 4B) is spatially variable, and enhanced in similar regions to the upwelling. Figure 10E compares the along-contour distribution of eddy transport and eddy kinetic energy for the -1.2 m SSH contour. There is a significant correlation (r = 0.71; Figure 10F), maximised at a lag distance of 126 km. When this peak distance is averaged over all contours with sea surface heights less than -0.6 m, the eddy kinetic energy is found to lag the eddy transport by 110 ± 20 km (where the uncertainty is the standard error in the mean), and has maximum lag correlations of 0.50 ≤ r ≤ 0.71. This distance is in accordance with an observational study of eddy heat fluxes (Foppert et al., 2017). For parameterisations of eddies that use eddy kinetic energy, this is an important result, as the downstream distance only covers two grid cells in a coarse 1° model, and is thus unlikely to negatively affect the accuracy of parameterised eddy transport. For eddy parameterisations used in models with a resolution finer than 1°, the spatial offset between EKE and eddy transport may need to be considered.

The same analysis can be performed at each hotspot for contours where there is a significant correlation between eddy kinetic energy and upwelling, yielding downstream distances of 50 ± 40 km at Macquarie Ridge, 70 ± 10 km at the Pacific-Antarctic Ridge, 100 ± 20 km at Drake Passage, 230 ± 60 km at the Southwest Indian Ridge and 80 ± 20 km at Kerguelen Plateau. Individual contours at each hotspot deviate from these mean distances, but eddy kinetic energy is persistently downstream of eddy transport.

We have shown that in ACCESS-OM2-01, eddy kinetic energy maxima are found downstream of upwelling hotspots. Since upwelling and the interfacial form stress conversion are highly correlated, the eddy kinetic energy is also found downstream of the baroclinic energy conversion process. There are several candidate mechanisms for this downstream offset that have been suggested in the literature, including the advection of eddies (and therefore eddy kinetic energy) (Thompson and Naveira Garabato, 2014), the advection of pressure perturbations (Montgomery potential flux in a layered framework) (Chang and Orlanski, 1993; Abernathey and Cessi, 2014; Chapman et al., 2015), the barotropic generation acting downstream of the baroclinic generation (Barthel et al., 2022), a time-delay between the instability process and eddy kinetic energy generation (Chen et al., 2016) (which via advection of eddy potential energy due to barotropic eddies results in a spatial offset), and bathymetric configuration and eddy-mean flow interactions (Foppert et al., 2017). Using the thickness-weighted energetics framework, the plausibility of two mechanisms, relating to the advection of eddy kinetic energy and the action of barotropic instability, can be investigated.

Eddy kinetic energy advection terms and the Reynolds stress conversion associated with barotropic instability can be identified in the eddy kinetic energy evolution equation (6). When integrated over the full depth of the ocean and between sea surface heights -0.6 m and -1.5 m, the advective terms contribute 0.51 GW to energy conversions, and the Reynolds stress conversion 2.8 GW. These conversions are small compared with the baroclinic eddy kinetic energy generation term, u i ″ · h i ∇ ϕ i ″ ¯ , which when integrated over the same region is three orders of magnitude larger at 1.3 TW. The small magnitudes of the advective and Reynolds stress conversion terms in this region are only partly due to cancellation of local positive and negative conversions – the local absolute value of the energy conversions integrated over the same SSH range is 0.18 TW for the advective terms, 0.17 TW for the Reynolds stress conversion and much larger at 10 TW for the baroclinic eddy kinetic energy generation term. Furthermore, there is no significant (to 95% confidence) correlation between the eddy kinetic energy and Reynolds stress conversion along SSH contours (Figure 10D). Correlation analyses hence do not support a mechanism whereby the barotropic instability facilitated eddy kinetic energy generation downstream of eddy transport. Our analysis thus suggests that neither the downstream advection of eddy kinetic energy, nor the downstream action of barotropic instability from baroclinic instability, are likely to be primary mechanisms for creating the spatial offset between upwelling and eddy kinetic energy. Our findings show a reduced magnitude and thus importance of the barotropic energy conversion compared to the idealised two-layer model analysis of Barthel et al. (2022) – the authors acknowledge their simplified set-up may overstimulate barotropic instability (Youngs et al., 2017; Barthel et al., 2017; Barthel et al., 2022). This difference suggests that the many-layered, global ocean model employed here experiences different dynamics to the idealised model, and thus that it is important to verify the findings of idealised models in more realistic global models.

We have not verified other proposed mechanisms for downstream eddy kinetic energy, but suggest that future work could use a high-resolution global ocean model to investigate the downstream advection of pressure perturbations using wave activity flux theory (Takaya and Nakamura, 2001; Chapman et al., 2015). The wave activity flux framework includes an ageostrophic Montgomery potential flux term that has been suggested in idealised studies (Chang and Orlanski, 1993; Chapman et al., 2015) and observations (Foppert, 2019) to dominate horizontal wave activity transport. Via this mechanism, pressure perturbations are quickly swept downstream where they convert eddy potential energy into eddy kinetic energy (Abernathey and Cessi, 2014). Our time-mean energy conversion analysis is also unsuitable for examining timescales between baroclinic instability and eddy kinetic energy generation, along the lines of Chen et al. (2016), and is a subject for future work.