Nonlinearly interacting entrainment due to shear and convection in the surface ocean

Large-eddy simulations were performed to investigate the entrainment buoyancy flux at the mixed layer base due to nonlinearly interacting shear-driven turbulence (ST) and convective turbulence (CT). The fluxes due to pure ST and pure CT were first evaluated, and their scalings were derived. The entrainment flux due to coexisting ST and CT was then evaluated and compared to the scaling-based fluxes due to the pure turbulences. It was found that nonlinear effects reduce the entrainment flux by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$30 \%$$\end{document}30% when the turbulent kinetic energy production rates of ST and CT are comparable. The mixing parameterization schemes used in ocean general circulation models (OGCMs) fail to accurately reproduce the mixing due to the pure turbulences and/or the nonlinear effects, and thus the mixed layer depth (MLD). Because analysis using global datasets suggests that nonlinear effects are large at the mid-latitudes, these results indicate that the nonlinear effects might be responsible for the deep MLD biases in OGCMs and that mixing parameterization schemes need to be improved to correctly represent ocean surface mixing due to shear and convection, as well as waves, in OGCMs.

Large-eddy simulations were performed to investigate the entrainment buoyancy flux at the mixed layer base due to nonlinearly interacting shear-driven turbulence (ST) and convective turbulence (CT). The fluxes due to pure ST and pure CT were first evaluated, and their scalings were derived. The entrainment flux due to coexisting ST and CT was then evaluated and compared to the scalingbased fluxes due to the pure turbulences. It was found that nonlinear effects reduce the entrainment flux by 30% when the turbulent kinetic energy production rates of ST and CT are comparable. The mixing parameterization schemes used in ocean general circulation models (OGCMs) fail to accurately reproduce the mixing due to the pure turbulences and/or the nonlinear effects, and thus the mixed layer depth (MLD). Because analysis using global datasets suggests that nonlinear effects are large at the mid-latitudes, these results indicate that the nonlinear effects might be responsible for the deep MLD biases in OGCMs and that mixing parameterization schemes need to be improved to correctly represent ocean surface mixing due to shear and convection, as well as waves, in OGCMs.
Vertical turbulent mixing induced by wind, surface cooling, and surface waves forms vertically uniform surface mixing/mixed layer (ML) in the upper stratified ocean. As the mixing deepens the ML, water below the ML is entrained into the layer, and ML water properties are changed. For example, the entrainment of colder and nutrient richer water changes temperature and the concentration of nutrients in the ML and then affects subsequent air-sea interaction 1,2 and primary production 3,4 , respectively. It is in fact suggested that the ML deepening due to the entrainment enhances the variability of Pacific Decadal Oscillation and then affects the regional and global climate systems 2 .
In ocean general circulation models (OGCMs), mixing parameterization schemes such as the Mellor-Yamada (MY) scheme 5 and the K-profile parameterization (KPP) scheme 6 are adopted to represent entrainment. Because the parameterization of the entrainment buoyancy flux at the ML base ( P b ) is key to developing the mixing parameterization, the flux has been investigated in several studies [7][8][9][10] . Studies of the atmospheric boundary layer (ABL) [11][12][13][14] are also useful because the entrainment process in the ocean surface boundary layer is dynamically similar to that in the ABL. Many ABL studies 15,16 focused on the entrainment process caused by convective turbulence (CT) because it is likely dominant in the ABL.
In the ocean, CT is considered to be dominant if the sea surface is severely cooled. Deep convection down to 1000 m depth and greater in Labrador and Greenland seas is an example. When CT is dominant, temperature, salinity and momentum become almost uniform in the vertical. Under purely convective forcing (no wind forcing), the entrainment buoyancy flux at the ML base ( P C b ) in the ocean is considered to be proportional to the surface buoyancy flux ( B f , which is defined as positive for sea surface cooling) and is expressed as Here, n ∼ = 0.2 9,13,14 when convection is shallow. If, on the other hand, convection is deep and turnover time of convection L MLD /W * is longer than 1/f, where W * [≡ (B f L MLD ) 1/3 ] is convective velocity scale, L MLD is ML depth (MLD), and f is the Coriolis parameter, the Earth rotation (the Coriolis acceleration term) acts and inhibits vertical velocity and hence CT 17 . Thus, the Earth's rotation decreases n at the convective Rossby number [18][19][20] .
If, on the other hand, the sea surface is weakly cooled as in low-latitude ( < 10 • ) regions, wind-induced sheardriven turbulence (ST) becomes dominant. When ST is dominant, temperature and salinity are well homogenized However, previous ABL studies showed that the ABL structure with coexisting ST and CT is different from the structure with each pure turbulence [24][25][26] , and thus CT is suppressed in the entrainment zone (corresponding to the ML base in the ocean) by ST 14 . These results indicate that the effects of nonlinear interaction between ST and CT may also be large in the surface ocean, and thus the entrainment buoyancy flux P b at the ML base cannot be expressed as a linear combination of P S b and P C b . Nevertheless, the nonlinear effects were not evaluated under realistic ocean surface forcing, and it is not clear when and where the nonlinear effects, if any, are large.
It is well known that recent OGCMs still have serious biases in simulating MLDs 27 . In recent decades, much attention has been directed to the shallow MLD biases in OGCMs and the effects of surface waves on these biases 27,28 . Surface waves interact with the wind-driven flow shear to form secondary circulations (Langmuir circulations) that induce turbulence and deepen the ML [29][30][31] . Non-breaking surface waves without wind-driven flow may also cause turbulence and deepen the ML 32,33 . Most OGCMs do not include these surface wave effects, and this is considered as one of the main reasons for the shallow MLD biases. Note, however, that deep MLD biases have also been found; the MLDs in OGCMs are sometimes greater than the observed values in regions such as the mid-latitudes (see Fig. 1 of Belcher et al. 28 and Fig. 11 of Tsujino et al. 34 ), even though most OGCMs omit surface wave effects. This fact clearly demonstrates that entrainment processes due to ST and CT also need to be re-investigated quantitatively.
The aim of this study is to evaluate the entrainment buoyancy flux at the ML base induced by nonlinearly interacting ST and CT and quantify the nonlinear effects in the surface ocean. Surface wave effects are not considered here to isolate this interaction processes from other complicated processes. To this end, large-eddy simulations (LESs) of the upper-ocean turbulence forced by uniform steady wind stress and/or cooling were performed; the configuration is described in "Methods" section. Uniform and steady surface forcing is used as a first step to understand nonlinear interaction between ST and CT. In "Results" section, we first evaluated the parameter dependences of the entrainment buoyancy flux due to pure ST and pure CT and derived the scaling of each type of turbulence. Then, we quantified the nonlinear effects by comparing the entrainment flux due to coexisting ST and CT with the fluxes due to each pure turbulence under the realistic ocean forcing parameters. The mixing parameterization schemes used in OGCMs were also tested to see whether the entrainment buoyancy flux due to nonlinear effects as well as the pure turbulence is reproduced. The global distribution of the intensity of the nonlinearity is presented using global datasets in "Discussion" section.

