## Abstract

Osborn’s method is commonly used to obtain diffusion coefficient *K* from the turbulent dissipation rate *ε*. This method is the relational expression . The dissipation flux coefficient is often set to a constant value of 0.2, but this study of LES revealed that varies greatly vertically in the bottom boundary layer because of the influence of the seabed. Consequently, the eddy diffusion coefficient is overestimated in the lower part of the bottom boundary, but it is slightly underestimated in the upper part. Therefore, Osborn’s method with constant cannot give the correct diffusivity. Furthermore, even if treating as a function of flux Richardson number , as defined originally by Osborn, the estimation is underestimated by the advection effect because of the influence of spatial nonuniformity. Energy budget analysis revealed that this defect can be improved using the extended flux Richardson number, which can be estimated by multiplying using a constant correction factor. Furthermore, we proposed two alternative estimation methods. For the first method, which estimates the relation between and the gradient Richardson number , can be expressed with instead of with a correction factor. We can estimate the reasonable diffusivity if we have current data supplementary to obtain . For the second method, can be expressed as a similarity function of the height above the bottom normalized by the Ozmidov scale. This method can provide an acceptable estimate of diffusivity without current data for several circumstances.

## 1. Introduction

Vertical diffusion is an important process in ocean science. Diffusive processes attributable to turbulence can transport vast amounts of nutrients from nearby seabed areas to the photic layer, which increases oceanic primary production (Tilzer and Goldman 1978; Cloern 1991; Furevik and Foldvik 1996; Cross et al. 2014). Heat transport caused by vertical diffusivity is also well known to be related to maintenance of the global meridional overturning circulation and to climate change (Munk 1966; Munk and Wunsch 1998). Direct ship observations for vertical diffusivity are extremely difficult. Recently, however, we can accomplish them indirectly from the turbulent dissipation rate, which is obtained using a microstructure profiler (Oakey 1982; Moum and Osborn 1986; Hibiya et al. 2007) or an expendable current profiler (Nagasawa et al. 2002; Hibiya et al. 2006). Hibiya et al. (2006) mapped the global distribution of vertical eddy diffusivity. Reportedly, this technique has been used also in coastal oceans. Peters (1999), Rippeth et al. (2001), and Tsutsumi and Matsuno (2012) discussed the turbulent diffusivity related to turbulence near the seabed in shallow oceans. To calculate the eddy diffusivity (turbulent diffusion coefficient) from the dissipation rate, one must know the relational expression linking these physical quantities.

Osborn (1980) proposed a very useful estimation formula for vertical diffusivity from the dissipation rate. The many researchers described above used this formula to obtain the vertical eddy diffusivity. This formula is obtained from the turbulent kinetic energy (TKE) budget. If we can assume the statistical stationarity and uniformity for turbulence, then the sum of TKE production rate *P* (defined by the first term enclosed by parentheses in the next equation) and buoyancy flux *B* (the second term) is expected to be balanced with the energy dissipation rate *ε* as

where the quantities with an attached prime are the anomaly (turbulent) components from horizontal mean values, the bar denotes a horizontal mean, *g* represents the gravitational acceleration coefficient, *ρ* stands for density, and expresses a constant reference density. When obtaining this relation, the terms of time derivative of TKE and terms of space derivative related to the advection and the pressure transport expressed in the flux form were eliminated from the energy budget equation based on the assumptions of statistical stationarity and uniformity, respectively.

Vertical diffusivity is related to the buoyancy flux *B* because the vertical diffusivity can be expressed as

Therefore,

where is the square of Brunt–Väisälä frequency. Strictly speaking, the density in the definition of is the potential density. The difference for the selection of density or potential density is not large for shallow seas. The ratio of the buoyancy flux and production rate is known as the flux Richardson number, which represents the ratio for potential energy conversion against kinetic energy production. The sum energy produced by *P* and *B* is consumed by heat energy (internal energy of water). After substituting these relations into (1), the diffusivity can be expressed as

