Role of Mixed‐Layer Instabilities in the Seasonal Evolution of Eddy Kinetic Energy Spectra in a Global Submesoscale Permitting Simulation

A submesoscale‐permitting global ocean simulation is used to study the upper ocean turbulence in high kinetic energy (KE) regions. Submesoscale processes peak in winter so that the geostrophic KE spectra tend to be relatively shallow in winter ( ∼k−2 ) with steeper spectra in summer ( ∼k−3 ). This transition in KE spectral scaling has two phases. In the first phase (late autumn), KE spectra show the presence of two spectral regimes: ∼k−3 power‐law in mesoscales and ∼k−2 power‐law in submesoscales. The first phase appears with the onset of mixed‐layer instabilities, which convert available potential energy into KE, and this process results in a flattening of KE spectra at submesoscales. However, KE spectra at longer wavelengths follow ∼k−3 scaling associated with a forward enstrophy transfer. In the second phase (late winter), KE produced through mixed‐layer instabilities is transferred to larger scales, and k−2 power‐law also develops in mesoscales.

) and viscous dissipation (high wavenumbers, 2.5 max E k k   ) alter the spectra significantly, so are outside our scope. In panels (a-e) and (k-o) we show a dashed line for 2 E k  scaling and dotted line for 3 E k  scaling. In panels (f-j) vertical lines indicate the wavenumbers at which E Ro crosses 0.5 (solid lines are for winter and dashed lines are for summer). Color shading represents 95% confidence limits assuming a decorrelation time of 7 days (Rocha, Chereskin, et al., 2016) (p) The differences between the summer and winter spectral slopes (linear fits were computed in the wavenumber range [1/50, 1/12.5] cycles/km) for rotational KE spectra. The error bars represent 95% confidence limits computed by employing the Jackknife resampling method on 1,000 samples per region (Thomson & Emery, 2014 submesoscale KE and vorticity magnitudes peaking in the winter season (Dong et al., 2020b;Sasaki et al., 2017).
Submesoscale flows are generally defined as (1) E  Rossby number ( E Ro) flows and, in order to have a quantitative estimate of the submesoscale wavenumber range, we compute scale-dependent E Ro using the following relation (see Supporting Information S1 for the full derivation): and 0 E f is the Coriolis frequency at the center of the domain. In this work, we identify spatial scales that have 0.5 E Ro  as submesoscales (vertical lines in Figures 2f-2j indicate wavenumbers at which E Ro crosses 0.5). As seen in Figure 2, E Ro crosses 0.5 at wavenumbers [1/50, 1/25] cycles/ km in all regions and there is little seasonal variability in these estimates. Thus, it would be reasonable to categorize length scales smaller than about 50 km as submesoscales. Note that this submesoscale range is for our reference only because ( ) E Ro k estimates may change depending on how one defines E Ro (see e.g., Li & Lindborg, 2018). Our only purpose here is to have a quantitative reference to distinguish between mesoscales and submesoscales.
As seen in Figures 2f-2j, the divergent spectra magnitudes in the upper-ocean can be as large as the rotational ones within the submesoscale range (scales smaller than 50 km) in the summer season. Nevertheless, in the winter season, upper ocean submesoscale flows are predominantly rotational, which indicates that they are in near geostrophic balance. Note that rotational and geostrophic spectra are not necessarily equivalent due to the presence of weak vertical velocities associated with the balanced flow (Wang & Bühler, 2020). However, we expect the correction to be insignificant in the upper ocean because K k K k is mostly smaller than 1, especially in winter (Figures 2f-2j). Thus, this technical distinction does not affect the conclusions of the paper and we further discuss this aspect in the next section. In the following, we consider just the rotational KE spectra as we are interested in geostrophic turbulence.

KE Spectra in the Upper Ocean and Associated Spectral Scaling
In summer, 3 E k  power-law in rotational KE spectra at mesoscales and submesoscales is thought to be associated with the forward enstrophy transfer as in classical 2D and QG turbulence theories (Charney, 1971;Kraichnan, 1967). Khatri et al. (2018) computed enstrophy fluxes using satellite altimetry measurements and output from an eddy-permitting climate model simulation and showed that the enstrophy is indeed transferred downscale at mesoscales on the ocean surface.
On the other hand, 2 E k  scaling in winter KE spectra is generally found to be consistent with expectations from MLI, in which perturbations grow by extracting PE from lateral buoyancy gradients at submesoscales and release KE, which is transferred upscale (Boccaletti et al., 2007;Callies et al., 2016;Uchida et al., 2017). Many works have computed inter-scale KE transfers using outputs from high-resolution simulations and observed an upscale KE transfer in submesoscales (Ajayi et al., 2021;Dong et al., 2020b;Sasaki et al., 2017). As a result of this upscale transfer, KE spectra may flatten as in 2D turbulence (Kraichnan, 1967) and this effect is evident in rotational KE spectra in winter (Figures 2a-2e).
In winter, we expect KE production at both the interior baroclinic Rossby deformation scale ( R E L ) due to mesoscale baroclinic instability and at the mixed-layer deformation scale ( ML E L ) due to MLI. Lilly (1989) first proposed to study geostrophic turbulence in the presence of two forcing length scales, and this approach may be useful in understanding the nature of upper ocean turbulence. Assuming that there is sufficient scale separation between R E L and ML E L , both forward enstrophy flux (due to mesoscale baroclinic instability) and inverse KE flux (due to MLI) may control the slope in geostrophic KE spectra in the upper ocean at intermediate wavelengths (larger than ML E L and smaller than R E L ). Hence, the effective slope at these intermediate wavelengths could range between 5/3 E k  and 3 E k  (Lilly, 1989;Maltrud & Vallis, 1991). One caveat to this assumption is that 2D turbulence theory is not directly applicable in the presence of two forcing scales. A simultaneous presence of multiple instability processes adds further complexity, and it requires theoretical analysis to understand how cascades and the corresponding scaling laws would change in such scenarios.

of 13
Such analysis is outside the scope of this paper. Nevertheless, there is evidence that superposition of a forward enstrophy and inverse KE fluxes may occur in the upper ocean (see Figure 5 in Khatri et al., 2018).

Spectra in the Ocean Interior
Thus far we have focused on the upper ocean spectra at 21 m depth. For comparison, in Figures 2a-2e we also show the mean rotational KE spectra at 410 m depth. Little seasonal variability is seen at this interior depth. In Figure 2p, we show the differences between the summer and winter spectral slopes (computed in the wavenumber range [1/50, 1/12.5] cycles/km) in rotational KE spectra as a function of depth. The seasonal variability is pronounced in the upper 50 m, where the slope differences are close to 1, and can be associated with the seasonality in the mixed-layer depth, which varies between 10-20 m in summer and 100-150 m in winter (de Boyer Montégut et al., 2004).
For completeness, we also assess seasonally averaged PE spectra ( ( ) E P k in Figures 2k-2o, see Supporting Information S1 for computational details), which do not indicate any seasonal variability in spectral slopes. However, spectral slope magnitudes in PE spectra can differ between the ocean surface and ocean interior (Callies & Ferrari, 2013;Callies et al., 2016).

Interplay Between Mesoscale Baroclinic Instability and Mixed-Layer Instability
In this section we examine how rotational KE spectra evolve from summer to winter in the upper ocean. In Figure 3a, monthly averaged rotational KE spectra for July, November, and February are shown in the Kuroshio Current region (KCR). As discussed in Section 3, the KE spectra in July follow close to 3 E k  scaling, which can be associated with the downscale enstrophy transfer process as in 2D turbulence (Khatri et al., 2018;Kraichnan, 1967). In contrast, the scaling is close to 2 E k  in February due to enhanced submesoscale activity (e.g., see Figure 1 in Dong et al., 2020b). In the transition month of November, the KE spectrum shows both spectral scalings, with 3 E k  at wavelengths larger than about 40 km and 2 E k  at smaller wavelengths. To understand why these distinct phases appear in rotational KE spectra, we first analyze the time series of domain-averaged available PE (APE) and enstrophy ( Figure 3b). APE peaks in autumn (Oct-Nov) and this APE slowly decays by March-April accompanied by an increase in the mean enstrophy in late winter. The increase in APE in autumn months is associated with a quick decline in the 2 E N magnitude (Figure 3c, see Supporting Information S1 for time series plots in other regions) and deepening of the mixed-layer, which is generally caused by surface cooling in autumn (Iwamaru et al., 2010). On the other hand, buoyancy variance 2 E b    does not show any strengthening in winter. Instead, 2 E b    tends to be strongest in summer.
Moreover, as shown in Figure 3e, submesoscale APE to KE conversion ( E w b   , where E w and E b are the vertical velocity and buoyancy anomalies computed by removing spatial linear trends from each data snapshot following Uchida et al. (2017)) and vertical velocity magnitudes also increase during winter. This increase is consistent with the MLI mechanism, which converts APE into KE leading to an increase in submesoscale KE and enstrophy Sasaki et al., 2017). As seen in Figure 3e, the E w b    magnitude gradually increases during autumn months (Oct-Nov) and peaks in winter months (Jan-Mar). Thus, we can divide the MLI activity time into two periods: MLI onset period (Oct-Nov, green patch in Figure 3e) and MLI peak period (Jan-Mar, orange patch in Figure 3e). The spectral scaling in rotational spectra appear different in these two periods ( Figure 3a) and we interpret the differences in terms of MLI and 2D/QG turbulence theories.
Note that, in addition to MLIs, frontogenesis and turbulent mixing due to boundary layer parameterization schemes can also lead to an increase in the vertical buoyancy flux E w b    (Jing et al., 2020;McWilliams, 2016). However, Dong et al. (2020b) computed spectral KE flux in the Kuroshio extension region using LLC4320 output and found that enhanced E w b    and inverse KE transfer in submesoscales occur simultaneously during winter, and this is consistent with the expected seasonality in the MLI activity. It is worth noting that mixed layer deformation radii have seasonal dependence and the LLC4320 simulation can resolve MLIs primarily in winter months (Dong et al., 2020a). In summer, when mixed layers are shallow, MLI wavelengths are smaller, and instability growth rates tend to be weaker (Callies et al., 2015).