Results
In this section, we first show simulated profiles of horizontally averaged flow, buoyancy, and TKE tendency terms in typical cases. Then, we evaluate the simulated entrainment buoyancy flux due to pure ST and pure CT to obtain their scalings ( P S b and P C b , respectively). Using the scalings, we evaluate the simulated entrainment buoyancy flux due to coexisting ST and CT ( P b ) and compare P b with P S b + P C b to quantify the nonlinear interaction effects between ST and CT. Finally, the mixing parameterization schemes are tested to evaluate the extent to which they reproduce the entrainment flux of the pure turbulences and the nonlinear effects.
Profiles of horizontally averaged velocity, buoyancy, and TKE tendency terms in the ML: Typical cases. First, the results of typical cases of pure ST, pure CT, and coexisting ST and CT are shown. In this subsection, the initial stratification ( N 0 = 2.0 × 10 −2 s −1 ), initial MLD ( L 0 = L D /4 , where L D is the domain size described in "Methods" section), and Coriolis parameter ( f = 10 × 10 −5 s −1 ) are unchanged. Figure 1 shows the temporal depth variation of the horizontally averaged current speed and buoyancy and vertical profiles of the TKE tendency terms [see Eq. (21) in "Methods" section] averaged over 4.0 < t/T f < 5.0 , where t is time, and T f ( = 2π/f ) is the inertial period. In the simulation of pure ST ( U 2 * = 1.0 × 10 −4 m 2 s −2 and B f = 0 ), the horizontal mean current speed oscillated with the inertial period (Fig. 1a). The entrainment of less buoyant water into the ML decreased the buoyancy in the ML, and the MLD ( L MLD ) increased with time (Fig. 1d). Here, the MLD was defined as the depth at which the buoyancy production ( P b ), which is a TKE tendency term [Eq. (21)], was minimum 9 . Shear production ( P s ) was the dominant source of the pure ST, and the buoyancy production ( P b ) and dissipation ( D s ) were the sink terms in the ML (Fig. 1g). The convergence of the vertical transport of the TKE ( P t ) is positive in the lower ML, but it is much smaller than P s in this case.
In the simulation of pure CT ( U 2 * = 0 and B f = 19.6 × 10 −8 m 2 s −3 ), the MLD increased with time, although the horizontal mean current speed was zero (Fig. 1b). The buoyancy in the ML decreased because of entrainment and surface buoyancy flux (Fig. 1e). Because there was no horizontal mean velocity shear, P s was zero (Fig. 1h). P b was positive in the upper ML, whereas it was negative in the lower ML. P t was negative in the upper ML and positive in the lower ML, resulting in downward transport of the TKE in the ML.
In the simulation of coexisting ST and CT ( U 2 * = 1.0 × 10 −4 m 2 s −2 and B f = 19.6 × 10 −8 m 2 s −3 ), the current became more vertically uniform in the ML than that in the pure ST simulation (Fig. 1a and c). P s was positive in the ML, as in the pure ST simulation ( Fig. 1g and i). In the lower ML, P t was also positive, as in the pure CT simulation, and the contribution of P t to the TKE tendency became larger at the MLD, in contrast to that in the pure ST simulation (Fig. 1g-i). The vertical profiles of P s and P t in this coexisting turbulence simulation differ  www.nature.com/scientificreports/ from the linear combinations of those in the pure turbulence simulations, indicating that ST and CT interact nonlinearly with each other, as described in previous studies [24][25][26] . This result implies that the entrainment buoyancy flux ( P b ), which corresponds to the buoyancy production rate ( P b ) at the MLD, due to coexisting ST and CT is also not represented by the linear combination of those fluxes induced by pure ST and pure CT.
TKE tendency terms at the ML base and their scalings for pure ST and pure CT. In this subsection, we evaluate the TKE tendency terms at the ML base for pure ST and pure CT to derive scalings of the entrainment buoyancy flux for these two pure cases ( P S b and P C b , respectively). The scalings are used to quantify the nonlinear effects in coexisting ST and CT in the next subsection.
Pure ST case. The parameter dependence of the TKE tendency terms at the MLD in 135 pure ST simulations are reported here. Figure 2 shows a scatter plot of Ro ( ≡ U * /fL MLD ) and each TKE tendency term at the MLD ( P s , P t , −P b , and −D s ) normalized by U 3 * /L MLD averaged over 2.5 < t/T f < 3.5 and 4.0 < t/T f < 5.0 . Here, Ro and the TKE tendency terms at the MLD ( P s , P t , −P b , and −D s ) were sampled every T f /40 and averaged for T f . All these normalized tendencies decrease with decreasing Ro , except the normalized P t at Ro greater than 3, where it is almost constant. These decreases are especially rapid at Ro < 3 . Although P t is much smaller than P s in the typical case (where Ro ∼ 8 ) in the previous subsection, it become comparable at Ro ∼ 3 and larger at Ro 3.
The symbols and colors in Fig. 2 show the differences in initial MLD ( L 0 ) and N 0 /f , respectively. Note that despite the large number of simulations (135), the normalized tendency terms at the same Ro with different L 0 and N 0 /f collapse onto a single line. Although previous studies 35,36 suggested that the TKE tendency terms without the Earth's rotation depend on is the vertical buoyancy difference across the MLD) and/or Fr ≡ U * /N 0 L MLD , we found that stratification had little effect on the TKE tendency terms in the present parameter range.
By least-square fitting, we obtained the scalings of P S s and P S t as (4) P S s = 0.33Ro exp − www.nature.com/scientificreports/ where P S s and P S t are scaling-based tendencies for ST (while P s , P t , −P b , and −D s are simulation-based tendency) and there functional forms were determined as described in "Methods" section. The scaling of P b (i.e., P S b ) was derived from the relationship between P t /P s and P b /P s . The flux Richardson number R f (≡ −P b /P s ) at the MLD seems to be constant if P s ≫ P t 37,38 , but it is expected to change as the rotation effect increases because P s and P t became comparable at Ro 3 ( Fig. 2a and b). From Fig. 3, the relationship between P t /P s and R f was obtained as that is, Finally, the scaling of D s (i.e., D S s ) was obtained as where we assumed almost steady TKE at the MLD (due to slow deepening of the ML). These relations reproduce well the simulated dependence of the TKE tendency terms on Ro (Fig. 2).

