Roles of Shear and Convection in Driving Mixing in the Ocean

Using field, numerical, and laboratory studies, we consider the roles of both shear and convection in driving mixing in the interior of the density‐stratified ocean. Shear mixing dominates when the Richardson number Ri < 0.25, convective mixing dominates when Ri > 1.0, and in the intermediate regime when 0.25 < Ri < 1.0 both shear and convection can contribute to mixing. For pure shear mixing the mixing efficiency Rif approaches 0.5, while for pure convective mixing the mixing efficiency Rif approaches 0.75. The diapycnal diffusivities for the two mechanisms are given by very different expressions. Despite these complexities, a simple mixing length model using the mean flow shear S provides robust estimates of diffusivity across the range 0 < Ri < 2. To account for the roles of both shear and convection over this range of Ri, we also formulate a modified version of the empirical KPP model for parameterizing ocean mixing in numerical models.

/ 0.25 Ri N S . Numerous studies have since extended this analysis and shown that, as long as Ri is not too large, shear flows are unstable when there is a viscous-diffusive fluid (Howland et al., 2018), when there is preexisting ambient turbulence in the fluid (Kaminski & Smyth, 2019), or when the shear is unsteady (Radko, 2019).
When Ri < 0.15 Venayagamoorthy and Koseff (2016) where B is the buoyancy flux. For Ri > 0.15 the estimates on the left and right hand sides of Equation 2 progressively diverge (see Venayagamoorthy & Koseff, 2016 their Figures 1 and 3 for details). However if Ri < 0.15, then the true diffusivity can be estimated from the buoyancy flux B, and hence Equation 1 can be simplified to the expression first proposed by Osborn (1980) where the flux Richardson number Ri f = B/(B + ϵ). The challenge in using Equation 3 is that Ri f must be independently known, and in practice Ri f is not constant (e.g., Bluteau et al., 2013 10.1029/2020GL089455 2 of 12 (Caulfield & Peltier, 2000) is the local mixing efficiency based on the re-sorted background density and, as seen from the last term in Equation 1, ℳ represents the irreversible diffusive flux.
While the model in Equation 1 is mechanism agnostic-there is no assumption about the mechanism causing mixing-it requires the ability to obtain a reference state of minimum potential energy by re-sorting the 3D density field. This can be accomplished either by re-sorting the 3D density field at a particular time step (as in a mumerical model), or by using a laboratory experiment with a closed control volume where (after mixing has occurred) the fluid adiabatically re-sorts itself as it comes to rest. As all ocean turbulence measurements are obtained from either a microstructure vertical profile or from moored time series at one location, the model in Equation 1 cannot be used with these turbulence measurements to compute the diffusivity.

Shear Mixing
The most commonly studied turbulent mixing mechanism in the stratified ocean is that driven by a vertically sheared mean flow. Miles (1961) and Howard (1961)  . Numerous studies have since extended this analysis and shown that, as long as Ri is not too large, shear flows are unstable when there is a viscous-diffusive fluid (Howland et al., 2018), when there is preexisting ambient turbulence in the fluid (Kaminski & Smyth, 2019), or when the shear is unsteady (Radko, 2019).

Geophysical Research Letters
While preceding work has demonstrated that both shear and convection mechanisms can drive diapycnal mixing in the interior of the ocean, a number of questions remain unresolved: (i) when does either mechanism dominate; (ii) what is the efficiency of mixing; (iii) what is the diapycnal diffusivity for each mechanism; (iv) can the two different mechanisms operate together; and (v) what are the implications for parameterizing ocean mixing? These questions are the focus of the present study.

Mixing Mechanisms
To quantify mixing in a density-stratified fluid, Winters et al. (1995) argued one must distinguish between the reversible and irreversible components of mixing. Using this framework,  generalized a model for the (true) diapycnal diffusivity 

K Ri
Ri where the (true) flux Richardson number Caulfield & Peltier, 2000) is the local mixing efficiency based on the re-sorted background density and, as seen from the last term in Equation 1, ℳ represents the irreversible diffusive flux.
While the model in Equation 1 is mechanism agnostic-there is no assumption about the mechanism causing mixing-it requires the ability to obtain a reference state of minimum potential energy by re-sorting the 3D density field. This can be accomplished either by re-sorting the 3D density field at a particular time step (as in a mumerical model), or by using a laboratory experiment with a closed control volume where (after mixing has occurred) the fluid adiabatically re-sorts itself as it comes to rest. As all ocean turbulence measurements are obtained from either a microstructure vertical profile or from moored time series at one location, the model in Equation 1 cannot be used with these turbulence measurements to compute the diffusivity.