where coefficient is known as the dissipation flux coefficient (Moum 1996), turbulent flux coefficient (Zhou et al. 2017), and is sometimes known as mixing efficiency by oceanographers, which is well known as Osborn’s method. Knowing this coefficient (i.e., ), one can obtain the vertical diffusivity from the dissipation rate and Brunt–Väisälä frequency. Reportedly, has a critical value for the steady turbulent state (Ellison 1957). Osborn (1980) reported adopting 0.15 for the critical flux Richardson number, which suggests . However, many researchers have adopted this coefficient = 0.2 and an equality relation rather than the inequality relation as to obtain the mixing coefficient.

Recently, some researchers have questioned the use of as a constant value (Lozovatsky and Fernando 2013; Mater and Venayagamoorthy 2014b,a). Davis and Monismith (2011), Dunckley et al. (2012), and Walter et al. (2014) showed wide variation of from observations. Barry et al. (2001), Shih et al. (2005), and Rehmann (2004) proposed an alternative power law for the relation between the diffusivity *K* and turbulent activity parameter , where cannot be a constant value.

One objective of this paper is confirmation of the applicability of Osborn’s method for bottom boundary layer turbulence developed under a geostrophic current from a large-eddy simulation (LES) study. We shall investigate two questions: whether the spatial and temporal uniformity on the derivation process of Osborn’s equation is correct or not, and whether the (i.e., ) has a uniform value or not. For this purpose, we shall investigate the budget of TKE, demonstrating the inapplicability of Osborn’s method. We propose alternative estimation formulas for diffusivity. We will propose obtained by considering correction for spatial nonuniformity, by considering the relation of and (gradient Richardson number), and depending on the vertical distance *z* from the seabed scaled with the Ozmidov scale.

Section 2 presents a description of the model configuration. Section 3 presents numerical results and examines the validity of Osborn’s formula. The TKE budget is also analyzed. In section 4, we propose alternate estimation formulas for diffusivity. The last section is devoted to a summary and discussion.

## 2. Model configuration

We investigate the turbulence generated near the ocean floor overlying geostrophic flow in a stratified ocean. Vertical mixing of the density attributable to the turbulence will be specifically examined.

An LES model, a parallelized large-eddy simulation model (PALM), developed by Raasch and Schröter (2001) at Hannover University, was used for this study. The model configuration is almost identical to that described by Wakata et al. (2017). This model solves Navier–Stokes equations with Boussinesq approximation. The prognostic equations are separated into a grid-scale (GS) component and a subgrid-scale (SGS) component using a local spatial filter, which differs from the familiar turbulence models such as the model and Mellor–Yamada model using the ensemble mean. The SGS TKE was solved according to Deardorff (1980), with minor changes. The explicit description of the GS large-eddy turbulent motions can be given by the numerical simulation. The top boundary condition is assumed to be slip. The Prandtl boundary layer (constant momentum flux layer) is adopted between the bottom surface and the first computational grid level, reflecting the no-slip condition for the bottom boundary condition. The bottom roughness height is 0.01 m. This value becomes a few centimeters for a barred beach (Garcez Faria et al. 1998). The computational area is 250 m × 250 m × 103 m; the grid size is 1 m × 1 m × 0.5 m. The horizontal boundary condition is assumed to be periodic, which means that the west-side (north side) boundary variables are linked to eastward-side (south side) ones. The horizontal domain is extended periodically. The eastward velocity , which is vertically constant and in a geostrophic balance, is set as the initial condition. The northward body forcing representing the sea surface height gradient is set in the momentum equation to balance with the Coriolis force related to the initial eastward velocity. The GS momentum equations are described here; variables attached with a bar are GS variables:

In this equation, suffixes represent numbers 1–3, are Cartesian coordinates, and are the velocity components . The summation is taken as 1–3 hereafter if the same suffix appears twice in each term (Einstein convention). Parameter is the Coriolis parameter with Earth’s angular velocity and latitude *ϕ*, is the reference density of seawater, *g* is the gravitational acceleration parameter, is Kronecker’s delta, and is Eddington’s epsilon (Levi-Civita symbol). The last term of the right side is the external forcing term balancing with the Coriolis force related to the initial eastward velocity . Without this term, the fluid starts to rotate because of the inertial oscillation. Here is the modified pressure with perturbation pressure and SGS turbulent kinetic energy *e*. The term is a stress term related to SGS motion, which is obtained as . The coefficient parameter is obtained from a 1.5-order closure after Deardorff (1980). Details are presented in a report of a study by Maronga et al. (2015). The initial potential temperature is assumed to be linearly increasing upward as

where is a coefficient and represents the depth to the ocean floor. The GS potential temperature is solved from the following equation:

The first term shows the advection attributable to the GS current. The second term on the right side shows the influence of SGS as , where is determined from the SGS closure model. The parameters for the respective cases are presented in Table 1. The standard case (case A in Table 1) is explained here. The latitude is set as 35°N. For the initial condition, the salinity is set to 35 in practical salinity units, and °C m^{−1} and °C for the potential temperature. The state equation for seawater proposed by Jackett et al. (2006) is adopted. Because we assumed homogeneous salinity, the salinity does not change over time. The density is determined only from the change of temperature. The velocity is set to 1 m s^{−1}. The small random disturbances of velocity with amplitude of 0.01 m s^{−1} are added to for the initial velocity to become the seed of shear instability. The same values for salinity and are used for all cases studies. The simulation was performed for 100 000 s.

## 3. Simulation results

We introduce the property of the turbulence that develops near the seabed under the uniform geostrophic eastward flow in this section. We show the time development of turbulence, and the horizontal mean field. The vertical profiles of the turbulent transfer terms in the energy budget are analyzed to evaluate the applicability of the Osborn’s method near the seabed.

### a. Time evolution of turbulence near the seabed

We first explain how the turbulence develops with time. The vertical cross section of east–west velocity for case A is shown in Fig. 1. Strong vertical shear caused by the bottom stress appears near the seabed (5000 s). Unstable waves appear by the shear instability (6000 s). This unstable wave breaks abruptly. The turbulence develops gradually near the seabed (6500–7000 s). The fine velocity structure of turbulence is apparent. This type of linear instability for the Ekman solution for the boundary layer was investigated based on results of linear stability analysis and numerical simulations (Leibovich and Lele 1985; Wakata 2010, 2013).

We show the horizontally averaged basic field hereafter. The horizontal mean TKE dissipation rate is presented in Fig. 2. After breaking of instability waves, the bottom boundary layer develops and the bottom boundary layer thickness increases. The Ekman layer depth is regulated by the Coriolis force, but it is still gradually eroding the upper stratified region and expanding upward because no restoration mechanism exists for the temperature field. The vertical profile of the velocity is shown in Fig. 3a. The east–west flow is defected from the geostrophic current under 50-m height above the bottom (hab). Also, the north–south flow appears there, which is related to the Ekman spiral’s structure near the seabed. The vertical profile of the potential temperature becomes steep lower than 50 m hab (i.e., the fluid is well mixed vertically by turbulence) (Fig. 3b). GS and SGS kinetic turbulent energy are large there (Fig. 3c). Turbulence develops with consumption of the energy of these basic fields. The external forces, maintaining the geostrophic flow related to the ocean surface slope, become an energy source for the kinetic energy.

### b. Osborn’s method

We next investigate the applicability of Osborn’s method. Aside from Osborn’s method, the vertical diffusivity is calculable directly from the buoyancy flux, which is obtainable from LES study. To evaluate Osborn’s method, we regard this value as a “truth” value for comparison.