Ro
Pure CT case. For pure CT, the TKE tendencies are suggested to depend on B f and Ro b ≡ W * /fL MLD [18][19][20] . Figure 4 shows a scatter plot of Ro b and each TKE tendency term (except P s , because P s = 0 in the pure CT cases) normalized by B f in 50 pure CT simulations. At Ro b > 3 , all the normalized tendency terms are almost constant, and P C b ∼ = −0.2B f , which is consistent with previous studies 9, 19 . At Ro b < 3 , they decrease with decreasing Ro b . Because the scaling of Wang 19 (Fig. 4b), the scaling was modified in this study. By least-square fitting, we obtained where P C t , P C b , and D C s are scaling-based tendencies for CT. Here, P C t + P C b + D C s � = 0 because the temporal change in the TKE is large (i.e., ongoing ML deepening occurs), in contrast to that in the pure ST cases.   Figure 5a shows P s /P S s (ratio of the simulated shear production with CT to the scaling-based shear production due to pure ST) as a function of P C b /P S b (a measure of CT relative to ST). Note that P s /P S s was excluded if P t /P s > 100 because P s has little effect on TKE tendency. P s /P S s increases with P C b /P S b , indicating that ST is intensified by CT. For example, P s becomes 20-30 times larger than P S s at P C b /P S b ∼ = 10 2 . Note, however, that as CT becomes more dominant, P t has more impact on the TKE tendency than P s [see the color representing the intensity of P t relative to P s ( P t /P s ) in Fig. 5a]. Figure 5b shows −(P s − P S s )/P b , a ratio of the increased shear production by CT ( P s − P S s ) to the entrainment buoyancy flux ( P b ), as a function of P C b /P S b . The impact of the increased shear by CT on the entrainment is found largest at around P C b /P S b = 10 1/2 . Previous studies 14,39 , on the other hand, suggested that CT is inhibited by ST. To see this effect, we plot P t /(P S t + P C t ) as a function of P C b /P S b in Fig. 5c. Here, P t is a measure of CT at the MLD, and the ratio of P t to P S t + P C t represents the intensity of the nonlinear interaction effects. P t was 40 % smaller than P S t + P C t at 10 −1/2 < P C b /P S b < 10 1/2 , indicating that ST and CT interact nonlinearly to decrease P t in this parameter range. This indicates that ST disrupts coherent structure of pressure and/or TKE and vertical velocity [see Eq. (21) in "Methods" section] associated with convective motion (CT). This P t decrease [ P t − P S t + P C t ] by the interaction contributes more to P b than the increased P s by the interaction (Fig. 5b and d). Consequently, P b became 30 % smaller than P S b + P C b at 1 < P C b /P S b < 10 3/2 (Fig. 5e).

