The LES used in this study was originally developed by Moeng (1984) and modified by P. Sullivan (e.g., Sullivan et al., 1996). The model had been applied to the equatorial ocean by Wang et al. (1996); Wang et al. (1998) and Wang and Müller (2002). The model employs a Fourier pseudospectral method in both horizontal directions and a second-order finite difference scheme in the vertical direction. The radiation conditions are applied to the bottom boundary, allowing downward propagating internal waves to leave the system (Klemp and Durran, 1983). Periodic boundary conditions are used in the horizontal directions.

The governing equations (Equations 1.1–1.4) are

where u = (u, v, w) is the velocity, p is the pressure (normalized by a reference density), α = 2.6×10⁻⁴ K⁻¹ is the thermal expansion coefficient, g is the gravitational acceleration vector, τ is the subgrid stress tensor, T is the potential temperature, q = (q1, q2, q3) is the subgrid heat flux, e is the subgrid-scale turbulent kinetic energy, ϵ is the turbulent kinetic energy dissipation rate, and K is a diffusion tensor. Detailed descriptions of discretization and subgrid-scale parameterization can be found in Sullivan et al. (1996) and Wang et al. (1996). Both vertical and horizontal components of the earth’s rotation are ignored.

The computational domain is 512 m × 512 m in the horizontal directions and 256 m in the vertical direction, respectively. The domain is discretized at Δx = Δy = Δz = 1 m. Such domain sizes and grid resolutions can resolve both the “long” scale (wavelength much larger than the size of turbulent eddies) oscillations/internal waves that are observed during field measurements (e.g., Moum et al., 1992) and the small overturning scales during the evolution of shear instabilities.

To investigate the relationship between turbulence strength and background variables, we designed 27 experiments with different initial conditions of velocity and temperature, resulting in different combinations of stratification, shear, and Ri.

In order to assess the properties of turbulent mixing induced by a sheared and density-stratified parallel flow, initial depth-dependent profiles for the horizontal velocity u0 and temperature T₀, and related N02, S02 and Ri0 are fixed to the idealized and dimensional profiles (Figure 1). The explicit expressions for u0 and squared buoyancy frequency N02 are given by

This flow structure follows that of Smyth and Peltier (1989). Here, corresponding A_i and B_j are the factors of u0 and N02, respectively. u0 decreases with depth slowly above 100 m with weak shear, and decreases with depth dramatically to nearly zero at about 150 m, resulting in large shear in-between with a squared shear peak of 2.24 × 10⁻⁴ s⁻² at 128 m. N02 has a maximum of 0.2 × 10⁻⁴ s⁻² also at 128 m; T₀ is obtained by integration of N02 using a T-dependent ocean state equation. Consequently, the minimum of Ri0, Rimin, reaches 0.2 at 128 m. Hereafter, the depth of 128 m is denoted as z0. Away from this stratified shear layer, Ri0 is very large: it increases from its minimum to larger than 2 above 100 m and below 150 m depths. Here, the layer with Ri0< 0.25 is defined as the initial unstable layer (hereafter IUL). Consequently, the flow within the IUL is unstable to KH instability, which ensures the generation of turbulence after small amplitude perturbations kickstart the instability (Smyth and Peltier, 1989; Smyth et al., 2005). We note that, though these profiles cannot fully capture all potential profiles of the unstable shear layers in the ocean, they represent a large part of the characteristics of shear-driven turbulence in the stratified ocean.

A total of 27 unstable flows, with different Rimin and IULs, have been designed to obtain sufficient turbulence properties. The parameters of the flows are listed in Table 1, where the aforementioned u0 and N02 profiles are for the experiment A₇B₇. Here we choose two other examples to depict our setting. For example, in experiment A₆B₆, the initial zonal velocity u0 in A₇B₇ is multiplied by a constant factor A₆ = 0.9, which results in a S02 that is 0.81 times that of experiment A₇B₇; at the same time, N02 is multiplied by B₆ = 0.81, thus the Ri0 profile of A₆B₆ remains the same as that of experiment A₇B₇. In experiment A₆B₅, the initial zonal velocity u0 in A₇B₇ is multiplied by a constant factor A₆ = 0.9, but N02 is multiplied by a constant factor B₅ = 0.64; as a result, the stratification weakens more than the shear, and the profile of Ri0 decreases, reaching a minimum of 0.05. Based on this rule, the other 24 profiles are designed. The factors to u0, named A₁–A₇, increase from 0.4 to 1.0, while the factors to N02, named B₁–B₇, increase from 0.16 to 1.0. The set of 27 profiles contains not only variable Ri0 ( Rimin ranging from 0.057 to 0.201) with constant N02 or S02, but also constant Ri0 with variable N02 or S02, making the resulting turbulence more ergodic and the afterward statistical analysis more flexible.