Hereafter, variables with an overbar denote a resolved-scale GS variables or grid mean representative values of SGS fluxes, whereas those with a single prime denote the SGS variables following the notation of Heinze et al. (2015). GS variables are divisible into the horizontal mean variable and deviation from it. The horizontal mean is denoted using angle brackets. Deviation from it is shown with a double prime. Then, a variable can be presented as .

The density flux caused by vertical motion can be estimated from the contribution from GS and SGS buoyancy fluxes (Fig. 4c). Therefore, the true diffusivity coefficient is obtained from

On the other hand the diffusivity of Osborn’s method is obtainable from (4). Figure 5 presents diffusivity estimated using Osborn’s method with and using the direct method. Osborn’s estimation is greatly overestimated near the seabed, although it is slightly underestimated in the upper part of the boundary layer. As discussed in a report by Osborn (1980), Osborn’s method is expected to give a possible maximum value for the diffusivity, but it is underestimated in the upper boundary layer. The difference of and arises from the nonuniform vertical profile of the dissipation rate (Fig. 4a). The dissipation rate is large near the seabed and small in the upper side of the boundary layer because the strong vertical shear near the seabed and small shear for the upper part as usually found in many observations (Simpson et al. 1996; Rippeth et al. 2003; Tsutsumi and Matsuno 2012; Endoh et al. 2016). This vertical profile entails a high diffusion coefficient in the lower part of the boundary layer for . The small Brunt–Väisälä frequency near the seabed, which appears in the denominator of the diffusion coefficient definition in (4), also gives rise to the large diffusion coefficient. By contrast, has a peak in the middle of boundary layer because the heat flux is small in the lower and upper parts of the boundary layer (Fig. 4c) because the vertical motion related to the heat flux is reduced near the seabed and at the top of the boundary layer by the influence of the seabed and the large stratification at the top of the boundary layer, respectively (Fig. 4b).

The dissipation flux coefficient can be obtained directly from the definition not using in (4) as

which is presented in Fig. 6 at each time. This is regarded as a “true” value for comparison. It tends to be zero near the seabed and to be larger than 0.2 in the upper part of boundary layer. It tends to be greatly varied above the boundary layer. The lower part of the vertical profile is elongated as the mixed layer develops, but the general view of its pattern is not varied largely for each time profile. Wide variation above the boundary layer results from the effects of internal gravity waves emitted from the bottom boundary layer rather than the local instability. Therefore, the use of constant for Osborn’s method is not correct.

### c. Energy budget

The TKE budget is discussed here to investigate the validity of the assumptions for introducing Osborn’s method: temporal and spatial uniformity. The budget analysis has been studied precisely for the PALM model by Heinze et al. (2015), who discussed the cloud-topped boundary layer for the atmosphere. Turbulence is produced by surface heating.

As described in the introduction, Osborn’s method depends on the energy budget equation. We describe the budget related to the sum of GS and SGS TKEs. The notation is presented by Heinze et al. (2015) as

where the left-hand side () represents the time derivative of the sum of GS and SGS TKEs, the first term () on the right-hand side is the TKE term representing generation caused by the vertical shear of horizontal mean basic flow, the second () represents baroclinic energy conversion, the third () is the advection, the fourth () denotes the pressure work, and the last () is the kinetic energy dissipation rate. Each term’s vertical profile is depicted in Fig. 7. The sum of and are depicted as a purple line (*T*), but is quite small (almost negligible). These are obtained after average time from to s.