Entrainment flux in the ocean mixing parameterization schemes used in OGCMs.
In the previous subsection, we suggested that the nonlinear interaction between ST and CT likely inhibits the entrainment at the ML base at 1 < P C b /P S b < 10 3/2 . To accurately simulate the ML-related processes using OGCMs, the mixing parameterization schemes should reproduce the entrainment buoyancy flux of the nonlinear interaction effects as well as those of the pure turbulences. To see the extent to which the schemes reproduce the entrainment flux of the pure turbulences and the nonlinear effects, one-dimensional (1D) simulations with the mixing parameterization schemes were performed, and the results are compared to the LES results in this subsection.
The mixing parameterization schemes tested in this study were the KPP scheme 6 , level 2.5 MY scheme 5,40,41 , and level 2.5 Nakanishi-Niino (NN) scheme [42][43][44] . The NN scheme is a modified version of the MY scheme and is used in ocean models 45 as well as atmospheric models 46,47 . The boundary and initial conditions were the same as those in the LES. The detail configuration for 1D simulations is described in "Methods" section. Figure 6 shows scatter plots of Ro ( ≡ U * /fL MLD ) and P S b normalized by U 3 * /L MLD for the 1D pure ST simulations (cf. Fig. 2c), Ro b (≡ W * /fL MLD ) and P C b normalized by B f for the 1D pure CT simulations (cf. Fig. 4c), and P C b /P S b and P b normalized by P S b + P C b in the 1D simulations of coexisting ST and CT (cf. Fig. 5c). Figure 6 also shows the temporal change in MLD in the 1D simulations using the typical parameters used in the LESs shown in Fig. 1. For pure ST, the normalized P S b s in the KPP and MY schemes show more scattering than those in the LES at Ro > 3 , suggesting that these schemes are likely affected by other factors such as Ri * and/or Fr , which had little effect on P C b in the LESs in the present parameter range (Fig. 6a). Because these normalized P s b s were also underestimated, the ML deepened less from t/T f = 1 to t/T f = 4.5 than it did in the LES (Fig. 6d). The NN scheme successfully reproduced the dependence of P S b on Ro with less scatter, but it overestimated P S b (Fig. 6a) and thus the MLD (Fig. 6d).
The scalings of P S b in these schemes were evaluated for later use. For simplicity, P S b was assumed to be proportional to Ro d , where d is constant. By least-square fitting, we obtained  www.nature.com/scientificreports/ For pure CT, on the other hand, the normalized P C b s in the KPP and NN schemes were similar to those of the LES, whereas those in the MY scheme were much smaller (Fig. 6b). As a result, the MLDs were well reproduced by the KPP and NN schemes and underestimated by the MY scheme (Fig. 6e). The decreases in P C b with decreasing Ro b were smaller in all of these schemes than in the LES (Fig. 6b), probably because they do not include the effects of Earth's rotation. Here, we assume P C b ∝ B f and derive the scalings of P C b in these schemes by least-square fitting as Pt / Ps < 10 3 < Pt / Ps < 100 10 < Pt / Ps 100 < Figure 5. Scatter plots of (a) P C b /P S b and P s /P S s , (b) P C b /P S b and (P s − P S s )/P b , (c) P C b /P S b and P t /(P S t + P C t ) , (d)  (Fig. 6c). The P b s in the MY and NN schemes are greater than P S b + P C b , indicating that they tend to represent the interaction effects in an opposite sense. The MLD was underestimated by the MY scheme (Fig. 6f) owing to the underestimations of P S b (Fig. 6a) and P C b (Fig. 6b), despite the failure to reproduce the nonlinear effects (Fig. 6c). On the other hand, the MLD was overestimated by the NN scheme (Fig. 6f) because of the overestimation of P S b (Fig. 6a) and the failure to reproduce the nonlinear effects (Fig. 6c). The KPP scheme successfully reproduced the nonlinear effects except at P C b /P S b > 10 (Fig. 6c), although the underestimation of P S b resulted in the underestimation of the MLD (Fig. 6f). Note that the above differences between the LES and mixing schemes would become smaller by tuning empirical parameters in the schemes, though it seems uneasy to reduce the differences in the pure and coexisting turbulence regimes simultaneously.