We note that, because our study is focused on the turbulent mixing in the interior ocean, the surface forcings, including both the wind-stress-induced friction velocity and the surface heat/buoyancy flux, are set to zero in the LES experiments. This avoids the influence of boundary forcing on turbulence just below the ML base (Zaron and Moum, 2009). In addition, large-scale forcing that represents the maintenance of the background flow via larger-scale motions (Wang et al., 1998) is not set, either. As such, each of our experiments documents a non-forced evolution of turbulence, which provides ‘pure’ KH instability-induced turbulence data; this contrasts with the observed turbulence which could result from more complex processes. Due to the absence of both the external forcing and large-scale forcing, the turbulence decays rapidly which usually lasts for less than 24 hours.

Since the instability develops from the IUL, the initial variables, N0¯, S0¯ and Ri0¯ are vertically averaged over the IUL. The thickness of IUL (IULT, denoted as h0) is also considered an important initial variable. They are used for subsequent calculation and parameterization.

As the energy sources and sinks of turbulent evolution (Figure 3a), within the TL and over the turbulent stage, the energy transferred via shear production ESPLES, the energy transferred via turbulent dissipation EϵLES and the energy transferred via buoyancy production EBLES can be directly calculated from the LES outputs by the following equations:

where tstart and tend are the start and end of the turbulent stage, and z₁ and z₂ are the upper and lower boundaries of the TL. Primes and tildes represent the deviations from the horizontal mean, and the horizontal average, respectively.

The 3D (horizontally over the domain and vertically over the TL) and temporally (over the turbulent stage) averaged TKE dissipation rate ϵLES, turbulent diffusivity κLES and buoyancy frequency NLES can be directly calculated from the LES outputs by

where the angle brackets represent the temporal average over the turbulent stage.

In this study, the LES-provided ϵLES, κLES, τLES, NLES, ESPLES, EϵLES and EBLES will be used as the “true” values, based on which the new parameterization would be built based on the initial variables like h0, N0¯, S0¯ and Ri0¯ (Figure 3b).

ESPLES comes from the release of MKE, which is approximately equal to the difference in kinetic energy between the initial unstable flow and the quiescent flow after turbulence. In Figures 1b, c, at the end of the turbulent stage, the strength of shear within the TL is reduced, and the mean Ri is close to about 0.25 at the center and boundary of the TL. This indicated that the flow now reaches a marginally stable state (Thorpe and Liu, 2009).

Based on this feature, Kunze et al. (1990) assumed that the shear in an unstable stratified shear flow would be reduced if turbulent fluxes raised the Ri = Ri0 < 0.25 to Ri = 0.25. Assuming that shear and stratification in the IUL are constant, the difference in kinetic energy between the initial unstable state and the final state of marginal instability is defined as the AKE ( Ka), and is calculated as Ka = h02(S2 – 4N2)/24. For nonlinear shear profiles like ours, Ka can be calculated as h02S02 – 4N02¯/24, where overbar represents the vertical average over the IUL. Though this computational approach scarifies certain physical fidelity compared to the rigorous numerical integration method, their values are consistent to a large extent (not shown); therefore, for calculation efficiency, we adopt this convenient method in the present study.

However, it is noted that the AKE ( Ka) is different from ESPLES. ESPLES represents the real released MKE throughout the TL turbulent processes, which is calculated from the LES results via (Equation 4) and is shown in Figure 3a. Ka is proposed for the purpose of parameterization, which represents the idealized amount of MKE released through the instability, without considering the complex energetics. For this reason, let ESPLES be expressed as

where the parameter λ1 is introduced in detail in the following framework of parameterization construction.