Shear Mixing
The most commonly studied turbulent mixing mechanism in the stratified ocean is that driven by a vertically sheared mean flow. Miles (1961) and Howard (1961)  . Numerous studies have since extended this analysis and shown that, as long as Ri is not too large, shear flows are unstable when there is a viscous-diffusive fluid (Howland et al., 2018), when there is preexisting ambient turbulence in the fluid (Kaminski & Smyth, 2019), or when the shear is unsteady (Radko, 2019).
When Ri < 0.15 Venayagamoorthy and Koseff (2016) where the flux Richardson number Ri f = B/(B + ϵ). The challenge in using Equation 3 is that Ri f must be independently known, and in practice Ri f is not constant (e.g., Bluteau et al., 2013 where the flux Richardson number Ri f = B/(B + ϵ). The challenge in using Equation 3 is that Ri f must be independently known, and in practice Ri f is not constant (e.g., Bluteau et al., 2013;Davis & Monismith, 2011;Holleman et al., 2016;Ivey et al., 2018;Ijichi & Hibiya, 2018;Shih et al., 2005;Walter et al., 2014). And it has been argued that Ri f cannot even be expressed as a function of a single dimensionless parameter (e.g., Mater & Venayagamoorthy, 2014b).
This problem of needing to first specify Ri f is circumvented with the model of Osborn and Cox (1972) who proposed that the vertical diffusivity for heat K θ can be estimated as If the temperature gradient    / z dominates the density gradient, then one can assume that K θ is equal to K ρ . Direct numerical simulation (DNS) studies of stratified shear flows have compared K θ with K ρ (e.g. Itsweire et al., 1993), and recent simulations have shown that K θ is within a factor or 2 (either low or high) of  * K for buoyancy Reynolds numbers Re b = ϵ/(νN 2 ) as large as 5 × 10 4 (Kirkpatrick et al., 2019;. Although restricted to cases with Re b < 250, the DNS study of unsheared stratified mixing by Taylor et al. (2019) also found that Equation 4 was a good estimator of  * K .
A third mixing model, proposed by Odier et al. (2009), is based on a Prandtl mixing length model for a shear flow and is given by When temperature dominates the density gradient, the mixing length scale L ρ is (Ivey et al., 2018)  And with L ρ now defined, the mixing efficiency Ri f in Equation 3 can be written as (Ivey et al., 2018) where   Using field and laboratory observations obtained across a range of length scales spanning three orders of magnitude, Ivey et al. (2018), also showed the mixing length scale L ρ ≈ 0.3L E , where the Ellison length scale The Ellison scale can be computed from time series measurements and has the advantage of not requiring direct estimates of turbulence quantities such as χ or ϵ from specialized instruments. The Ellison scale does, however, require sampling fast enough to resolve the larger energy containing scales and records that are not influenced by internal waves or mean flows. Equation 5 can thus be approximated by Almost all numerical studies have focused on shear flows when Ri < 0.25 and are unstable with active mixing (e.g., Shih et al., 2005;Smyth et al., 2001). But what happens when Ri > 0.25? As we show below, mixing still occurs in both field observations and numerical modeling and we evaluate the predictive capability of the shear-driven mixing models in Equations 3, 4, and 8 when Ri > 0.25. We also determine when the mechanism of convective instability starts to become important.

Convective Mixing
By convective mixing we refer here to mixing that occurs in the interior of the density-stratified water column due to a local instability in the mean density gradient. This process is distinct from mixing driven by an externally imposed buoyancy flux at the free surface (as considered by Sohail et al., 2018, for example, where processes such as surface cooling, evaporation or salt rejection can occur). A common feature of interior convective mixing observed in either laboratory experiments (Wykes & Dalziel, 2014) or in idealized numerical simulations (Chalamalla & Sarkar, 2015;Gayen & Sarkar, 2011;Puthan et al., 2019) is the often efficient mixing with  * 0.5 f Ri . As discussed above, as internal waves always induce some background vertical velocity shear S, internal wave forcing does not drive pure convective instabilities in the ocean.
A good example of pure convective instability, however, is the laboratory experiments of Barry et al. (2001). They conducted a zero mean shear mixing experiment by horizontally oscillating a vertically oriented square-bar mesh grid in a rectangular tank filled with a density stratified fluid (see their Figure 1). The tank was filled with a linearly salt-stratified fluid, the initial density profile was measured, and the fluid was then continuously agitated by rapidly moving the grid back and forth throughout the entire tank (but without creating any mean flow). Once the grid was stopped, the fluid was allowed to come to rest and the final density profile was then measured. The diffusivity  * K was calculated from the increase in background potential energy after mixing, and N * was taken as the average of the two (near-linear) before and after density profiles. The moving grid was instrumented with a force transducer and the dissipation ϵ in the fluid was calculated directly. The experiments varied ϵ, N * , and fluid viscosity ν but, as S = 0, the experiments all had The buoyancy Reynolds number reached values as large as Re b = 10 5 , comparable to the highest values seen in the ocean (e.g., Gargett et al., 1984). Barry et al. (2001) found from their measurements that the diapycnal diffusivity was given by is independent of the large scale shear S which, as noted by Mater and Venayagamoorthy (2014a), is a key difference between pure convective and shear mixing.
While Barry et al. (2001) did not discuss or provide estimates of mixing efficiency * f Ri , these can now be calculated from their measurements using the framework proposed by . In particular, substituting Equation 9 into Equation 1 yields

Geophysical Research Letters
the predictive capability of the shear-driven mixing models in Equations 3, 4, and 8 when Ri > 0.25. We also determine when the mechanism of convective instability starts to become important.

Convective Mixing
By convective mixing we refer here to mixing that occurs in the interior of the density-stratified water column due to a local instability in the mean density gradient. This process is distinct from mixing driven by an externally imposed buoyancy flux at the free surface (as considered by Sohail et al., 2018, for example, where processes such as surface cooling, evaporation or salt rejection can occur). A common feature of interior convective mixing observed in either laboratory experiments (Wykes & Dalziel, 2014) or in idealized numerical simulations (Chalamalla & Sarkar, 2015;Gayen & Sarkar, 2011;Puthan et al., 2019) is the often efficient mixing with  * 0.5 f Ri . As discussed above, as internal waves always induce some background vertical velocity shear S, internal wave forcing does not drive pure convective instabilities in the ocean.
A good example of pure convective instability, however, is the laboratory experiments of Barry et al. (2001). They conducted a zero mean shear mixing experiment by horizontally oscillating a vertically oriented square-bar mesh grid in a rectangular tank filled with a density stratified fluid (see their Figure 1). The tank was filled with a linearly salt-stratified fluid, the initial density profile was measured, and the fluid was then continuously agitated by rapidly moving the grid back and forth throughout the entire tank (but without creating any mean flow). Once the grid was stopped, the fluid was allowed to come to rest and the final density profile was then measured. The diffusivity  * K was calculated from the increase in background potential energy after mixing, and N * was taken as the average of the two (near-linear) before and after density profiles. The moving grid was instrumented with a force transducer and the dissipation ϵ in the fluid was calculated directly. The experiments varied ϵ, N * , and fluid viscosity ν but, as S = 0, the experiments all had The buoyancy Reynolds number reached values as large as Re b = 10 5 , comparable to the highest values seen in the ocean (e.g., Gargett et al., 1984). Barry et al. (2001) found from their measurements that the diapycnal diffusivity was given by is independent of the large scale shear S which, as noted by Mater and Venayagamoorthy (2014a), is a key difference between pure convective and shear mixing.
While Barry et al. (2001) did not discuss or provide estimates of mixing efficiency * f Ri , these can now be calculated from their measurements using the framework proposed by . In particular, substituting Equation 9 into Equation 1 yields For a fixed Pr, Equation 10 demonstrates that Ri L L , a very different length scale dependence than that found in Equation 7 for a shear flow.
Using Equations 9 and 10, we can estimate the range of , respectively. Equations 9 and 10 thus predict wide ranges in mixing efficiencies and diffusivities resulting from convective mixing in the ocean.

Mixing Mechanisms in the Ocean
The shear and convective mixing mechanisms discussed above are thus quite different, as shown by the differing expressions for both the mixing efficiencies in Equation 7 versus Equation 10 and the diffusivities in Equation 3 versus Equation 9. The diffusivity in Equation 9 is derived from a laboratory experiment with zero mean shear and hence Ri is infinite, but in the ocean there is inevitably some shear and so Ri will almost always be finite. Are the diffusivities observed in the ocean therefore ever consistent with the prediction for pure convective mixing given by Equation 9 shear mixing dominate ocean mixing? Can the two mechanisms operate together? To address these questions, we examine below two sources of information: field observations collected on the Australian North West Shelf (NWS) and numerical simulations of an oscillating stratified flow over a smooth bottom.

Field Observations
The field observations come from a 34 m long mooring bottom-anchored on the NWS in 105 m of water. The site is energized by strong and nonlinear internal tidal motions on the NWS. Details of the measurements are given in Bluteau et al. (2016a), but briefly the mooring was deployed from 4-22 April 2012 and was instrumented with: (i) 22 temperature sensors (including Seabird Electronics SBE56 and SBE39 sensors) sampling respectively at 0.5 and 10 s; (ii) a conductivity-temperature sensor (Seabird Electronics, SBE37) sampling at 15 s; (iii) a temperature-pressure sensor sampling at 10 s (SBE39, Sea-Bird Electronics); and (iv) an upward-looking 300 kHz acoustic-Doppler current profiler (ADCP, RDInstruments) measuring velocities in 1 m bins over 1-min averages (see Table 1  Estimates with L E > L z were removed from further analysis, leaving us with observations of the density stratified fluid at a depth of around 100 m, but free from any surface or bottom influence. Each MTP included a motion pack to recover the instruments' motion and correct the velocity spectral observations using the techniques described by Bluteau et al. (2016b). We fitted the ADV and FP07 spectral observations (computed with 8.53 min segments) over their respective inertial subranges to yield estimates ϵ and hence χ using the inertial subrange fitting methods described by Bluteau et al. (2011Bluteau et al. ( , 2017. With χ determined, the diffusivity K θ was then determined from Equation 4. Whether the mixing mechanism is shear or convective in nature, for the reasons discussed above in Section 2.1 we assume that K θ is a good estimator of K ρ * .

Numerical Simulations
Using We chose a computational grid with an effective grid resolution in the mean flow direction of x + = (xu * /ν) = 60 (where u * was the bottom friction velocity) and in the cross flow direction of y + = (yu * /ν) = 30. The horizontal grid resolution in both the mean and transverse flow directions was uniform throughout the domain. In the vertical direction, a variable grid resolution was used with strong clustering near the bottom with z + = (zu * /ν) = 2 to resolve the viscous sublayer. This resolution was sufficient to resolve the near-wall eddies that carry the Reynolds stress. Far from the bottom in the weakly sheared interior the vertical grid scale was relaxed to z + = (zu * /ν) = 20. As discussed in Gayen et al. (2010), our numerical simulation thus resolved the flow near the bottom and did not require any near-wall closure model. A dynamic eddy viscosity and eddy diffusivity model was used for subgrid scale closure well away from the bottom in the interior.
After just one tidal cycle, the flow started to form a three-layer structure (see Figure 1a of Gayen et al., 2010). The three-layer structure consisted of a strongly sheared (but nearly well-mixed layer) bottom layer of height h m , a second weakly sheared (but linearly stratified) region in the interior, and a third relatively thin (and strongly sheared) pycnocline region separating the bottom layer from the linearly stratified interior. The height of the bottom mixing layer h m was defined as the height where the local density gradient first reached 0.1 of the interior density gradient. The gradient Richardson number Ri was generally very small near the bottom and increased with height. The height above the bottom where Ri first reached 0.25 was defined as h 0.25 (Gayen et al., 2010 found h m ≈ 0.9h 0.25 ), and the value of Ri then continued to increase with height above h 0.25 . A classical logarithmic mean velocity profile was only seen during a small portion of the tidal cycle. All mixing estimates were calculated once the simulation had reached a statistical steady state, which took about 10-12 tidal cycles. Over the next five tidal cycles, all mean and turbulence properties at a particular instant and height were horizontally averaged in the (x, y) plane, a spatial averaging process denoted by <>. Data were taken in the vertical over the range z + = 200 − 800. The Ellison scale was calculated as z and the vertical eddy diffusivity was calculated using two methods. First, to allow direct comparison with the analysis methods used in processing the field data above, we calculated the diffusivity from the buoyancy flux B as K ρ = <B / N 2 >. Second, we re-sorted the instantaneous density fields to get the irreversible diffusive flux ℳ and hence also calculated Geophysical Research Letters z + = (zu * /ν) = 2 to resolve the viscous sublayer. This resolution was sufficient to resolve the near-wall eddies that carry the Reynolds stress. Far from the bottom in the weakly sheared interior the vertical grid scale was relaxed to z + = (zu * /ν) = 20. As discussed in Gayen et al. (2010), our numerical simulation thus resolved the flow near the bottom and did not require any near-wall closure model. A dynamic eddy viscosity and eddy diffusivity model was used for subgrid scale closure well away from the bottom in the interior.
After just one tidal cycle, the flow started to form a three-layer structure (see Figure 1a of Gayen et al., 2010). The three-layer structure consisted of a strongly sheared (but nearly well-mixed layer) bottom layer of height h m , a second weakly sheared (but linearly stratified) region in the interior, and a third relatively thin (and strongly sheared) pycnocline region separating the bottom layer from the linearly stratified interior. The height of the bottom mixing layer h m was defined as the height where the local density gradient first reached 0.1 of the interior density gradient. The gradient Richardson number Ri was generally very small near the bottom and increased with height. The height above the bottom where Ri first reached 0.25 was defined as h 0.25 (Gayen et al., 2010 found h m ≈ 0.9h 0.25 ), and the value of Ri then continued to increase with height above h 0.25 . A classical logarithmic mean velocity profile was only seen during a small portion of the tidal cycle. All mixing estimates were calculated once the simulation had reached a statistical steady state, which took about 10-12 tidal cycles. Over the next five tidal cycles, all mean and turbulence properties at a particular instant and height were horizontally averaged in the (x, y) plane, a spatial averaging process denoted by <>. Data were taken in the vertical over the range z + = 200 − 800. The Ellison scale was calculated as z and the vertical eddy diffusivity was calculated using two methods. First, to allow direct comparison with the analysis methods used in processing the field data above, we calculated the diffusivity from the buoyancy flux B as K ρ = <B / N 2 >. Second, we re-sorted the instantaneous density fields to get the irreversible diffusive flux ℳ and hence also calculated  

Field Observations
The values of both Ri and K θ varied strongly during the 10 days over which the MTPs collected turbulence observations (Figure 1). Around 65% of the data had Ri < 0.25 and 95% of the data had Ri less than 0.6 (Figure 1b). Around 90% of the events had diffusivities K θ between 10 −4 m 2 s −1 and 10 −2 m 2 s −1 (Figure 1c). The largest Ri observed was Ri ≈ 1.25, but mixing events with such large values of Ri were rare. For Ri < 0.02 bin-averaged diffusivities K θ > 10 −2 m 2 s −1 and in this limit the observations had Re b > 10 5 ( Figure 1a). As Ri increased, the diffusivities then progressively decreased with increasing Ri (Figure 1a). Zaron and Moum (2009) suggested a modified version of the original KPP model of Large et al. (1994), where this decay was exponential with Ri with an exponent of −10. As shown in Figure 1a, this appears a reasonable description for small Ri < 0.15, but for most of the Ri range a lower value for the exponent is required, and the second line shown for comparison in Figure 1a) has an exponent of −2. This choice of exponents is motivated by the numerical results discussed below.
In Figure 2a we compare these observed diffusivities K θ with the predicted diffusivities for the mixing length model in Equation 8. While individual data points are scattered, once binned the observations are consistent with the predictions within the error bounds shown over the observed range of Ri. In Figure 2b we compare these observed diffusivities with the predictions for convective mixing in Equation 9. For Ri close to 0, Equation 9 clearly underpredicts the observed K θ , but as Ri increases the binned observations approach the predicted value for convectively dominated mixing. The challenge is there are very few estimates in the field data with Ri > 0.75, reflecting the practical difficulty of making reliable estimates of the mean velocity gradients to accurately estimate high Ri values in the field. This constraint with the available field data is one motivation for examining independent data from a numerical model, as discussed below.

Numerical Simulations
The results of the numerical simulations are shown in Figure 3. Overall around 85% of the data have

Field Observations
The values of both Ri and K θ varied strongly during the 10 days over which the MTPs collected turbulence observations ( Figure 1). Around 65% of the data had Ri < 0.25 and 95% of the data had Ri less than 0.6 (Figure 1b). Around 90% of the events had diffusivities K θ between 10 −4 m 2 s −1 and 10 −2 m 2 s −1 (Figure 1c). The largest Ri observed was Ri ≈ 1.25, but mixing events with such large values of Ri were rare. For Ri < 0.02 bin-averaged diffusivities K θ > 10 −2 m 2 s −1 and in this limit the observations had Re b > 10 5 ( Figure 1a). As Ri increased, the diffusivities then progressively decreased with increasing Ri (Figure 1a). Zaron and Moum (2009) suggested a modified version of the original KPP model of Large et al. (1994), where this decay was exponential with Ri with an exponent of −10. As shown in Figure 1a, this appears a reasonable description for small Ri < 0.15, but for most of the Ri range a lower value for the exponent is required, and the second line shown for comparison in Figure 1a) has an exponent of −2. This choice of exponents is motivated by the numerical results discussed below.
In Figure 2a we compare these observed diffusivities K θ with the predicted diffusivities for the mixing length model in Equation 8. While individual data points are scattered, once binned the observations are consistent with the predictions within the error bounds shown over the observed range of Ri. In Figure 2b we compare these observed diffusivities with the predictions for convective mixing in Equation 9. For Ri close to 0, Equation 9 clearly underpredicts the observed K θ , but as Ri increases the binned observations approach the predicted value for convectively dominated mixing. The challenge is there are very few estimates in the field data with Ri > 0.75, reflecting the practical difficulty of making reliable estimates of the mean velocity gradients to accurately estimate high Ri values in the field. This constraint with the available field data is one motivation for examining independent data from a numerical model, as discussed below.

Numerical Simulations
The results of the numerical simulations are shown in Figure 3. Overall around 85% of the data have 10.1029/2020GL089455 the oscillating flow was reversing and the vertical shear S was small. In Figure 3a, we plot K ρ versus Ri and in Figure 3b we plot  * K versus Ri. What is striking is that K ρ appears a very good representation of  * K . For very low Ri, both diffusivities approach 10 −4 m 2 s −1 . As Ri increases, both diffusivities decrease progressively toward 10 −6 m 2 s −1 as Ri approaches 2. Also shown in Figures 3a and 3b, are exponential decay lines for the diffusivities as a function of Ri. Consistent with the field data, an exponential decay with an exponent of −10 is reasonable only for small Ri. For 0.15 < Ri < 2 a much smaller exponent is clearly required, and for comparison we show a line with an exponent of −2. This is a simple first suggestion and clearly more data is required, particularly in the regime where Ri > 1.25. In summary, the field observations ( Figure 2a) and the numerical simulations (Figure 3b) thus show a very similar decrease of diffusivity with increasing Ri, although the numerical results extend to much higher Ri than the field observations.
Given the very similar dependence of K ρ and  * K with Ri in Figures 3a and 3b, we restrict our analysis to K ρ in the remaining panels. In Figures 3c and 3d, we compare the diffusivities K ρ from the simulations with those predicted from the same two mixing models evaluated above with the field observations. The predictions from the mixing length model in Equation 8 are in good agreement with the numerical diffusivities out to the observed upper bound of Ri = 2 (Figure 3c). Figure 3d compares the simulation diffusivities K ρ with the predicted diffusivities for convective mixing in Equation 9. For small Ri, the convective model underestimates the diffusivities by as much as a factor of 5. When Ri > 1, however, more than 80% of the diffusivities are within a factor of 0.5 (high or low) of the predictions of the convective model in Equation 9.
IVEY ET AL.

Discussion and Conclusions
The field observations and the numerical simulations are very consistent in demonstrating that the Richardson number Ri provides a useful differentiation between the shear and convective mixing mechanisms (from hereon we drop the superscript * for simplicity of notation). Our results show that when Ri < 0.25, the mixing is dominated by shear, the mixing efficiency Ri f = Ri f (L O , L S , L ρ ), and Ri f is in the range from 0 to 0.5. When Ri > 1.0, the mixing is dominated by convection, the mixing efficiency Ri f = Ri f (L O ,L K ), and Ri f is in the range from 0 to 0.75. In the intermediate case when 0.25 < Ri < 1.0, our results imply that both the shear and convective mechanisms can contribute to mixing.
The field and numerical results both show there is a progressive two order of magnitude decrease in K ρ as Ri increases from 0 to 2. However, the field and numerical results show that K ρ is not determined solely by Ri: at the same value of Ri, the field values of K ρ are nearly two orders of magnitude larger. For both the shear IVEY ET AL.
10.1029/2020GL089455 9 of 12 There is a factor of 10 2 difference in the Re b between the field observations in Figure 1 and the numerical results in Figure 3. The much higher value of Re b in the field seems the likely cause for the difference in the magnitude of K ρ for the same Ri.
Our analysis shows that while Ri is a useful parameter in distinguishing between shear and convective mixing, Ri alone cannot be used to predict diffusivity in the ocean. The popular KPP closure model of Large et al. (1994) makes four assumptions: mixing only occurs due to shear; K ρ depends only on Ri; in the limit when Ri → 0 the diffusivity K ρ → K 0 = 5 × 10 −3 m 2 s −1 ; and when Ri > 0.7 the diffusivity K ρ = 0. None of these assumptions are supported by our results.
Our results can, however, be used to formulate a modified version of the KPP model similar to that suggested by Large et al. (1994) (a decaying polynomial dependence on Ri) and Zaron and Moum (2009) (a decaying exponential dependence on Ri). Consistent with both the field and numerical results presented here, the following empirical expression accounts for both shear and convective mixing over a wide range of Ri where Re b0 is the buoyancy Reynolds number in the limit when Ri → 0. For the field data, Re b0 ≈ 10 5 and K 0 ≈ 10 −2 m 2 s −1 (Figure 1b)-only a factor of 2 above the K 0 value proposed by Large et al. (1994)-while for the numerical data Re b0 ≈ 10 3 and K 0 ≈ 8 × 10 −5 m 2 s −1 (Figure 3a). Equation 11 requires as input a value of Re b0 in order to determine K 0 . As a first step, we propose adopting the same K 0 = 5 × 10 −3 m 2 s −1 suggested by Large et al. (1994) and then using Equation 11 over the full range of the computed Ri. The exponent of -10 is reasonable only for Ri < 0.15 and a lower exponent is needed as Ri increases. There are few data points at large Ri and the exponent of −2 is just a first suggestion and clearly more data, particularly at Ri > 1.25, is required to refine this estimate. And more generally, our proposed formulation in Equation 11 needs to be tested further with measurements of K 0 and estimates of Re b0 before use in modeling applications.
The simple mixing length model in Equation 8 appears to be a reliable predictor of the observed diffusivities K ρ in both the field observations and the numerical results reported here. In practice, this agreement occurs over a broader range of Ri than might be anticipated from mixing length theory, and even appears to be a good predictor of mixing at high Ri when mixing is dominated by convection. Whether the mechanism driving the mixing is shear or convection and irrespective of the value of Ri, it appears Equation 8 can be applied to mooring data to get long time estimates of K ρ . Such long-term estimates of K ρ would be invaluable in evaluating the description of mixing within circulation models.