Discussion
In this section, we estimate the global distribution of the intensity of the nonlinear effects in the real ocean. The parameters Ro , W 3 * /U 3 * , and Ro b were calculated from the U * , B f , and L MLD in observed and reanalysis data in autumn and winter (see "Methods" section) and are shown in Fig. 7. (Note that the U 3 * /L MLD -normalized P S b and P C b depend only on these parameters.) Here, data at B f < 0 were excluded from the analysis. (These data were often found near the equator.) Figures 7a and b show that Ro is larger due to smaller f and amounts to 10 1/2 (= 3.2) or more in tropical regions ( < 20 • ). On the other hand, W 3 * /U 3 * is small ( < 10 1/2 ) at lower latitude than 15 • and higher latitude than 40 • and large ( > 10 1/2 ) at mid-latitude ( 15 • -40 • ) ( Fig. 7c and d). As a result, ST is more dominant than CT ( P C b /P S b < 1 ) at the low and high latitudes ( Fig. 7g and h). At higher latitude than 45 • N in the North Atlantic, W 3 * /U 3 * is also large (Fig. 7c and d), and hence P C b /P W b > 10 3/2 ( Fig. 7g and h), especially in winter. However, Ro b < 1 in winter ( Fig. 7e and f) suggests that deep water formation in the North Atlantic is inhibited by Earth's rotation 17 . Note that at mid-latitudes, P C b /P S b ranges from 1 to 10 3/2 , indicating that the nonlinear interaction between ST and CT is expected.
To estimate geological distribution of the expected nonlinear interaction intensity ( Fig. 7i and j), we used Fig. 8, where scatter plots of Ro and W 3 * /U 3 * from the observed (Fig. 8a) and simulated (Fig. 8b) data are shown. In Fig. 8b, the intensity of the nonlinear effects [ P b /(P S b + P C b ) ] averaged over bins on (Ro, W 3 * /U 3 * ) space is also shown. Here, the observed (Ro, W 3 * /U 3 * ) at a certain grid point (Fig. 7a-d and 8a) was converted to P b /(P S b + P C b ) using simulated relation between (Ro, W 3 * /U 3 * ) and bin-averaged P b /(P S b + P C b ) in Fig. 8b, and this was considered as observed P b /(P S b + P C b ) . The cross-hatching in Fig. 7i and j represents the region where Ro and W 3 * /U 3 * are outside of the simulated parameter range. These figures show that our simulations covered most of the observed pairs of Ro and W 3 * /U 3 * , except in the Southern Ocean, where surface cooling is weak relative to wind stress ( W 3 * /U 3 * < 0.1 ) and ST is expected to strongly dominate CT. The observed P b /(P S b + P C b ) is less than 0.8 at mid-latitudes and 0.6 at some region between 15 • and 25 • , indicating that the nonlinear interaction between ST and CT probably suppresses the entrainment at the ML base there. Because some mixing parameterization schemes such as the MY and NN schemes cannot reproduce the nonlinear effects, the fact that the schemes did www.nature.com/scientificreports/ not successfully reproduce the effect of nonlinear interaction between ST and CT might explain the mid-latitude deep MLD biases in winter observed in some OGCMs as seen in Fig. 1 of Belcher et al. 28 and Fig. 11 of Tsujino et al. 34 .
[More than half OGCMs in the Coupled Model Intercomparison Project phase 6 evaluated by Tsujino et al. 34 adopted schemes similar to the MY and NN schemes (1.5 or higher order turbulence closure schemes), though each of the schemes used slightly different parameterizations and/or tuning parameters from those of the MY and NN schemes.] This result suggests that mixing parameterization schemes need to be checked and improved (if necessary) to correctly represent ocean surface mixing due to ST and CT.
In this study, we found the nonlinear interaction between ST and CT is expected large in mid-latitude ML in the ocean. However, the interaction mechanism remains to be investigated in more detail. We also found that the KPP, MY, and NN schemes do not well represent pure ST mixing. Because ST likely plays more role in the ocean than in the atmosphere, we consider that this issue should not be overlooked in ML mixing schemes. Effects of wave and/or time-variying forcing as well as heterogeneous background environment (such as ocean front) also need to be considered simultaneously for realistic ML simulation in the OGCMs. These will be studied in future.