Furthermore, another parameter λ2 = EϵLES/ESPLES is introduced so that. is parameterized. Similarly, a third parameter λ3 = –EBLES/ESPLES is introduced. Under an assumption that the input energy is either transferred to potential energy or internal energy over the whole turbulent stage, λ3 is naturally equal to 1 – λ2 (Figure 3a). Then, the mean EϵLES and EBLES can be expressed by considering the turbulent timescale τLES as:

Here, λ1, λ2 and λ3 are variables to be parameterized by Ri0¯.

In sum, for each experiment, Ka is calculated by initial variables, while the energy transferred via shear production ( ESPLES) is calculated from LES outputs. By equating the parameterized ESPLES to the Ka (Equation 8), λ1 is derived for each experiment. Through regression analysis, λ1 is parameterized as a function of Ri0¯. Similarly, parameterizations for λ2, λ3, NLES and τLES are obtained. Ultimately, parameterized expressions for the turbulent diffusivity and dissipation rate are derived (Figure 3b).

A dataset of observations is collected to verify our parameterization. First, turbulence activity was measured during the Tropical Instability Wave Experiment (TIWE) in the fall of 1991 at 0°, 140°W (Lien et al., 1995). During this experiment, two overlapping time series of measurements were made from two independent ships, Wecoma and Moana Wave, so the method and data can be compared and validated. 3918 casts and 2072 casts of microstructure temperature, conductivity, and shear measurements in the upper 200 m were made using the profiler CHAMELEON and the Advanced Microstructure Profiler (AMP). The horizontal velocity was measured by the ship-mounted Acoustic Doppler Current Profilers (ADCPs) with the vertical resolution of approximately 8 m.

The TKE dissipation rate ϵobs is estimated by the method of sensing small-scale shears from the free-falling profilers (Moum et al., 1995). κobs is calculated as γϵobs/N2, where γ is taken as a common value of 0.2 (Lien et al., 1995; Zaron and Moum, 2009). Because of the occasional necessity of repairs and delays caused by other operational difficulties, the time series of profiles was unevenly sampled. To simplify the calculation, all data were averaged hourly with the vertical resolution of 1 m. In the next subsections, we will further process this dataset to comply with our EPP scheme.

Researchers usually directly apply a parameterization to observed hydrologic data to evaluate its performance. However, we note that it is difficult to fairly evaluate the performance this way. Firstly, observational data often lack the precise background variables that are required to initialize the potential turbulent events, unlike the well-controlled initial conditions in LES experiments. The so-called background fields may also have undergone the influence of prior turbulence. Secondly, the observed mixing coefficients (such as ϵ and κ) are subject to other larger-scale forces, such as advection and shear production, which is also unlike the freely developed turbulence as seen in LES. Lastly, turbulence observed at an observational site may originate from remote locations rather than local instability.

To compare EPP with observations, some turbulent events are picked out (marked by the white square in Figures 4a, b). As shown in Figures 4c, d, such turbulent events resemble LES experiments. Enhanced turbulence follows a fluid state within an IUL with Ri ∈ (0, 0.25).

The turbulent kinetic energy (TKE) and its dissipation rate, ϵ, are important metrics describing the development and decay of turbulence. The temporal variability of horizontally averaged TKE and ϵ over each model layer, and the identified TL and turbulent stage for experiment A₇B₇ are shown in Figure 2. τLES is about 4 hours. TKE increases rapidly in the domain of IUL during the onset of the turbulent stage. TKE surges to a peak (~ 5 × 10⁻³ m² s⁻²) within about 1 hour in the domain of IUL. Then TKE declines gradually to the background value over 3 hours. ϵ increases and maintains a value ranging 4 × 10⁻⁹ – 4 × 10⁻⁸ m² s⁻³ during the first 2 hours of the turbulent stage. After this, ϵ surges to more than 1.6 × 10⁻⁷ m² s⁻³ and keeps for about 2 hours. Although a short period with high ϵ is excluded from the turbulent stage, most of the turbulent characteristics have been captured. ϵ for the 27 experiments show a similar evolution as described above.

In addition, due to the vertical penetration of the turbulence, the TL becomes thicker rapidly. The thickness of the TL at the end of the turbulent stage is defined as the turbulent penetration thickness (hereafter TPT, denoted as H). H (90 m) is approximately 4 times larger than h0 (22 m) in experiment A₇B₇ (Figure 2). Strong vertical turbulent momentum and buoyancy fluxes therein result in a decrease in the temperature and velocity above the z0 and an increase of them below the z0 (Figure 1a). This variation has been shown in the direct numerical simulation results of Smyth and Winters (2003) and the CROCO ocean model results of Penney et al. (2020).

