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 [Formula: see text] 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.

Introduction

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 () is key to developing the mixing parameterization, the flux has been investigated in several studies[7-10]. Studies of the atmospheric boundary layer (ABL)[11-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 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 () in the ocean is considered to be proportional to the surface buoyancy flux (, which is defined as positive for sea surface cooling) and is expressed asHere, [9,13,14] when convection is shallow. If, on the other hand, convection is deep and turnover time of convection is longer than 1/f, where is convective velocity scale, 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-20].

If, on the other hand, the sea surface is weakly cooled as in low-latitude () regions, wind-induced shear-driven turbulence (ST) becomes dominant. When ST is dominant, temperature and salinity are well homogenized in the vertical while momentum is not, because the turbulent kinetic energy (TKE) in the ML is produced by the vertical shear of the horizontal velocity at each depth. At the ML base, the TKE is converted to the potential energy, and the entrainment buoyancy flux () is expressed aswhere is the friction velocity and m is a proportional coefficient[8]. It should be noted that the depth scale of the wind-induced shear is given by (turbulent Ekman scale)[21] and the shear becomes almost zero at the depth below in neutrally stratified fluid. This indicates that the shear at the ML base is weakened (intensified) if is greater (smaller) than . Consequently, m in Eq. (2) should depend on the Rossby number [10,22]. However, the dependence of m on has not been well investigated.