The residual part is shown as a light blue line (*R*). The energy balance in (10) is not sufficiently attained quite near the seabed, which might be affected by the inclusion of the Prandtl boundary layer model for the first grid from the seabed. In other areas, the energy balance is quite well established. Heinze et al. (2015) noted the presence of a large residual part. It is an important difficulty deciding the part we should include this into. First, we had a similar problem of a large residual. However, results have demonstrated that this problem is solvable by taking a quite small time difference interval for the numerical simulation. We adopted s instead of the standard one s in PALM. This residual derives from the numerical error of the time derivative because the time derivative of TKE () is quite small but the time derivative on each grid point can be large because of the vital turbulent fluid motion. Energy production is balanced mainly with dissipation, and is partly balanced with baroclinic conversion and advection terms. The advection term *T* is understood not to be smaller than the baroclinic conversion term *B*, which indicates that the presumption for the Osborn’s method is not satisfied and that some correction is necessary. Some observations also suggest a similar energy budget balance. Shaw et al. (2001) and Reidenbach et al. (2006) show the main balance of production and dissipation. Feddersen et al. (2007) and Walter et al. (2014) show a slight excess of dissipation in the bottom mixing layer, which suggests the presence of advection *T*.

## 4. Alternative estimation formulas for vertical eddy diffusivity

In this section, we propose improvement of Osborn’s method and alternative estimation methods.

### a. Estimation based on the flux Richardson number

For nonuniformly distributed turbulence like that at the bottom boundary layer, Osborn’s method is not applicable, as noted in the previous section (Fig. 5). Therefore, we will return to the Osborn’s original definition of (4), in which the energy efficiency is not constant, but is instead a function of the flux Richardson number . Here, we estimate and evaluate the diffusivity obtained from this original definition.

The flux Richardson is obtainable from the LES results. It is defined as the ratio of the shear production term and buoyancy conversion as follows:

Figure 8 shows the flux Richardson number, which is not a constant value such as a critical Richardson number = 0.15 (Ellison 1957), but which decreases downward, indicating that the fluid motion is unstable near the seabed because of the shear instability. The term scatters widely above the boundary layer, where is not related to local instability but to nonlocal forcing. The vertical propagation of the internal gravity wave emitted from the top of the boundary layer was found after shear instability occurred at the first stage (not shown). The tendency of this large variation is somewhat exaggerated because of the assumption of vertical uniform mean flow of which the shear appears in the denominator of the definition of . However, this trend did not change greatly even when numerical experiments were conducted for the mean flow with slight shear. Both and are varied widely. The diffusivity rate [] estimated from this flux Richardson number is shown in Fig. 9a, which gives a better estimate than Osborn’s method with . However, it is apparently a slight underestimate. This discrepancy can be understood by reconsidering the energy budget. Because of the temporal and spatial uniformity, the time derivative term of TKE and the spatial derivative term expressed with a flux form were dropped when the equation of Osborn’s estimation method was derived. The former corresponds to *M*; the latter corresponds to *T* in Fig. 7. Although it is correct that *M* is small, *T* has a magnitude that cannot be ignored compared with the baroclinic conversion term *B*. Keeping the *M* and *T* terms, Osborn’s method would be modified as

Along Ivey and Imberger (1991), the generalized Richardson number is introduced as . Then, the dissipation coefficient can be expressed as without approximation. The generalized Richardson number and the flux Richardson number are shown in Fig. 10. Confined in the mixed layer, a linear relation is well established. The regression coefficient *α* is 1.19 (). Lozovatsky and Fernando (2013) obtained from atmospheric boundary layer data.

Then, the diffusive rate is obtainable from (12) as

Figure 9b presents the result of this modified method: The estimation line (red) almost overlaps the true line (black). The inadequacy related to the spatial nonuniformity can be corrected by multiplying a correction factor by .

### b. Estimation based on the gradient Richardson number

We can obtain vertical diffusivity accurately if we are able to obtain the flux Richardson number from some observation. However, to ascertain means that we already know the buoyancy flux (i.e., is already known and the estimation procedure in itself is no longer necessary). Although is difficult to observe, there is similar nondimensional number to the flux Richardson number; that is, the gradient Richardson number defined as