However, there is a significant difference in τLES across the 27 experiments, owing to the varied initial conditions (Figures 5a, b). τLES varies between 7920 s and 38880 s, and increases with increasing N0−1 which is often quoted for turbulence generation and dissipation (e.g., Moum, 1996b). Although τLES (i.e., 13680 s) for experiment A₇B₇ is much longer than the timescale N0−1 (i.e., 250 s), it is comparable with the value of Smyth et al. (2005). When N0 ≈ 1 × 10⁻² s⁻¹ (their Figure 1), their duration is about 5000 s and 50 times N0−1. It can also be found that τLES appear to increase with increasing S0−1, which is consistent with the results of Watanabe and Nagata (2021).

The TPTs of 27 experiments range from 20 m to 120 m, and their variations are significantly correlated with the Ri0¯ (Figure 5c). A linear regression of TPT on Ri0¯ can explain 96% of the variance. However, the IULT is also related to Ri0¯, thus the ratio of TPT to IULT, η, is a good index representing the penetration intensity of turbulence. In Figure 5d, η is a monotone-decreasing function of Ri0¯. TPTs can reach nearly 5 times the IULTs when the Ri0¯ is about 0.1, which is comparable to the values of 2–3 in Smyth et al. (2005) and 5 in Penney et al. (2020). Different from the theoretical results of Kunze (2014), turbulent entrainment can also occur even for Ri ∈ (0.03, 0.25).

To parameterize ESP with Ka, we calculated λ1 of each experiment, and found that λ1 varies between 1.5 and 25 (Figure 6a). When the velocity after simulated turbulent mixing is close to the prescribed idealized velocity of marginal instability, λ1 is small and close to 1; however, the more they differ from each other, the more MKE is released and the value of λ1 is larger than 1.

Calculation based on the LES results reveals that λ2 varies between 0.38 and 0.69 (Figure 6b). The larger λ2 is, the larger proportion of TKE is dissipated into the internal energy. It is noted that λ3 is intrinsically the flux Richardson number Rf, and λ3/λ2 is another measure of mixing efficiency γ which is the ratio of –EB to Eϵ (Smyth et al., 2001; Inoue and Smyth, 2009). We found that the values of γ range from about 0.4 to 1.4, which are larger than the commonly used value of 0.2; the underlying reason is that the calculation of γ includes the development stage of turbulence where γ is believed large (Smyth, 2020). Actually, γ can be larger than 1 when the flux Richardson number R_f is large, as seen in many numerical simulations and oceanic measurements (Moum, 1996a; Smyth et al., 2001; Salehipour and Peltier, 2015).

Finally, to parameterize ESP, Eϵ and EB with Ka, λ1 and λ2 should be parameterized at first. We found that λ1 and λ2 are good functions of Ri0¯ (Figures 6a, b). λ1 and λ2 both increase with increasing Ri0¯. λ1 increases exponentially with Ri0¯ while λ2 increases almost linearly with Ri0¯. Based on the data, we give the following fitting functions,

The regressions are shown in Figures 6a, b, respectively. Coefficient values and confidence intervals are listed in Table 2, and the residuals are presented in the Supplementary Material. λ1P can explain about 98% of the variance of λ1 (R² = 0.98), while λ2P can explain about 47% of the variance of λ2 (R² = 0.47). The low R² value indicates that nonlinear processes (e.g., vortex pairing) significantly modulate λ2. Nevertheless, Ri0¯ remains an important control parameter.

In Equation 9.1, timescale τLES need to be parameterized. Previously, σ =(S – 2N)/4 was used by Kunze et al. (1990) to estimate the inverse timescale for the growth of small amplitude billows based on linear stability analysis (Hazel, 1972), while N2/S was used to estimate the inverse timescale for the dissipation stage after the inception of shear turbulence based on the laboratory data (Thorpe, 1973). Polzin (1996) indicated that durations for turbulent events of observations during the North Atlantic Tracer Release Experiment encompass both growth and dissipation timescales. Considering that τLES is the full duration of turbulent evolution, which includes both the growth and decay stages based on nonlinear numerical simulations, we parameterize the turbulent timescale τLES as a linear combination of the two mentioned timescales,