of 13
Moreover, instabilities are more likely to be damped out in shallow mixed layers. Thus, MLIs do not energize submesoscales in summer as strongly as in winter, and we expect a seasonal cycle in the MLI activity. This seasonality is observed in LLC4320 output. However, it is not clear how sensitive MLI wavelengths are to the model resolution. For example, wintertime MLI wavelengths may be smaller in a higher resolution model than the LLC4320 configuration. A comparison of different submesoscale-resolving simulations is required to investigate this aspect. is the complex conjugate), (e) domain-averaged APE to KE conversion rate E w b    at wavelengths smaller than 50 km (blue) and vertical velocity squared (red) time series. Panels (f-j) rotational KE spectra (k-o) Vortex energy spectra (p-t) bounds of geostrophic KE spectra, ( ) V E K k , in different regions at 21 m depth. 95% confidence limits are shown in color shading assuming a decorrelation time of 7 days in (a, f-o) and 2.5 days in (d) . The vertical gray lines in (f-t) are shown to indicate the change in spectral slope

Interpreting the Dual Power-Laws in Late Autumn
The coexistence of two power-laws in rotational KE spectra is seen during the MLI onset period (indicated with the green patch in Figure 3e). In Figure 3d, we show the monthly averaged wavenumber spectra of APE to KE conversion (   * ( ) E w b k   , plots for other months can be found in the Supporting Information S1).
In July,   * ( ) E w b k   at wavelengths smaller than 100 km ( 0.01 E k  cycles/km) fluctuates around zero (see Figure 3d), which supports our inference that APE to KE conversion takes place only at mesoscales through interior baroclinic instabilities. Also, a net zero APE to KE conversion at submesoscales in July indicates the presence of an enstrophy inertial range at 10-100 km wavelengths, since an inertial range can only be defined over a range of scales without any direct sources or sinks of energy (Vallis, 2017). Thus, it is reasonable to assume that 3 E k  scaling in July KE spectrum is associated with the enstrophy inertial range.
In November and February, positive values of   * ( ) E w b k   show that APE to KE conversion due to MLIs happens in the submesoscale range (at 0.02 E k  cycles/km in Figure 3d). Since MLIs energize submesoscales, we expect flattening of KE spectra at submesoscales, which is seen in the November KE spectrum in Figure 3a. At wavelengths larger than about 50 km, rotational KE spectra follow 3 E k  scaling, and this is mostly associated with the forward enstrophy transfer in mesoscale eddies as expected in 2D turbulence (Khatri et al., 2018;Kraichnan, 1967). The simultaneous presence of two power-laws in rotational KE spectra is evident in all regions (to help visualize the spectral slope break, we plot 3 ( ) Figures 3f-3j). Thus, we interpret the presence of two power-laws as the following: 3 E k  power-law is in accord with the forward enstrophy transfer due to eddies arising from mesoscale baroclinic instability, and 2 E k   scaling is due to KE injection into submesoscales through MLIs. Frontogenesis and SQG-type dynamics can also result in 2 E k  scaling at submesoscales in rotational KE spectra (Badin, 2013;Ragone & Badin, 2016). However, 2 E b    time series in Figure 3c suggests that frontal structures are relatively weak in fall and winter months than in summer, and fronts are less likely to be the cause of 2 E k  spectral scaling at submesoscales. Moreover, there is no inertial range at wavelengths corresponding to 2 E k  scaling because APE to KE conversion acts as a source for KE in the whole submesoscale range (Figure 3d).
As mentioned in the previous section, rotational KE spectra are not equal to geostrophic KE spectra since rotational spectra can have contributions from ageostrophic motions. Thus, to verify that the two power-laws indeed exist in the geostrophic KE spectra, we first employ the linear wave-vortex decomposition method (Bühler et al., 2014) to diagnose geostrophic (or vortex) energy spectra ( ) is the total KE spectrum and ( ) V E E k is the geostrophic component of the total energy (KE + PE) spectrum. As seen in Figures 3k-3o, the dual power-law signature is also present in vortex energy spectra. Since both KE and PE contribute to vortex energy spectra, the magnitudes of vortex energy spectra and rotational KE spectra differ significantly. However, both geostrophic PE and KE spectra follow the same power-laws in the enstrophy inertial range under QG approximation (Callies & Ferrari, 2013). Thus, a comparison between spectral exponents can still be made. However, we must note that the assumption of vertical homogeneity made in the linear wave-vortex decomposition method (Bühler et al., 2014) is violated in the upper ocean, which could make our diagnosed ( ) V E E k invalid. Therefore, we offer a more stringent analysis by examining the theoretical upper and lower limits of the geostrophic KE spectra (derivation can be found in the Supporting Information S1). In Figures 3p-3t, we present the lower bound, , for the geostrophic KE spectra. The gaps between the bounds are too wide in the GSR and KCR regions for us to make any interpretations about spectral slopes; however, in the other three regions (Figures 3r-3t), the bounded range still displays a coexistence of two power-laws.
In this direction, Choi et al. (2019) also reported the coexistence of two spectral regimes in rotational KE spectra (transition from 3 E k  to 5/3 E k  power-law at about 50 km) computed from geostationary ocean color satellite measurements (for April month, when MLIs are relatively weak) in the East/Japan Sea region. This region sees a persistent presence of submesoscale filaments and strong buoyancy gradients. Hence, Choi et al. (2019) attributed the change in the spectral exponent to surface-QG dynamics at length scales smaller than 50 km since surface-QG theory predicts 5/3 E k  power-law for KE spectrum in the buoyancy cascading 9 of 13 inertial range (Pierrehumbert et al., 1994). However, in the regions considered in this work, we did not find sufficient evidence that supports surface-QG behavior, and the coexistence of two power-laws in geostrophic KE spectra can be understood in terms of seasonal variability in the MLI activity. Even so, we cannot rule out that surface-QG signature may be visible in other parts of the global ocean. Although the MLI mechanism can be studied in a layered SQG system, MLI and SQG dynamics differ considerably (Blumen & Wu, 1983;Callies et al., 2016). In SQG dynamics, a mesoscale strain field creates strong submesoscale fronts on the ocean surface. On the other hand, MLIs convert APE stored in mesoscale buoyancy gradients into KE and energize submesoscale flows in the mixed layer. Hence, these two processes affect surface KE spectra in different ways while possessing similarities in the flow structure.