which is obtained more easily from the observation than . The gradient Richardson number is depicted by the red line in Fig. 8. The ratio affords the turbulent Prandtl number. This ratio is apparently not largely valid. The scatter diagram for and is shown in Fig. 11a. If is smaller than 0.2, then an almost linear relation prevails between and , but it scatters largely over 0.2 of . However, if we restrict the plotting in the boundary layer of Fig. 11b, there is a good linear relation . The regression analysis yields , as estimated under 30 m hab. This relation was investigated from laboratory experiments (Rehmann and Koseff 2004), numerical simulations (Rohr et al. 1984; Shih et al. 2005; Stretch et al. 2010; Venayagamoorthy and Koseff 2016), and observations (Davis and Monismith 2011). Several estimation equations have been proposed for this relation from the LES study. Kitamura et al. (2013) expressed it as using a Butterworth function. In the region for small , the linear relation apparently holds. Taking the Taylor expansion for this relation, it can be approximated as . Zilitinkevich et al. (2008) proposed the relation as , which can be approximated as . Those results are consistent with the results presented herein.

The relation between and is also represented in the scatter diagram in Fig. 12. The regression coefficient *β* is 1.79 . Mater and Venayagamoorthy (2014a) also showed the relation in light of the closure turbulent model (Yamada 1975; Mellor and Yamada 1982) from the data of the homogeneous turbulence (Shih et al. 2005 and the turbulence of atmospheric boundary layer (Lozovatsky and Fernando 2013). Their coefficient of the first-order Taylor expansion becomes 1.28, which is slightly smaller than the present result of 1.79.

The diffusivity coefficient is obtainable from the gradient Richardson number as

Figure 13 presents diffusivity estimated from the gradient Richardson number at each time for case A, which provides a very good approximation. However, the value was obtained from one case study. To investigate the degree to which this formula is applicable for various circumstances, some experiments were conducted with different conditions (Fig. 14): case B is for weak current, case C for weak stratification, case D for low latitudes, and case E for high latitudes, as listed in Table 1. Every case study shows good coincidence between and . If one observes the gradient Richardson number and the turbulent dissipation rate simultaneously, one can obtain an acceptable vertical diffusivity coefficient. Lozovatsky and Fernando (2013) obtained a similar expression for the dissipation flux coefficient as by analyzing meteorological data from Salt Lake City, Utah. Taking the Taylor expansion of (15), the coefficient of the first order of becomes 1.79, which is approximately their result of 1.7. Recently, Zhou et al. (2017) noted the same idea on the replacement of and for definition for stratified plane Couette flow, but they do not discuss the correction factor *β*.

### c. Estimation based on distance from the seabed

To obtain the diffusivity coefficient from the gradient Richardson number as in section 4b, velocity data are required. This requirement also imposes a heavy burden on ship observations. In addition, many already existing dissipation rate data without velocity data are inapplicable. Here, we try to derive an alternative method without using velocity data.

The dissipation flux coefficients for the standard case are presented in Fig. 6. They are small near the seabed because mixing up the density is difficult. Actually, the vertical density gradient is already weak (i.e., the density is already mixed up near the seabed, which functions as a thermal insulation wall). In addition, the vortex becomes smaller and the ability of mixing is weakened, because the turbulent vortex feels the seabed as a rigid wall and the fluid motion is suppressed there. Therefore, it is reasonable to infer that the dissipation flux coefficient is a function of distance *z* from the seabed. Examining for each time for case A in Fig. 6, the drawn lines are apparently scattered, but looking inside the bottom boundary layer, it is apparently monotonically increasing with z. This property is probably related to the vortex size. Therefore, again, the dissipation flux coefficient was drawn on the coordinate , standardized by the Ozmidov scale (Fig. 15). This Ozmidov scale is defined by on which scale the buoyancy force and the inertial force can be balanced. This scale is the maximum overturning length scale permitted in a stratified fluid. We use values for *ε* and *N* at each depth. Therefore, is a function of *z*. The redrawing line seems to converge to one line: the dissipation flux coefficient can be expressed using a similar function of .

This similarity function can be approximated by the following rational function:

In that equation, , which are obtained using the least squares method for . The real height corresponding to is 48.75 m at 3000 s for case A, which can sufficiently cover all depths of the bottom boundary layer.

The dissipation rate is determined using this as shown below:

The estimated dissipation rate from this method is shown in Fig. 16, which is slightly worse than that of previous methods, but which is well acceptable. To draw it above , we extrapolate it tentatively by adopting the constant value of for . A value is adopted practically over *z* corresponding to . The dissipation rate is quite small there, and tends to be a small value above the bottom boundary layer. We cannot see any signal in this horizontal axis scale there. The is widely deviated. Therefore, this value has no physical meaning and is useful only for convenience. The determining mechanism of the dissipation flux coefficient in a space freed from the shear instability is beyond the scope of this research, which particularly addresses the bottom boundary layer.

To confirm the applicability of this estimation from (16) for other circumstances, we investigate many cases presented in Table 1. The for each case and the approximated line from (16) are depicted in Fig. 17, in which the approximate line can represent each case well. The diffusivity coefficient estimated from (17) is presented in Fig. 18, which can sufficiently express the profile of the original diffusivity.

To estimate , it was obtained from the simulation result up to 100 000 s. For confirmation, we extended the simulation of standard case A up to 200 000 s. The vertical profile of (not shown) was redrawn. The result is quite similar to that shown in Fig. 17. The vertical structure of is steady for the time-independent forcing maintaining the geostrophic current. However, currents and their forcing in the ocean are generally varied with time. We supplemented a simulation in which the geostrophic current varied with time similarly to tidal flow. We assume that the tidal flow oscillates east–west with amplitude of 0.25 m s^{−1} (case F). The cycle corresponds to an tide with a period of 12.14 h. The tidal current is driven by the body force (Wakata 2013; Wakata et al. 2017). An analysis was conducted for data sampled at intervals of 1000 s from to s to calculate . During this period, the height of the bottom boundary layer varied repeatedly, 12 times, from almost 0 to about 30 m hab. Although the turbulent mixed layer thickness changes greatly, the vertical profile of is considerably steady (Fig. 19). The approximation function is expressed as shown below:

Representation with polynomials would be convenient to apply to observation data.

## 5. Summary and discussion

The method for estimating eddy diffusion from the dissipation rate in the bottom boundary layer was investigated from the study of LES. Results showed that the most commonly used estimation method reported by Osborn (1980) with constant dissipation flux coefficient of gives overestimation in the lower part of the bottom boundary layer and slight underestimation in the upper part. It cannot provide good estimation, which is consistent with Walter et al. (2014), who showed overestimation in the bottom boundary layer from the dissipation flux coefficient measured directly in the nearshore environment of southern Monterey Bay, California.

Osborn’s method was derived based on the TKE balance. When deriving this formula, two assumptions were made. One is temporal and spatial uniformity for the statistical distribution of the turbulence. According to this assumption, we can eliminate the terms with the time derivative and with the flux form, such as the advection term in the TKE equation. The other assumption is that included in defined in (4) was replaced with a constant value that is a critical value for shear instability. Consequently, becomes a constant value.

We evaluated the validity of the first assumption by examining the magnitude of the term neglected in the process of deriving Osborn’s method in the TKE equation. The neglected term is of about the same order as the buoyancy flux term. It cannot be ignored. The LES study showed that in the bottom boundary layer is not constant and that it varies greatly with depth. The second notion, that is constant, was also unacceptable.

The difficulty caused by the first assumption can be avoided, as correction is done simply by multiplying by a correction factor.

The difficulty presented by the second assumption is more substantial. Although can be readily obtained with LES, the derivation is difficult in a general numerical model such as a Reynolds averaged Navier–Stokes (RANS). Some empirical relational expressions are often used in the closure model adopted in RANS model. It is also extremely difficult to obtain expressions directly from ship observations. The present LES study shows that a linear relation exists between and in the bottom boundary layer. Using this relation, can be expressed as . Many case studies have demonstrated that this form can provide good estimation for various situations.

The reason for the great change in at the boundary layer might be related to the turbulence intensity. Very strong turbulence exists near the seabed, decreasing rapidly upward. From laboratory experiments and numerical simulations, it is noteworthy that several regimes have been found for the relation between diffusion and dissipation rate depending on the turbulence strength. The existence of those regimes means that changes depending on if . When the turbulence is small (i.e., the relation of Osborn is applicable). However, (Shih et al. 2005) or (Barry et al. 2001) holds when the turbulence becomes stronger. Recently, these relations have been studied by many researchers (Davis and Monismith 2011; Lozovatsky and Fernando 2013; Bluteau et al. 2013; Bouffard and Boegman 2013; Walter et al. 2014). The variety of might result from the fact that different regimes exist in the narrow area adjacent to the seabed.

An alternative method using no current data was proposed. Results show that can be expressed as a similarity function with the vertical coordinate normalized by the Ozmidov scale. This function was approximated by a rational function or a polynomial function. The universal form of was obtained. Results confirmed that this method can provide good estimation for various situations.

Recently, Scotti and White (2016) and Zhou et al. (2017) have applied the Monin–Obukhov similarity (M-O) theory for Couette flow turbulence, and Scotti and White (2016) can express the flux Richardson number as a function of the distance from the seabed scaled with the M-O length. Their simulation assumed a Couette flow without Coriolis force, with buoyancy flux imposed at the top and bottom boundaries. The entire area can belong to the M-O boundary layer because the current is driven by boundaries. In the present study, the current is forced by body force such as pressure force and Coriolis force. Strong shear is concentrated near the seabed. Furthermore, the buoyancy flux is zero at the boundary. Therefore, the use of M-O layer theory is inapplicable for the bottom ocean boundary layer. The M-O layer might be present in the thin layer under the Ekman layer as a logarithm bottom boundary layer.

For this study, we use the Ozmidov scale for normalizing the vertical coordinate. The for each time in Fig. 6 has originally similar vertical profiles. The difference derives from the elongated lower part, where buoyancy flux is small because of the weak stratification (i.e., already mixed up). As for developing the boundary layer, this lower part is elongated. Scaling of the vertical coordinate with the Ozmidov scale shrinks those parts. Therefore, a universal profile might be attained. This universality derives from the basic arrangement of the strong shear area concentrated near the seabed where TKE spreads upward and a moderate stratification area where TKE can be converted to potential energy. In this sense, this method will not be applicable in the case in which the heat flux such as horizontal advection or the surface flux is added. However, for many ocean bottom boundaries, it will be possible to provide a simple estimation method for the vertical diffusivity coefficient that is more suitable than the Osborn’s method with constant .

Recently, the functional dependency of *R*_{f} on several parameters were reviewed and examined (Gregg et al. 2018; Monismith et al. 2018). The present study showed that near the bottom boundary layer, the distance from the seafloor normalized by the Ozmidov scale can be an important scaling parameter for Γ (hence *R*_{f}). We must emphasize one more point: the method of estimating the vertical diffusion coefficient obtained from this study can only be applicable in the turbulent bottom boundary layer over a flat seafloor. In the ocean interior, further consideration of the effects of nonisotropic propagation of internal gravity waves might be necessary and is deferred to future studies.

## Acknowledgments

Valuable comments and help in the use of PALM from Dr. S. Raasch of Hanover University are greatly appreciated. Discussions of this work with Dr. T. Matsuno, Dr. T. Endoh, Mr. T. Nogami, and Mr. T. Matsuo of Kyushu University are also appreciated. The author thanks the anonymous reviewers, who provided constructive and thoughtful comments that greatly improved the paper. Computations were conducted mainly using computer facilities at the Research Institute for Information Technology, Kyushu University.

## REFERENCES

_{2}critical latitude in the Barents Sea

## Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).