where f and g are determined to be 0.04044 and 0.02861 by two-variable linear regression. The parameterization of τLES, i.e., τP, explains about 50% of the variance of τLES (R² = 0.50) as shown in Figure 6c. This expression is simple and easy to be used in the parameterization scheme. Accordingly, Equation 9.1 becomes:

Equation 12 is the energy-constrained parameterization for the TKE dissipation rate induced by the KH instability for the IUL where Ri ∈ (0, 0.25), which is represented by the original background variables.

In the meanwhile, it is found that the stratification NLES in Equation 9.2 has a significant linear relationship with the initial value N0¯, i.e., NP  = hN0¯ (Figure 6d). where h is determined to be 0.6761. NP explains about 82% of the variance of NLES (R² = 0.82). In addition, λ3 is parameterized as 1 – λ2P. Accordingly, Equation 9.2 becomes:

Equation 13 is the energy-constrained parameterization for κ of the KH instability-induced turbulence for the IUL where Ri ∈ (0, 0.25), which is represented by the original background variables.

The last key property of the present parameterization lies in the vertical extension of the turbulent mixing. An important information obtained from both previous studies and the present LES results is that the shear instability-driven turbulence is not confined within the IUL, but extends to the neighboring layers. This phenomenon represents the release of accumulated energy from a potentially unstable fluid system. The TPT should represent the outer boundary of the system where the energy can be extracted. Thinking from this way, it is necessary to parameterize κ within the TPT, rather than at the grid points of Ri ∈ (0, 0.25) only like in the previous parameterizations.

The suitable parameterization for this issue includes two steps. The first step is to identify TPT and represent it with original background variables, while the second step is to redistribute κ vertically within TPT.

Firstly, as described in section 3a, the ratio of TPT to IULT, η, varies from 2 to 5; furthermore, it can be fitted as a linear function of Ri0¯ (Figure 5d)

where i and j are determined as –19.61 and 6.58. ηP explains about 82% of the variance of η (R² = 0.82). In practice, IULT can be determined by the grid spacing of the oceanic numerical model. Given that η can be obtained by Ri0¯, the TPT can also be obtained easily, providing the layers where the parameterization should be exerted.

Secondly, the vertical distribution of κ should be provided. The previous discussion and parameterization mainly focused on the mean κLES averaged vertically over the TPT under the assumption of 3D homogeneity of turbulence (Kaltenbach et al., 1994; Shih et al., 2005). Whereas, κ can be also calculated layer by layer, which can provide the vertical pattern of κ within the TL. Specifically, for each experiment, the κ profile is calculated first, and then it is normalized by κLES to get κ*, in the meanwhile, z is normalized by TPT to get z*. Finally, the normalized profiles for 27 experiments are averaged, which is shown in Figure 7. It is found that κ* reaches its maximum value which is about 2 times the vertical average at the center of TPT where Ri0 is the minimum. Out of the deeper and shallower boundaries of TPT, κ* rapidly decreases to almost 0. Observed diffusivity profiles resulting from KH billow breakdown in the Changjiang Estuary closely match this vertical distribution (Tu et al., 2024). This pattern of normalized profile can be described by the fitting function:

where z* represents the normalized depth of TPT. The fitted profile is very close to the actual mean profile. The profile of κ is obtained by multiplying κEPP in Equation 13 by κ*. Till now, we have finished building the new parameterization for the shear instability-driven vertical mixing in the interior ocean (Equations 13, 15).

Overall, given the new parameterization is based on an energy-constrained framework and provides the vertical diffusivity profile, it is named the energy-constrained profile parameterization, and abbreviated as EPP.

It should be noted that, although the EPP schemes are based on the same set of data as ϵLES and κLES, thus are non-independent, they are constructed according to the theoretical framework of energy constraint, rather than by simply fitting to ϵLES and κLES. Therefore, ϵLES and κLES of LES can be used to test our scheme.