KE Spectra in Late Winter
It is expected that APE to KE conversion through MLIs at submesoscales results in an upscale KE transfer (Ajayi et al., 2021;Callies et al., 2016;Dong et al., 2020b) and an inverse KE transfer tends to flatten the KE spectrum. Since   * ( ) E w b k   peaks in the wavenumber range [1/50−1/10] cycles/km (Figure 3d), MLI wavelengths lie between 10 E 50 km, which agrees with the estimates provided in Dong et al. (2020a). Thus, an upscale KE transfer from MLI wavelengths is required for the appearance of the relatively shallow 2 E k  scaling at wavelengths longer than 50 km in geostrophic KE spectra. This suggests that 2 E k  spectral scaling in mesoscales in the February KE spectrum is partly due to the upscale KE transfer from MLIs (also see discussions in Dong et al., 2020b;Sasaki et al., 2017). Note that the process of energizing mesoscales through the inverse KE transfer can take a few weeks to a couple of months (Qiu et al., 2014). We hypothesize that this delay is the main reason why power-law scalings in rotational KE spectra appear different in late autumn and late winter (Figure 3a), even though MLIs are active during the whole period (see Figure 3e).
In general, we expect the spectral exponent in the geostrophic KE spectrum throughout mesoscales and submesoscales to be controlled by an inverse KE flux (due to MLI) and a forward enstrophy flux (due to interior baroclinic instability). The temporal evolution of KE spectral slope, APE, and enstrophy in Figures 3a-3c confirms the roles of these spectral fluxes in different months. In summer, spectra show 3 E k  scaling associated with the forward enstrophy flux as in 2D/QG turbulence (Charney, 1971;Kraichnan, 1967). In late winter, MLIs occur at submesoscales (Figure 3d), and both inverse KE and forward enstrophy fluxes can contribute to shaping the KE spectrum. Hence, we do not expect an inertial range at 10 -100 E km wavelengths in winter, and a theoretical power-law scaling for the KE spectrum is not possible. The appearance of 2 E k  power-law in late winter KE spectra could be a mere coincidence. Nevertheless, the flattening of the KE spectrum in winter is expected due to the inverse KE transfer from MLIs and is robust across the oceanic regions we analyzed. We provide a schematic in Figure 4 explaining the different phases of geostrophic KE spectrum evolution in the upper ocean.