Methods
Simulations and data. Numerical model and experimental configurations for large-eddy simulations. The LES model used in this study is the same as that used in Ushijima and Yoshikawa 48,49 . The governing equations are the momentum equation, continuity equation, and advection-diffusion equation of buoyancy under the incompressible, f-plane, Boussinesq, and rigid-lid approximations. Subgrid-scale parameterization follows the method described by Deardorff 50 and Maronga et al. 51 . At the surface, constant wind stress ( ρ 0 U 2 * , where ρ 0 = 1.0 × 10 3 kg m −3 is the reference water density) and buoyancy flux ( B f ) were imposed. We also imposed subgrid-scale shear production at the surface. At the bottom, the free-slip condition and no-buoyancy flux condition were imposed. The lateral boundaries were periodic in both directions. The initial stratification  www.nature.com/scientificreports/ ST and CT, L D was set to 4L P73 (1 + 5.1B f /U 2 * N 0 ) , where L P73 ≡ U * / N 0 f is the MLD scale characterizing the wind-induced ML in the stratified ocean under the Earth's rotation 49,52 , whereas L D was set to 17 B f /N 2 f in the simulations of pure CT. We performed several LESs with quarter-grid spacing but the same domain size and found that the TKE tendency terms obtained in the simulations with higher resolution were almost the same as those obtained with the original resolution. The dependence on the resolution is discussed in detail in Supplementary Information. The integration was continued for 5T f , where T f = 2π/f is the inertial period. Observed and reanalysis data. Data were analyzed for autumn (October, November, and December in the northern hemisphere and April, May, and June in the southern hemisphere) and winter (January, February, and March in the northern hemisphere and July, August, and September in the southern hemisphere), when ST and CT are typically expected to coexist. The climatology of the ML temperature ( T ML ) and salinity ( S ML ) as well as the MLD ( L MLD ) of the mixed layer Argo dataset, gridpoint value (MILA-GPV) 53 were used. The surface fluxes were the 6-hourly momentum fluxes ( τ x , τ y ), net heat flux ( H f ), and freshwater flux ( E − P ) from the National Centers for Environmental Prediction (NCEP) data 54 for 2001-2010, where E and P are the evaporation rate and precipitation rate, respectively. The shortwave radiation was assumed not to penetrate below the surface for simplicity, and the evaporation rate was estimated from the latent heat flux with the latent heat vaporization of water 55 . These fluxes were converted to the friction velocity U * [= (τ 2 x + τ 2 y ) 1/4 /ρ 1/2 0 ] and buoyancy flux B f [= −αgH f /ρ 0 C a + βg(E − P)S ML ] . Here, ρ 0 (= 1.0 × 10 3 kg m −3 ) and C a (= 4.0 × 10 3 J kg −1 C −1 ) are the reference density and heat capacity of water, respectively. The thermal expansion rate ( α ) and haline contraction rate ( β ) were calculated from T ML and S ML using the equation of state for seawater 56 . The momentum flux, buoyancy flux, and MLD were averaged into seasonal (three-month) climatological data points. The horizontal resolution of the data was 5 • × 5 • .
Analysis. TKE tendency terms in the LES. The TKE tendency terms were calculated as where u i represents the velocity components (u, υ, w) in the x i direction, x i (i = 1, 2, 3) denotes the Cartesian coordinates (x, y, z), b is buoyancy, π = p + 2ρ 0 e/3 is modified pressure, p is pressure, and s ij ≡ (∂u i /∂x j + ∂u j /∂x i )/2 . The subgrid-scale kinetic energy (e), eddy viscosity ( ν ), eddy diffusivity ( κ ), and dissipation rate ( ε ) were calculated using sub-grid scale parameterization 50,51 . The overbar represents the horizontal average, and the prime indicates anomalies from the horizontal average. In the above equation, P s , P b , P t , and D s represent the rates of shear production, buoyancy production, convergence of vertical transport of the TKE, and dissipation of the TKE, respectively. (Note that they are a function of z.) Functional forms of the P S s and P S t scalings. In the pure ST simulations, the normalized P s is almost linearly proportional to Ro at large Ro , but the slope increases for smaller Ro (Fig. 2a). Under neutral stratification, the vertical shear of the horizontal velocity (the Ekman velocity shear) decreases with depth (|z|) as exp(−|z|/L EKD ) , where L EKD ≡ U * /f is the depth of the turbulent Ekman layer 21 . Therefore, the vertical shear at the MLD is expected to be proportional to exp[−cL MLD /(U * /f )] = exp(−c/Ro) , where c is a constant. Consequently, we assume that the normalized P s is proportional to Ro exp(c/Ro).
On the other hand, the normalized P t at Ro < 3 decreases with decreasing Ro but does not vary significantly with Ro at Ro > 3 (Fig. 2b). Because the variation of this normalized P S t with Ro is similar to that of the ∂ ∂t www.nature.com/scientificreports/ normalized P C t with Ro b in the pure CT simulations, we assume that the normalized P S t has the same functional form as the normalized P C t [Eq. (9)].