To evaluate the EPP scheme, the parameterized ϵEPP and the original ϵLES calculated from LES, are shown in Figure 8a. The parameterized values compare remarkably well to the values of LES. To be specific, ϵEPP explain about 81% of the variance of ϵLES. 96% of the samples show a discrepancy within a factor of 2 for ϵEPP, while about 70% of the samples show a discrepancy only within a factor of 1.5. The parameterized diffusivity κEPP and the LES-calculated diffusivity κLES are compared in Figure 8b. κEPP explains about 88% of the variance of κLES. 96% of the parameterized κEPP are within a factor of 1.5 to κLES.

Overall, the parameterized coefficients in the EPP scheme are in good agreement with the data calculated by LES, both in magnitude and variability.

The EPP scheme is also tested with independent observational data collected from the TIWE experiments (Lien et al., 1995). As described in section 2.6, 33 turbulent events similar to category 1 below the boundary layer (ML and deep cycle layer) are identified for the TIWE data (white squares in Figure 4).

When applying EPP to the IUL with Ri ∈ (0, 0.25) that is below the boundary layer in observations, the TKE dissipation rate and diffusivity are calculated according to the schematic in Figure 3b. Firstly, the temperature, salinity, and velocity u0 within the IUL can be used to calculate the initial variables, such as N0¯, S0¯ and Ri0¯. Next, λ1P, λ2P, λ3P, Ka, τP, NP and ηP can be obtained with Equations 10, 11, 14. Using these variables and Equations 12, 13, the vertical average ϵEPP and κEPP are calculated. Secondly, h0 is regarded as the sum of grid spacings of IUL. Multiplying h0 by the ratio ηP, the TPT and grids within which the turbulence can penetrate are obtained. Finally, by multiplying κEPP by the normalized profile from Equation 15, the κEPP profile is obtained. It is worth noting that, possibly due to the influence of other forcings, the stratification within turbulent events is not always smaller than the initial value N0¯ as in LES. Therefore, we use the observed buoyancy frequency as NP.

Applying the EPP to the 33 turbulent events, it is seen that the parameterized values can basically capture the magnitude and amplitude of the observed ϵobs and κobs (Figure 9). The agreements between ϵEPP and ϵobs, and between κEPP and κobs are both within a factor of 10 for about 88% and a factor of 5 for about 70%, respectively. Nonetheless, the EPP shows advantages compared to the widely used previous schemes, which is discussed in the next subsection.

Several previous parameterizations, including RSP, ZM (Zaron and Moum, 2009), PP81, P88 and KPP are compared, which shows an overall better performance of EPP. REV parameterization is used for ZM while shear instability mixing component of KPP is adopted for comparison. Through the comparison, we also analyze the underlying mechanisms why EPP performs better.

As shown in Figure 8a, if taking the LES as metrics, ϵRSP overestimates ϵLES, and has a larger scatter than ϵEPP, while only 59.3% of κRSP are within a factor of 2, compared to the 100% of ϵEPP. This may be due to the large inverse timescale employed in the RSP scheme, which results in a large ϵRSP, but its smaller γ partially offsets this overestimation. κKPP, κPP81 and κZM all underestimate the κLES by a factor of 4–18 on average, while κP88 overestimates them by a factor of 18. KPP, PP81, and ZM all prescribe distinct yet nearly invariant diffusivity when the Ri< 0.25 based on observed averages. Their empirical rigidity neglecting turbulence-scale dynamics in idealized LES experiments. For P88, the overestimation stems from its mathematical formulation where κP88 tends toward infinity as Ri approaches 0. The EPP seems to best fit the data, this is because this scheme fully considers the dynamical variables in addition to Ri0¯.

As for the parameterizations of the observations (Figure 9a), more of ϵRSP underestimates ϵobs than ϵEPP, but about 72% of ϵRSP approximate ϵobs well. 73% of the parameterized κRSP are within a factor of 10, while 58% of the samples are within a factor of 5. Compared with ϵEPP, more of κRSP underestimate κobs. κZM, κKPP, κPP81 and κP88 overestimate κobs by more than a factor of 10 on average, and they fail to capture the variability of κobs. The fact that variability of turbulent diffusivity depends not only on Ri0 but also on other variables such as shear and stratification (Richards et al., 2021) is obviously missed in these parametrizations. κKPP and κPP81 overestimate κobs, indicating that KPP and PP81 must adjust their parameters according to the observations at different conditions to achieve the best parameterization, but it is almost impossible to experience all different conditions. In contrast, the EPP and RSP is relatively adaptive to the observations and LES data.