Discussion and Conclusions
In this study, we used output from a 1 / 48 E  global ocean simulation (LLC4320) to examine the behavior of upper ocean geostrophic turbulence in different high eddy activity regions in Western Boundary Currents and the Southern Ocean. In agreement with previous studies, we found a strong seasonality in submesoscale turbulence and the associated geostrophic kinetic energy (KE) spectra (Callies & Ferrari, 2013;Sasaki et al., 2017). Specifically, rotational KE spectra tend to follow 3 E k   power-law in summer consistent with two-dimensional (2D) and quasi-geostrophic (QG) turbulence (Charney, 1971;Kraichnan, 1967), and 2 E k   power-law in winter. It is proposed that mesoscale and submesoscale turbulence in the upper ocean geostrophic flows can be understood in terms of QG turbulence and mixed-layer instability (MLI) mechanism. The spectral exponent in the KE spectrum is affected by the downscale enstrophy transfer due to interior baroclinic instability at mesoscales as well as the upscale KE transfer due to MLIs at submesoscales.
We described two distinct physical phases in the seasonal transition of upper ocean rotational KE spectra from 3 E k  scaling in summer to 2 E k  scaling in winter. The first phase occurs in late autumn, during which MLIs start to develop and convert available potential energy into KE at submesoscales leading to a flattening of KE spectra at wavelengths smaller than about 50 km. Rotational KE spectra at relatively larger length 10 of 13 scales still follow 3 E k  scaling associated with the forward enstrophy transfer due to mesoscale baroclinic instabilities. Hence, rotational KE spectra show a simultaneous presence of two power-laws in late autumn: 3 E k  at wavelengths larger than about 50 km and 2 E k  at smaller length scales. In the second phase (late winter), MLIs are at their peak and KE generated through MLIs at submesoscales is transferred to mesoscales through an upscale (inverse) KE transfer process. As a consequence of this upscale KE transfer, 2 E k  spectral scaling develops in rotational KE spectra at wavelengths of 10-100 km (see also Dong et al., 2020b;Sasaki et al., 2017). The KE transfer from submesoscales to mesoscales may take a few months after the onset of MLI (Qiu et al., 2014), and we hypothesize that this time lag is why we measure distinct phases in the temporal evolution of rotational KE spectra. We further compared rotational KE spectra against geostrophic total energy spectra computed using the linear wave-vortex decomposition method (Bühler et al., 2014) and reached the same conclusions. We present a schematic in Figure 4 that summarizes the physical mechanisms involved with these seasonal transitions. Note that it is required to compute both KE and enstrophy fluxes to confirm the superposition of inverse KE and downscale enstrophy fluxes at submesoscales. Moreover, geostrophic KE spectra computations can be further improved by considering the k  scaling associated with the downscale enstrophy cascade in wavelengths smaller than baroclinic instability wavelength (BI), and an inverse KE cascade is expected at longer wavelengths. (b) In late autumn, mixed-layer baroclinic instabilities (MLI indicates the corresponding wavelengths) extract energy from lateral buoyancy gradients and energize submesoscales and flatten the KE spectrum at submesoscales. (c) KE generated through MLIs is transferred upscale to mesoscales. In this scenario, both upscale KE flux and downscale enstrophy flux affect the spectral slope, and a shallower   k 2 scaling develops in the mesoscale KE spectrum. Black arrows denote the directions of the spectral fluxes of KE and enstrophy (note that the dimensions of these fluxes are not the same). Smallscale dissipation mainly occurs at the model grid scale and, at large-scales, several factors contribute to KE arrest, for example, bottom drag, Rossby waves.