In autumn and winter, especially at mid-latitudes, ST and CT typically coexist in the ML. Some previous studies[9,10,23] assumed for simplicity that the entrainment buoyancy flux at the ML base due to coexisting ST and CT () can be expressed as the linear combination of the entrainment buoyancy fluxes due to pure ST and pure CT, that is,However, previous ABL studies showed that the ABL structure with coexisting ST and CT is different from the structure with each pure turbulence[24-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 at the ML base cannot be expressed as a linear combination of and . 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-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 ( and , respectively). Using the scalings, we evaluate the simulated entrainment buoyancy flux due to coexisting ST and CT () and compare with 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 (), initial MLD (, where is the domain size described in "Methods" section), and Coriolis parameter () are unchanged.

Figure 2

Scatter plots of and (a) , (b) , (c), and (d) in the pure ST simulations () averaged over and . Each TKE tendency term was normalized by . Symbols represent the initial MLD (). Colors represent initial stratification (). Dashed lines are the scalings derived in this study [Eqs. (4), (5), (7), and (8)].

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 , where t is time, and () is the inertial period. In the simulation of pure ST ( and ), 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 () increased with time (Fig. 1d). Here, the MLD was defined as the depth at which the buoyancy production (), which is a TKE tendency term [Eq. (21)], was minimum[9]. Shear production () was the dominant source of the pure ST, and the buoyancy production () and dissipation () were the sink terms in the ML (Fig. 1g). The convergence of the vertical transport of the TKE () is positive in the lower ML, but it is much smaller than in this case.

Figure 1

Time-depth variation of the horizontally averaged (a)–(c) current speed and (d)–(f) buoyancy (B) and (g)–(i) the vertical profiles of the TKE tendency terms averaged over in the (a), (d), (g) pure ST (), (b), (e), (h) pure CT (), and (c), (f), (i) coexisting ST and CT () simulations. Initial stratification and initial MLD are and , respectively. Solid lines in (a)–(f) represent the MLD.

In the simulation of pure CT ( and ), 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, was zero (Fig. 1h). was positive in the upper ML, whereas it was negative in the lower ML. 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 ( and ), the current became more vertically uniform in the ML than that in the pure ST simulation (Fig. 1a and c). was positive in the ML, as in the pure ST simulation (Fig. 1g and i). In the lower ML, was also positive, as in the pure CT simulation, and the contribution of 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 and in this coexisting turbulence simulation differ 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-26]. This result implies that the entrainment buoyancy flux (), which corresponds to the buoyancy production rate () 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 ( and , 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 () and each TKE tendency term at the MLD (, , , and ) normalized by averaged over and . Here, and the TKE tendency terms at the MLD (, , , and ) were sampled every and averaged for . All these normalized tendencies decrease with decreasing , except the normalized at greater than 3, where it is almost constant. These decreases are especially rapid at . Although is much smaller than in the typical case (where ) in the previous subsection, it become comparable at and larger at .

The symbols and colors in Fig. 2 show the differences in initial MLD () and , respectively. Note that despite the large number of simulations (135), the normalized tendency terms at the same with different and collapse onto a single line. Although previous studies[35,36] suggested that the TKE tendency terms without the Earth’s rotation depend on (where is the vertical buoyancy difference across the MLD) and/or , 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 and aswhere and are scaling-based tendencies for ST (while , , , and are simulation-based tendency) and there functional forms were determined as described in "Methods" section. The scaling of (i.e., ) was derived from the relationship between and . The flux Richardson number at the MLD seems to be constant if [37,38], but it is expected to change as the rotation effect increases because and became comparable at (Fig. 2a and b). From Fig. 3, the relationship between and was obtained asthat is,Finally, the scaling of (i.e., ) was obtained aswhere 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 (Fig. 2).

Figure 3

Scatter plot of and the flux Richardson number in the pure ST simulations () averaged over and . Symbols and colors are the same as in Fig. 2. Dashed line is the regression line [Eq. (6)].

Figure 4

Scatter plots of and (a) , (b), and (c) in the pure CT simulations () averaged over and . Each TKE tendency term was normalized by . Symbols and colors are the same as in Fig. 2. Dashed lines are the scalings derived in this study [Eqs. (9)–(11)], and solid line in (b) is the scaling of[19] [].

Figure 5

Scatter plots of (a) and , (b) and , (c) and , (d) and , and (e) and in the coexisting ST and CT simulations. Symbols are the same as in Fig. 2. Colors represent . Black circle and bar show the averaged value and standard deviation, respectively.

Pure CT case

For pure CT, the TKE tendencies are suggested to depend on and [18-20]. Figure 4 shows a scatter plot of and each TKE tendency term (except , because in the pure CT cases) normalized by in 50 pure CT simulations. At , all the normalized tendency terms are almost constant, and , which is consistent with previous studies[9,19]. At , they decrease with decreasing . Because the scaling of Wang[19] [] slightly overestimated at (Fig. 4b), the scaling was modified in this study. By least-square fitting, we obtainedwhere , , and are scaling-based tendencies for CT. Here, because the temporal change in the TKE is large (i.e., ongoing ML deepening occurs), in contrast to that in the pure ST cases.

Nonlinear effects of ST and CT interaction at the ML base

Using the scalings of the TKE tendency terms for pure ST ( [Eq. (4)], [Eq. (5)], and [Eq. (7)]) and pure CT ( [Eq. (9)] and [Eq. (10)]) and simulation-based tendency of , and for 450 simulatinos of coexisting ST and CT, the nonlinear effects of ST and CT interaction at the MLD are evaluated. Figure 5a shows (ratio of the simulated shear production with CT to the scaling-based shear production due to pure ST) as a function of (a measure of CT relative to ST). Note that was excluded if because has little effect on TKE tendency. increases with , indicating that ST is intensified by CT. For example, becomes 20–30 times larger than at . Note, however, that as CT becomes more dominant, has more impact on the TKE tendency than [see the color representing the intensity of relative to () in Fig. 5a]. Figure 5b shows , a ratio of the increased shear production by CT () to the entrainment buoyancy flux (), as a function of . The impact of the increased shear by CT on the entrainment is found largest at around .

Previous studies[14,39], on the other hand, suggested that CT is inhibited by ST. To see this effect, we plot as a function of in Fig. 5c. Here, is a measure of CT at the MLD, and the ratio of to represents the intensity of the nonlinear interaction effects. was smaller than at , indicating that ST and CT interact nonlinearly to decrease 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 decrease [] by the interaction contributes more to than the increased by the interaction (Fig. 5b and d). Consequently, became smaller than at (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 . 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-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 () and normalized by for the 1D pure ST simulations (cf. Fig. 2c), and normalized by for the 1D pure CT simulations (cf. Fig. 4c), and and normalized by 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 s in the KPP and MY schemes show more scattering than those in the LES at , suggesting that these schemes are likely affected by other factors such as and/or , which had little effect on in the LESs in the present parameter range (Fig. 6a). Because these normalized s were also underestimated, the ML deepened less from to than it did in the LES (Fig. 6d). The NN scheme successfully reproduced the dependence of on Ro with less scatter, but it overestimated (Fig. 6a) and thus the MLD (Fig. 6d).

Figure 6

Same as (a) Fig. 2c, (b) Fig. 4b, and (c) Fig. 5e but for the results of the LES (red) and 1D simulations with the KPP (blue), MY (green), and NN (orange) schemes. Temporal variation of the MLD in the simulations with the typical parameters of (d) pure ST, (e) pure CT, and (f) coexisting ST and CT. The parameters used in (d)–(f) are the same as those in Fig. 1a–c, respectively. Time is normalized by . Colors represent the LES (red) and the KPP (blue), MY (green), and NN (orange) schemes.

The scalings of in these schemes were evaluated for later use. For simplicity, was assumed to be proportional to , where d is constant. By least-square fitting, we obtainedFor pure CT, on the other hand, the normalized 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 with decreasing 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 and derive the scalings of in these schemes by least-square fitting asThe nonlinear interaction effects between ST and CT in these mixing parameterization schemes were quantified using these scalings [Eqs. (12)–(17)] (Fig. 6c). The s in the MY and NN schemes are greater than , 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 (Fig. 6a) and (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 (Fig. 6a) and the failure to reproduce the nonlinear effects (Fig. 6c). The KPP scheme successfully reproduced the nonlinear effects except at (Fig. 6c), although the underestimation of 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.

Figure 8

Scatter plots of and calculated from (a) observed data and (b) LES data. Symbols in (a) and (b) represent seasons and initial MLD (), respectively. Color represents latitude in (a) and bin-averaged in (b). Solid lines are contour lines of , and dash-dotted line in (a) is contour line of .

Discussion

In this section, we estimate the global distribution of the intensity of the nonlinear effects in the real ocean. The parameters , , and were calculated from the , , and in observed and reanalysis data in autumn and winter (see "Methods" section) and are shown in Fig. 7. (Note that the -normalized and depend only on these parameters.) Here, data at were excluded from the analysis. (These data were often found near the equator.) Figures 7a and b show that is larger due to smaller f and amounts to or more in tropical regions (). On the other hand, is small () at lower latitude than and higher latitude than and large () at mid-latitude (–) (Fig. 7c and d). As a result, ST is more dominant than CT () at the low and high latitudes (Fig. 7g and h). At higher latitude than in the North Atlantic, is also large (Fig. 7c and d), and hence (Fig. 7g and h), especially in winter. However, 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, ranges from 1 to , indicating that the nonlinear interaction between ST and CT is expected.

Figure 7

Seasonal (three month) averages of (a), (b) , (c), (d) , (e), (f) , (g), (h) , and (i), (j) in (a), (c), (e), (g), (i) autumn and (b), (d), (f), (h), (j) winter estimated from observations in 2001–2010. in (i) and (j) were evaluated from the observed and the simulated bin-averaged relationship between , and shown in Fig. 8b.

To estimate geological distribution of the expected nonlinear interaction intensity (Fig. 7i and j), we used Fig. 8, where scatter plots of and from the observed (Fig. 8a) and simulated (Fig. 8b) data are shown. In Fig. 8b, the intensity of the nonlinear effects [] averaged over bins on space is also shown. Here, the observed at a certain grid point (Fig. 7a-d and 8a) was converted to using simulated relation between and bin-averaged in Fig. 8b, and this was considered as observed . The cross-hatching in Fig. 7i and j represents the region where and are outside of the simulated parameter range. These figures show that our simulations covered most of the observed pairs of and , except in the Southern Ocean, where surface cooling is weak relative to wind stress () and ST is expected to strongly dominate CT. The observed is less than 0.8 at mid-latitudes and 0.6 at some region between and , 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 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 (, where is the reference water density) and buoyancy flux () 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 (, where B is the horizontally averaged buoyancy) was zero () from the surface () to the initial MLD () and constant () from to the bottom of the ocean (, where is the domain length).

To evaluate the buoyancy entrainment flux due to pure ST, pure CT, and coexisting ST and CT, simulations were performed with several values of the momentum flux (), surface buoyancy flux (), initial stratification (, and , corresponding to temperature changes of 0.078, 0.31, 1.3, 5, and in ), initial MLD (), and Coriolis parameter (f). In the simulations of pure ST, we set , and (corresponding to wind speeds at 10 m height of 7, 10, and ), , , and , and , and (corresponding to latitudes of , and ). In the simulations of pure CT, we set , , and (corresponding to surface cooling of 25, 50, 100, 200, and ), and , and . A total of 135 and 50 simulations were performed for pure ST and pure CT, respectively. To examine the nonlinear interaction between ST and CT, a total of 450 simulations with , and , , and , and , and , and were performed. These parameters are determined from the typical values of the observed climatologies in autumn and winter.

The governing equations were discretized using the second-order finite-difference scheme and integrated in time using the second-order Runge–Kutta scheme. The number of grid cells was , and the grid spacing was uniform. The domain size, , was varied according to the friction velocity (), buoyancy flux (), initial stratification (), and Coriolis parameter (f). In the simulations of pure ST and coexisting ST and CT, was set to , where is the MLD scale characterizing the wind-induced ML in the stratified ocean under the Earth’s rotation[49,52], whereas was set to 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 , where is the inertial period.

Experimental configurations for one-dimensional simulations with mixing parameterization schemes

The governing equations for one-dimensional (1D) simulations are momentum equation and diffusive equation of buoyancy,where (U, V) is the horizontal velocity components, B is buoyancy, and and are the eddy viscosity and diffusivity, determined in the KPP[6], MY[5,40,41], or NN[42-44] schemes, respectively. The boundary condition at the surface and the bottom, domain depth (), and the number of the vertical grid cells, are same as those in the LES. A total of 135, 50, and 450 simulations for pure ST, pure CT, and coexisting ST and CT were respectively performed for each 1D experiment with different mixing schemes with same momentum flux (), buoyancy flux (), initial stratification (), initial MLD (), and Colioris parameter as those in the LES.

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 () and salinity () as well as the MLD () of the mixed layer Argo dataset, gridpoint value (MILA-GPV)[53] were used. The surface fluxes were the 6-hourly momentum fluxes (), net heat flux (), and freshwater flux () 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 and buoyancy flux . Here, and are the reference density and heat capacity of water, respectively. The thermal expansion rate () and haline contraction rate () were calculated from and 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 .

Analysis

TKE tendency terms in the LES

The TKE tendency terms were calculated aswhere represents the velocity components in the direction, denotes the Cartesian coordinates (x, y, z), b is buoyancy, is modified pressure, p is pressure, and . 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, , , , and 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 and scalings

In the pure ST simulations, the normalized is almost linearly proportional to at large , but the slope increases for smaller (Fig. 2a). Under neutral stratification, the vertical shear of the horizontal velocity (the Ekman velocity shear) decreases with depth (|z|) as , where is the depth of the turbulent Ekman layer[21]. Therefore, the vertical shear at the MLD is expected to be proportional to , where c is a constant. Consequently, we assume that the normalized is proportional to .

On the other hand, the normalized at decreases with decreasing but does not vary significantly with at (Fig. 2b). Because the variation of this normalized with is similar to that of the normalized with in the pure CT simulations, we assume that the normalized has the same functional form as the normalized [Eq. (9)].

Supplementary Information 1.

Supplementary Information 2.