NOVEMBER 2008

TAYLOR AND SARKAR

2535

Stratification Effects in a Bottom Ekman Layer JOHN R. TAYLOR

AND

SUTANU SARKAR

Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla, California (Manuscript received 29 October 2007, in final form 8 July 2008) ABSTRACT A stratified bottom Ekman layer over a nonsloping, rough surface is studied using a three-dimensional unsteady large eddy simulation to examine the effects of an outer layer stratification on the boundary layer structure. When the flow field is initialized with a linear temperature profile, a three-layer structure develops with a mixed layer near the wall separated from a uniformly stratified outer layer by a pycnocline. With the free-stream velocity fixed, the wall stress increases slightly with the imposed stratification, but the primary role of stratification is to limit the boundary layer height. Ekman transport is generally confined to the mixed layer, which leads to larger cross-stream velocities and a larger surface veering angle when the flow is stratified. The rate of turning in the mixed layer is nearly independent of stratification, so that when stratification is large and the boundary layer thickness is reduced, the rate of veering in the pycnocline becomes very large. In the pycnocline, the mean shear is larger than observed in an unstratified boundary layer, which is explained using a buoyancy length scale, u /N(z). This length scale leads to an explicit * buoyancy-related modification to the log law for the mean velocity profile. A new method for deducing the wall stress based on observed mean velocity and density profiles is proposed and shows significant improvement compared to the standard profile method. A streamwise jet is observed near the center of the pycnocline, and the shear at the top of the jet leads to local shear instabilities and enhanced mixing in that region, despite the fact that the Richardson number formed using the mean density and shear profiles is larger than unity.

1. Introduction The surface wind-driven Ekman layer forms when frictional terms contribute to the leading-order momentum balance, leading to an ageostrophic flow. Because the velocity gradients, and hence the size of the viscous stresses, are strongest at the sea surface, the ageostrophic component of the flow decreases with depth, leading to a turning of the mean velocity profile known as the Ekman spiral. The bottom Ekman layer, formed as a mean current flow over the seafloor, is directly analogous to the surface wind-driven Ekman layer. In a surface Ekman layer, the Ekman transport (the vertical integral of the Ekman layer velocity) is directed 90° to the right of the wind stress in the Northern Hemisphere. In a bottom Ekman layer, the Ekman transport is also directed to the right of the bottom stress. However, because the stress at the seafloor is in the opposite

Corresponding author address: Sutanu Sarkar, Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla, CA 92093. E-mail: [email protected] DOI: 10.1175/2008JPO3942.1 © 2008 American Meteorological Society

direction of the geostrophic current, the Ekman transport in the bottom boundary layer is 90° to the left of the geostrophic current, and the Ekman spiral turns counterclockwise with increasing depth in the Northern Hemisphere. Also like the surface mixed layer, the bottom boundary layer is a site of intense turbulent dissipation and mixing (Garrett et al. 1993; Thorpe 2005). Together with the surface mixed layer, the bottom boundary layer is a major “hotspot” for diapycnal mixing in the ocean (Thorpe 2005). It has long been hypothesized that mixing of the density field in bottom boundary layers may be important to abyssal mixing through along-isopycnal advection out of the boundary layer and boundary layer detachment (Munk 1966; Armi 1978). In addition, through frictional loss and Ekman pumping, the bottom boundary layer provides an important momentum sink for deep currents and mesoscale eddies. Three-dimensional numerical simulations of stratified bottom Ekman layers have been described previously by several authors but, in all cases, stratification was applied with a cooling heat flux at the lower wall, akin to the stable atmospheric boundary layer. Cole-

2536

JOURNAL OF PHYSICAL OCEANOGRAPHY

man et al. (1992) performed a direct numerical simulation (DNS) of a turbulent Ekman layer with a friction Reynolds number of Re ! 340 and a constant heat * flux at a smooth, no-slip lower wall. They compared the stratified Ekman layer to a previous study of an unstratified Ekman layer (Coleman et al. 1990) and found that the surface heat flux limits the transport of turbulent kinetic energy into the outer layer and broadens the Ekman spiral. Shingai and Kawamura (2002) considered the same flow at a friction Reynolds number Re ! 428.6. They found that the boundary layer thick* ness defined in terms of either the momentum or the buoyancy flux decreases sharply with the application of a surface heat flux. In general, these studies imply that when a strong, stable stratification is applied to the wall under an Ekman layer, stratification acts to suppress the turbulent production in the boundary layer, increasing the turning angle and decreasing the boundary layer height. In the ocean, with the exception of isolated hotspots, the seafloor can be assumed to be adiabatic. Changes in density then affect the Ekman layer through the stratification associated with the ambient water. Because a mixed layer can be found near the seafloor throughout most of the ocean, stratification can be expected to affect the boundary layer in a much different manner than in a typical stable atmospheric boundary layer, where a heat flux is often present at the ground. By comparing the expected unstratified turbulent Ekman layer depth to the thickness of bottom mixed layers in field data, it is apparent that stratification often limits the boundary layer height. In an unstratified Ekman layer, the boundary layer height is expected to be approximately h ! 0.5 u /f, where u is the friction ve* * locity and f is the Coriolis parameter. Typical midlatitude values, say u ! 1 cm s"1 and f ! 10"4 s"1, imply * an unstratified Ekman layer depth of about h ! 50 m. In many cases, especially in a coastal environment where the stratification is typically large, the observed mixed layer heights are often much smaller than 50 m. For example, Perlin et al. (2007) observed a mixed layer thickness of about 10 m in an Ekman layer over the Oregon shelf. It is not fully clear how the Ekman layer structure changes when the Ekman layer height is limited by stratification, and this will be one focus of the present study. The bottom boundary layer plays an important role in the drag induced on mean currents and mesoscale eddies. To obtain an accurate prediction of the ocean state, numerical ocean models must represent this loss of momentum. However, due to computational restrictions on the grid size, this is not straightforward. To accurately apply the no-slip boundary condition at the

VOLUME 38

seafloor, a numerical model must resolve the viscous sublayer. Because the viscous sublayer in the ocean is thin, O(0.1–10 cm) (Caldwell and Chriss 1979), numerical models are clearly unable to resolve this region, and an approximate boundary condition must be used. It is common practice to model the seafloor stress, which then provides a Neumann boundary condition for the horizontal momentum equations. For example, the Regional Ocean Modeling System (ROMS) provides three methods for modeling the bottom stress based on the velocity at the lowermost grid cell: linear and quadratic drag coefficients and a law-of-the-wall using a specified bottom roughness. A general form for the bottom stress using the linear drag coefficient #1 and the quadratic drag coefficient #2 is

!w,x ! $#1 % #2&u2 % $ 2'u, "0 !w,y ! $#1 % #2&u2 % $ 2' $, "0

$1'

where (w,x and (w,y are the zonal and meridional components of the bottom stress, respectively. Typical values of the bottom stress coefficients are #1 ! 2 ) 10"4 m s"1 and #2 ! 0 for linear bottom drag and #1 ! 0 and #2 ! 2 ) 10"3 for quadratic bottom drag (Haidvogel and Beckmann 1999). The friction velocity u ! &(w/*0 is often used in * scaling arguments for both bulk and turbulent properties, such as the boundary layer height, turbulent dissipation, Reynolds stresses, etc. Despite its first-order relevance, there remains some uncertainty about how to estimate the wall stress from observational data, especially when the velocity profile is affected by wall roughness and stratification. A common method used to evaluate the friction velocity is the so-called profile method (Johnson et al. 1994), which utilizes the classical law-of-the-wall (see, e.g., Pope 2000). Using this method, the friction velocity and wall roughness z0 are determined by fitting the observed velocity profile to the logarithmic profile |u| !

u

!"

z * ln , % z0

$2'

where |u| is the horizontal velocity magnitude and + ! 0.41 is the von Kármán constant. Johnson et al. (1994) applied this method to estimate the friction velocity at the bottom of the Mediterranean outfall plume. They found that estimates of the friction velocity depended on how far away from the bottom the fit was applied, even when restricted to the bottom mixed layer. This implies that the standard law-of-the-

NOVEMBER 2008

2537

TAYLOR AND SARKAR

wall was not valid throughout the bottom mixed layer. Dewey et al. (1988) compared several methods for estimating the bottom stress using microstructure profiles over a continental shelf. They found that the wall stress estimated using the profile method was a factor of 4.5 larger than that estimated using the dissipation method, with the friction velocity given by u ! "!"z#1#3, *

"3#

where again $ ! 0.41 is the von Kármán constant. The authors speculated that this discrepancy may be the result of form drag induced by local bedforms. Stahr and Sanford (1999) obtained velocity and dissipation measurements in the North Atlantic deep western boundary current and similarly found that the wall stress estimated from the profile method was 3 times larger than the wall stress estimated using the dissipation method. Perlin et al. (2005) also found that the profile method gave large values of the wall stress compared to the dissipation method and proposed that the elevated mean shear could be explained by the influence of stratification on the boundary layer structure. They proposed that when the local stratification is sufficiently large, stratification, not distance from the wall, limits the size of the largest turbulent eddies. To quantify this hypothesis, they used an empirical function involving the Ozmidov scale LOz ! (%/N3)1/2, where % is the turbulent dissipation rate. Estimates of the wall stress using this method, which will be referred to here as the “modified law-of-the-wall,” gave much better agreement with the dissipation method. Because direct field observations of dissipation require instruments capable of capturing finescale velocity fluctuations, it would be desirable to have an accurate method for estimating the wall stress from more commonly observed quantities such as the velocity and density profiles. Weatherly and Martin (1978) considered a stratified Ekman layer using a one-dimensional numerical model and the Mellor–Yamada level-II closure to parameterize turbulent mixing. They assumed that the flow outside the boundary layer was uniformly stratified, steady, and in geostrophic balance; the same assumptions that will be made in the present study. Near the wall, a mixed layer formed, which was separated from the outer layer by a strongly stable pycnocline. They found that the thickness of the bottom Ekman layer was strongly limited by the presence of a stable stratification outside the boundary layer. When the outer layer buoyancy frequency was N&/f ! 200, they found that the angle made by the surface stress relative to the free-stream velocity was '0 ! 27°, nearly twice the

value in an unstratified boundary layer of '0 ! 15°. Most of the change in the turning angle occurred in the pycnocline, and the flow in the mixed layer was nearly unidirectional. Perlin et al. (2007) observed a lower turning angle of 15° ( 5° in a stratified bottom Ekman layer over the Oregon shelf. In agreement with the simulations by Weatherly and Martin (1978), they found that nearly all of the Ekman transport occurs in the relatively thin bottom mixed layer. Weatherly and Martin (1978) proposed a scaling for the height of a stratified turbulent Ekman layer given by hWM ! A

u * f

!

1)

N2$ f2

"

*1#4

,

"4#

where A # 1.3 as determined empirically by the onedimensional model results. In the limit of an unstratified Ekman layer, Eq. (4) gives a significantly larger height than the conventional value of h ! 0.4 * 0.5+. This is a result of the choice by Weatherly and Martin to use the criteria that the turbulent kinetic energy becomes zero at the top of the boundary layer. Choosing the location where the velocity magnitude becomes equal to the free stream would have yielded a boundary layer height nearly a factor of 3 smaller than that obtained using the turbulent kinetic energy (TKE) criteria when the outer flow is unstratified. A somewhat different scaling law was proposed by Zilitinkevich and Esau (2002): hZE ! CR

u

* f

!

1)

C2RCuN N$ f C2S

"

*1#2

,

"5#

where the constants CR ! 0.5 and CuN /C2S ! 0.6 were determined by fitting to data from a large eddy simulation (LES). The functional form in Eq. (5) was found to be an adequate fit to field observations by Zilitinkevich and Baklanov (2002). When a stable stratification is found outside of a turbulent well-mixed region, internal waves generated by the interaction between the turbulent eddies and the stratification are possible. Turbulence-generated internal waves have been found in laboratory experiments and numerical simulations in a wide variety of flows: shear layers, gravity currents, boundary layers, etc. Taylor and Sarkar (2007a) observed turbulence-generated internal waves in numerical simulations of a stratified Ekman layer over a smooth wall at a modest Reynolds number. They found that the vertical energy flux associated with the upward-propagating internal waves was small compared to the integrated boundary layer dissipation but was of the same order as the integrated buoyancy flux. Because the buoyancy flux is respon-

2538

JOURNAL OF PHYSICAL OCEANOGRAPHY

sible for the transfer of turbulent kinetic energy to the potential energy field, this implies that the energy radiated by internal waves is comparable to that used to mix the background density field. As in many previous studies, they found that the waves propagating through the outer layer were associated with a relatively narrow band of frequencies leading to vertical propagation angles between 30° and 60°. Because the waves are generated in a turbulent region with a wide range of spatial and temporal scales, this result is remarkable. The authors described a model based on viscous ray tracing that was used to predict the decay in amplitude of a wave packet after it had traveled a given distance from the source. They found that this relatively simple linear model was able to capture many characteristics of the observed frequency spectrum of the internal waves in the outer layer, including the range of dominant propagation angles. This paper will be organized as follows: the governing equations and physical approximations are discussed in section 2, and the numerical method used to evolve the governing equations is presented in section 3. Results from the simulations will be separated into several sections: evolution of the mean density and velocity profiles will be discussed in section 4, boundary layer turbulence will be discussed in section 5, methods for estimating the friction velocity from field data will be evaluated in section 6, and turbulence-generated internal waves will be considered in section 7.

2. Formulation The turbulent Ekman layer considered here is formed when a steady flow in geostrophic balance encounters a nonsloping, adiabatic lower wall, as illustrated in Fig. 1. The primary objective of this study is to consider a controlled environment in which we can examine the influence of the outer layer stratification on a turbulent Ekman layer. Therefore, we make several simplifying approximations. The free stream is assumed to be in geostrophic balance and aligned with the x axis. The seafloor is represented by a nonsloping, rough surface. The roughness elements are too small to be resolved directly by the grid, but their effect is parameterized through a near-wall model. The lateral boundaries are periodic, which is consistent with the assumption that the flow is statistically homogeneous in the horizontal plane and that the domain is large enough so that the flow is decorrelated over a distance equal to the domain size. In outer units, D ! U"/f, the domain size is 0.108D # 0.108D # 0.0405D in the x, y, and z directions, respectively. If we assume that U" ! 0.0674 m s$1 and f ! 10$4 rad s$1, then D ! 674 m and the domain size is 72.8 m # 72.8 m # 27.3 m. As a result

VOLUME 38

FIG. 1. Schematic of computational model. Dimensional parameters can be obtained by assuming U" ! 0.0674 m s$1 and f ! 10$4 rad s$1. The domain size is 72.8 m # 72.8 m # 27.3 m. Three values of outer layer stratification are considered: N" /f ! 0, 31.6, and 75.

of the horizontal periodicity, the mean vertical velocity must be zero. Therefore, the features of oceanic boundary layers owing to Ekman pumping/suction driven by large-scale horizontal gradients in the outer flow will not be present here. The goal of an LES is to accurately solve a low-passfiltered version of the Navier–Stokes equations. Because the filtered version of the nonlinear advection terms involves the total and filtered velocity fields, we are left with fewer equations than unknowns, the wellknown turbulence closure problem. A model is then necessary to write the residual stresses in terms of the filtered quantities. After normalizing with the freestream velocity U", the length scale, D ! U"/f, and the outer layer density gradient, d%/dz", the LES-filtered governing equations can be written as 1 !u & u " !u ! $ !p$ & f kˆ # 'ıˆ $ u ( !t #0 1 #$ $ Ri% kˆ & ! 2 u $ ! " &, #0 Re% 1 !#$ & u " !#$ ! !2#$ $ ! " ", !t Re%Pr

'6( '7(

!"u!0

'8(

where the Reynolds number, Richardson number, and Prandtl number are defined by Re% !

U%D , '

Ri% ! $

g d# D2 , #0 dz% U 2%

Pr !

' , (

'9(

%0 is the constant density used to define the momentum using the Boussinesq approximation, ) is the molecular kinematic viscosity, * is the molecular diffusivity, and # and " are the subgrid-scale stress and density flux, re-

NOVEMBER 2008

2539

TAYLOR AND SARKAR

TABLE 1. Relevant physical parameters. Re,, N,/f, Pr, and z0 are input parameters and the other parameters are outputs of the numerical model. Re, 7

4.55 . 10 4.55 . 107 4.55 . 107

Re* 1.08 . 105 1.09 . 105 1.12 . 105

Ri* 0 1000 5625

N,/f 0 31.6 75

spectively. The density and pressure have been decomposed into a plane average plus a fluctuation; that is, ! " #!$ % !&. The hydrostatic pressure gradient and the plane-averaged buoyancy force are in balance and do not appear in Eq. (6). An alternative normalization can be carried out using the friction velocity u " '(w/!0. * This leads to the friction Reynolds and Richardson numbers: u ! u g d# !2 N2$ * * Re " , Ri " ) , ! " . " 2 * * " #0 dz$ u2 f f * *10+ Relevant input and output nondimensional parameters are listed in Table 1. The dimensional parameters U,, N,, f, -, and the roughness length scale z0 are inputs to the simulations. Note that Re, is the same for each simulation, but the drag coefficient depends on the outer layer stratification, and hence the friction Reynolds number varies between each case. We have performed simulations at three different values of Ri , * equivalent to changing the free-stream density gradient. For comparison with oceanographic conditions, observations of the bottom boundary layer over the Oregon shelf by Perlin et al. (2007) provide estimates of Re " 4 . 104 ) 8 . 105 and N,/f " 75. Therefore, * both the Reynolds numbers and stratification levels considered in the present study are comparable with the field data of Perlin et al. (2007). To provide dimensional scalings for our simulations, we will use U, " 0.0674 m s)1, f " 10)4 s)1, and - " 10)6 m2 s)1, which, assuming that u /U, " 0.049, yields / " u /f 0 33 m. * * The applied roughness length scale z0 " 0.16 cm is consistent with observations by Perlin et al. (2005), who found that z0 " 0.05–2 cm, depending on the method used to infer the wall stress.

3. Numerical methods Simulations have been performed using a computational fluid dynamics solver developed at the University of California, San Diego. The algorithm and numerical method are described in detail in Bewley (2008). Because periodic boundary conditions are applied in the horizontal directions, derivatives in these directions are

Pr 5–10

u*/U, 0.0488 0.0490 0.0497

z0//

20

4.80 . 10 4.78 . 10)5 4.71 . 10)5 )5

15.4° 18.9° 24.8°

computed using a pseudospectral method, while derivatives in the vertical direction are computed with second-order finite differences. Time stepping is accomplished with a mixed explicit/implicit scheme using third-order Runge–Kutta and Crank–Nicolson. It can be shown that the numerical scheme ensures the discrete conservation of mass, momentum, and energy. To prevent spurious aliasing due to nonlinear interactions between wavenumbers, the largest 1/3 of the horizontal wavenumbers are set to zero, the so-called 2/3 dealiasing rule (Orszag 1971). To prevent the formation of spurious energy near the grid scale, a low-pass spatial filter is applied to the velocity and temperature fields. A fourth-order compact filter with a sharp wavenumber cutoff (Lele 1992) is applied in the vertical direction every 10 time steps. An open boundary condition is implemented at the top of the computational domain to prevent spurious reflections of upward-propagating internal gravity waves. The combination of a radiation boundary condition (Durran 1999) and a sponge-damping region that is used here was also used by Taylor and Sarkar (2007a), who found that only 6% of the wave energy was reflected back from the open boundary. The subgrid-scale stress tensor ! and the subgridscale density flux " in Eqs. (6) and (7) are evaluated using the dynamic Smagorinsky model. The Smagorinsky coefficients, C and C!, are evaluated using the dynamic procedure as formulated by Germano et al. (1991). As in Germano et al. (1991), the Smagorinsky coefficients have been averaged over horizontal planes. After averaging, if either of the Smagorinsky coefficients is negative, it is truncated to zero. This effectively prevents backscatter of energy from the subgrid scales to the resolved scales and helps to ensure numerical stability (Armenio and Sarkar 2002). While it is more computationally intensive than other methods, the dynamic Smagorinsky model has been chosen here because it avoids the empirical specification of the Smagorinsky coefficient and has been shown to perform well for wall-bounded and density-stratified flows (e.g., Armenio and Sarkar 2002; Taylor et al. 2005). To ensure accuracy of the solution obtained using LES, it is generally necessary to resolve the turbulent scales responsible for a substantial portion, say, 150%,

2540

JOURNAL OF PHYSICAL OCEANOGRAPHY

of the turbulent kinetic energy. Near walls, this criterion becomes increasingly stringent as the Reynolds number increases. To avoid the need to resolve the very small turbulent motions near the lower wall, we have used a near-wall model to estimate the wall stress based on the resolved velocity at the first grid point. We have used a model proposed by Marusic et al. (2001), which Stoll and Porte Agel (2006) found works well for high Reynolds number, rough wall boundary layers. To compensate for known deficiencies in the Smagorinsky LES model at very large Reynolds numbers, the near-wall model was augmented by a novel adaptive stochastic forcing procedure. The forcing amplitude is small, less than 6% of !w/!t, and limited to the first nine grid points near the wall. This technique is described in detail and validated in Taylor and Sarkar (2007b). The flow has been discretized using 128 grid points in the horizontal directions and 201 points in the vertical direction. The horizontal grid spacing in wall units was # "# x $ "y $ 1856 before dealiasing. The minimum vertical grid spacing, which occurs at the lower wall, was "z# (1) $ 121, while the maximum grid spacing is "# z $ 1075, or "z /D $ 4.86 % 10&4. The first grid point, where the near-wall model is applied, was located 60.7 wall units away from the wall, which is within the expected logarithmic region. The grid spacing at the top of the mixed layer is always sufficient to resolve the local Ellison scale, defined as

LE $

'!"2(1#2 . d)!*#dz

'11(

It has been found by Taylor and Sarkar (2008) that the Ellison scale provides a reasonable estimate for the scale of the turbulent eddies that are primarily responsible for entrainment in to the bottom mixed layer. Resolving the Ellison scale at the top of the boundary layer is therefore necessary in order to capture the growth of the mixed layer due to turbulent entrainment (Taylor and Sarkar 2008). The vertical grid spacing and the Ellison scale for each of the stratified simulations presented here are shown in Fig. 2. The asterisks mark the location where the local Reynolds stress is 10% of its maximum value, which provides a good estimate for the boundary layer height. To obtain an initial condition for the velocity field, a low Reynolds number unstratified simulation was conducted until the flow reached a steady state. The velocity field from this simulation was then interpolated onto a finer grid and a simulation at a higher Reynolds number was continued until all transients had decayed. The stratified simulations were initialized using the steady-

VOLUME 38

FIG. 2. Vertical grid spacing and the Ellison scale. Asterisks mark the top of the boundary layer, defined as the location where the local Reynolds stress is 10% of its maximum value.

state unstratified velocity field and an undisturbed, piecewise-linear temperature profile. The initial temperature profile has a relatively thin mixed layer with a thickness of about z ! 0.04+. This was done so that the stochastic forcing, mentioned above, would not produce spurious internal waves by forcing a stratified region. For z , 0.08+, the temperature profile was set to - $ zd-/dz|., and for 0.04 / z/+ / 0.08, a linear profile was used to make the integrated heat content over the domain the same as if - $ zd-/dz|. everywhere. Each stratified simulation was run for about tf $ 1.5 to allow the unstratified turbulence levels to adjust to the imposed stratification. This creates a better initial velocity field, and the temperature field was then reset to the piecewise-linear profile. The time at which the temperature field was reset will be referred to as t $ 0.

4. Mean boundary layer structure The time history of the plane-averaged temperature gradient is shown in Fig. 3. After the flow is initialized, the mixed layer near the wall grows rapidly. As the flow develops, a strongly stratified pycnocline forms above the mixed layer. It can be shown that if there is no net heat flux into the domain (which is a good approximation given that the molecular heat flux through the top of the domain is small), then a pycnocline is necessary to maintain heat conservation after the formation of a mixed layer. At the start of each stratified simulation, two distinct pycnoclines are visible for a short time, one above the bottom mixed layer and another above a

NOVEMBER 2008

2541

TAYLOR AND SARKAR

FIG. 3. Evolution of the plane-averaged temperature gradient.

region where the residual motions from the initial velocity field mix the temperature field. Eventually these pycnoclines merge, and the resulting single pycnocline gradually moves away from the wall as the mixed layer grows. In the case with N!/f " 75, the pycnocline has a larger stratification, 4 times the outer layer value, compared to the case with N!/f " 31.6. After about tf " 6, the temperature fields in both stratified simulations reach a quasi-steady state where the temperature gradient in the pycnocline does not grow further and the mixed layer growth is relatively slow. Two averaging methods will be used for velocity- and temperature-dependent fields. The Reynolds average, denoted by angle brackets, will be taken as an average over a horizontal plane and in time. To remove any bias due to inertial oscillations, averages are taken over one inertial period after the flow has reached a quasi-steady state. However, because the mixed layer thickness continues to increase in time, albeit slowly, temporal averages are not appropriate for the thermal fields and would, for example, lead to a smearing out of the pycnocline. Temperature will therefore not be averaged in time, but the plane average of the temperature-dependent fields will be taken at the time corresponding to the center of the Reynolds average window, unless otherwise noted. The plane-averaged temperature profiles after the flow that have reached quasi-steady state are shown in Fig. 4a. Both profiles correspond to a time t " 9.4/f after

initializing the temperature field. The components of the Reynolds-averaged horizontal velocity are shown in Fig. 4b. Compared to the temperature profiles, it is apparent that most of the Ekman transport is confined to the mixed layer. The Ekman layer height is reduced by the outer layer stratification at a level that is quantitatively consistent with the scaling law of Zilitinkevich and Esau (2002), given in Eq. (5). The magnitude of the cross-stream velocity increases when the outer layer is stratified, resulting in a broadening of the Ekman spiral, as shown in Fig. 5. It is interesting to note that although the density gradient is zero near the lower wall, the surface turning angle #0 " tan$1(%y /%x) increases with the outer layer stratification, as listed in Table 1. The angle of Ekman veering, defined by

! " tan$1

! "

&"' , &u'

(12)

is shown as a function of z/* in Fig. 6. When the flow is unstratified, the turning of the mean velocity occurs gradually throughout the boundary layer. The rate of turning in the mixed layer does not depend strongly on the outer layer stratification, and when the boundary layer height decreases significantly, the rate of turning, d#/dz, becomes very large in the pycnocline. Despite the fact that the boundary layer height changes significantly, the Ekman transport is nearly independent of the outer layer stratification. An expres-

2542

JOURNAL OF PHYSICAL OCEANOGRAPHY

VOLUME 38

FIG. 4. (a) Instantaneous plane-averaged temperature profiles. Thick lines show the profiles at the center of the velocity- averaging window; thin lines show the profiles at the start and end of the averaging window. (b) Plane and time-averaged horizontal velocity components, averaged over one inertial period.

sion for the Ekman transport can be found by integrating the mean streamwise momentum equation. Assuming that the mean flow is steady, and neglecting the molecular viscosity,

!

!

0

u2 * !"" dz # cos$#0%. f

$13%

As seen in Table 1, both u /U&, and '0 increase with * N&/f. Their effect partially cancels, so that the Ekman transport normalized by U& and f is nearly independent of N&. Specifically, evaluation of the right-hand side of Eq. (13) yields a transport of 0.107 m2 s(1, 0.105 m2 s(1, and 0.101 m2 s(1 for N& /f # 0, 31.6, and 75, respectively, where the transport has been made dimensional by taking U& # 0.0674 m s(1 and f # 10(4 rad s(1. The same result could be obtained by numerically integrating the cross-stream velocity profiles. As has been found in previous studies, the large density gradient at the top of the boundary layer coincides with an increase in the mean shear. The temperature gradient normalized by the outer layer value is shown in Fig. 7a, and the Reynolds-averaged velocity and shear profiles are shown in Figs. 7b–c. A peculiar feature of the mean velocity profile when N&/f # 75 is that

the mean velocity and the mean shear are maximum at the same location near the center of the pycnocline. It can be shown that this is the result of the rapid rate of veering that occurs in the pycnocline in this case. Using the definition of the Ekman veering angle, the square of the mean shear can be written as

" # " # " # d!u" dz

2

)

d!"" dz

2

#

d!|u|" dz

2

) !|u|"2

" # d# dz

2

,

$14%

where | u| # *u2 ) +2. At the center of the pycnocline, the first term on the right-hand side is zero because the velocity magnitude |u| is maximum at this location.

FIG. 5. Mean velocity hodograph showing the Ekman spiral.

NOVEMBER 2008

2543

TAYLOR AND SARKAR

cases are in reasonably good agreement with the unstratified logarithmic law. Deviations from the logarithmic velocity profile in the cases when stratification is present can be seen clearly by plotting the normalized velocity gradient !#

| |

"z d!u" . u dz *

$16%

The quantity ' can be interpreted as the ratio of the observed mean velocity gradient to that expected from the logarithmic law and is shown in Fig. 9. When the outer layer is stratified, the mean shear in the pycnocline increases significantly compared to the log-law value. It is worth noting that deviations from the lawof-the-wall begin well within the mixed layer. The consequences for this observation will be discussed in section 7.

5. Boundary layer turbulence

FIG. 6. Ekman veering angle.

However, because the Ekman veering rate is very large at this location, the second term leads to a local maximum of the mean shear in the pycnocline. Most of the enhanced shear in the lower portion of the pycnocline is due to the spanwise shear, while the streamwise shear dominates in the upper portion of the pycnocline. Coleman (1999) showed that in an unstratified turbulent Ekman layer, the mean velocity magnitude should follow the classical logarithmic law (or law-ofthe-wall). The law-of-the-wall is expected to hold in a region far enough from the wall where viscosity can be neglected, but near enough to the wall so that the boundary layer depth is not felt. In this case, the only relevant length scale must be the distance from the wall, so that

| ddz!u" | # l* , u

$15%

where l # &z and & is the von Kármán constant that is empirically found to be about 0.41. The Reynoldsaveraged horizontal velocity magnitude is plotted in Fig. 8 on a semilogarithmic scale. Very near the wall, all

To understand the increase in mean velocity and mean shear in the pycnocline, it is necessary to consider the influence of stratification on the turbulent eddies. The turbulent kinetic energy in the mixed layer scales with the friction velocity u , which as we have seen does * not change significantly with the addition of an outer layer stratification. In addition, the density gradient in the pycnocline appears to scale with the outer layer stratification. Therefore, as the stratification in the outer layer increases, the kinetic energy associated with eddies at the top of the mixed layer remains about the same while the potential energy required for an eddy to overturn increases. Therefore, when the stratification becomes large enough, turbulence is inhibited at the top of the boundary layer, and as a result, the growth of the mixed layer is strongly limited. One consequence of the damping of turbulence by stratification is a decrease in the turbulent stresses, !u(w(" and !) (w(". The decrease in boundary layer height with increasing stratification is very apparent from the Reynolds stress profiles shown in Fig. 10. Here, the Reynolds stress includes both the resolved and subgrid-scale contributions. Note that above the boundary layer, the Reynolds stress approaches a small, nonzero value owing to the presence of turbulencegenerated internal gravity waves. It is evident from Fig. 10 that the Reynolds stresses change more rapidly with height when the outer layer stratification is stronger. Changes in the Reynolds stress with height result in a momentum flux that must be balanced by other terms at a steady state. The steady-state plane-averaged horizontal momentum equations can be written as

2544

JOURNAL OF PHYSICAL OCEANOGRAPHY

VOLUME 38

FIG. 7. (a) Temperature gradient, (b) horizontal velocity magnitude, and (c) mean shear.

!!u" !2!u" !!u"w"" , #0#$ % f !#" % $ !t !z dz2

&17'

!!#" !2!#" !!#"w"" . #0#$ $ f &!u" $ U%' % $ !t !z dz2

&18'

We have found that the Reynolds-averaged eddy viscosity is cogradient. Hence, in the upper portion of the Ekman layer where d!("/dz ) 0, the spanwise Reynolds stress is positive as evidenced in Fig. 10. Because both components of the Reynolds stress become nearly zero in the outer layer, there must be a region near the top of the Ekman layer where d!( *w*"/dz ) 0. When the leading-order y-momentum balance is between the Reynolds stress and Coriolis terms, the region with d!(*w*"/dz + 0 leads to !u" + U,. Therefore, a zonal jet is an inherent feature of a steady-state turbulent Ekman layer. The turbulent kinetic energy budget can be found by dotting u* into the perturbation momentum equations and taking the plane average, to give

where k # !u*i u*i "/2 is the turbulent kinetic energy. Reading from left to right, the terms on the right-hand side of Eq. (19) can be identified as the turbulent transport, pressure transport, viscous diffusion, production, dissipation, buoyancy flux, subgrid transport, and subgrid dissipation, respectively. The leading terms in the turbulent kinetic energy budget for N,/f # 75 are shown in Fig. 11. Near the wall, the leading-order balance is between production and dissipation and does not differ significantly from the unstratified case as shown in the inset. In the upper portion of the mixed

!k ! 1 1 !2k 1 ! #$ !w"u"i u"i " $ !w"p"" % !t 2 !z !z &0 Re !z2 * $ !Sij"!u"i u"j " $

$

! !u"' " % !z i 31

1 Re

*

!

! " 'ji

"

!u"i !u"i $ Ri !w"&"" * !xj !xj

!u"i , !xj

&19'

FIG. 8. Reynolds-averaged horizontal velocity magnitude.

NOVEMBER 2008

2545

TAYLOR AND SARKAR

FIG. 9. Nondimensional velocity gradient, , " (-z/u*) d#u%/dz.

layer, the turbulent transport appears as a source term representing the advection of turbulent eddies toward the pycnocline. Pressure transport, buoyancy flux, and dissipation act as energy sinks in the upper mixed layer. In the pycnocline, starting at about z/! " 0.15, the turbulent transport and buoyancy flux decrease as stratification suppresses turbulent motion while pressure transport becomes the dominant source term. When the pressure transport is positive the vertical energy flux, # p$w$% increases with height, consistent with an internal wave field that is gaining energy. This is direct evidence of the generation of internal waves by the interaction between boundary layer turbulence and a stable stratification. The properties of these waves will be discussed in section 6. The instantaneous temperature field is shown for N&/f " 75 in Fig. 12 in an x–z plane at t " 11.2/f. Figure 12b shows an enlarged version of the box drawn in Fig. 12a. Isotherms are drawn every 0.025!d#'%/dz&. The mixed layer is very homogeneous with disturbances rarely exceeding the contour level. Turbulence-generated internal waves are visible as disturbances of the isotherms in the outer layer, which is generally statically stable with the notable exception of the region just above the pycnocline. Density overturns are visible both below and above the pycnocline. As highlighted in Fig. 12b, the overturns above the pycnocline are reminiscent of Kelvin–Helmholtz billows associated with a negative mean shear. The irreversible mixing resulting from these overturns may be responsible for the reduced temperature gradient above the pycnocline as seen in Fig. 7a. The stability of the shear above the pycnocline can be

FIG. 10. Reynolds stress profiles.

examined using the gradient Richardson number. The mean gradient Richardson number, defined as

#Rig% "

#N2% " , )d#u%!dz*2 + )d#$%!dz*2 #S2% ()g!"0*d##%!dz

)20*

FIG. 11. TKE budget at the top of the mixed layer and pycnocline for N&/f " 75. Inset shows the TKE budget near the wall.

2546

JOURNAL OF PHYSICAL OCEANOGRAPHY

VOLUME 38

FIG. 12. Isotherms projected onto an x–z plane for the case with N&/f ' 75 (a) full computational domain and (b) zoom of boxed region near the pycnocline.

as shown in Fig. 13a. Heaving of the pycnocline, which is visible in Fig. 12, causes significant variations in the pycnocline height. Because we are interested in the shear and stratification a small distance above the pycnocline, averaging over a constant height could be problematic. Therefore, the average operator in Eq. (20) has been taken with respect to a constant distance from the center of the pycnocline (identified by the maximum temperature gradient). It has been shown that a stratified flow can develop linear shear instabilities if Rig is less than 0.25 somewhere in the flow (Miles 1961; Howard 1961). In our simulations !Rig" # 0.25 in the mixed layer, while !Rig" becomes very large in the outer layer where the stratification is large compared to the mean shear. Because !Rig" $ 1 above the pycnocline in both cases, it is perhaps surprising that local occurrences of Rig # 0.25 are not uncommon even above the pycnocline. Figure 13b shows the probability of the local Rig # 0.25. In the mixed layer (where the distance from the pycnocline is large and negative), the probability of Rig # 0.25 is very high, as expected. The probability of

Rig # 0.25 drops to nearly zero in the pycnocline, but a local maximum in the probability occurs above the pycnocline. A coincident local maximum in the probability of local density overturns can also be seen above the pycnocline (not shown). About 70% of the overturns occurring above the pycnocline are associated with a local shear that is larger in magnitude than the planeaveraged value. This suggests that while the mean shear is not unstable based on the typical gradient Richardson number criteria, local shear instabilities may drive the overturns observed in this region. To illustrate the appearance of overturns and unstable shear profiles in the region above the pycnocline, Fig. 14 shows an instantaneous x%z slice through a small section of the flow when N&/f ' 75. The shading indicates du(/dz, lines show isotherms at an interval of 0.01)d!*"/dz&, and circles indicate locations where Rig # 0.25. The perturbation shear appears to be closely associated with undulations in the pycnocline height. As we have seen, occurrences of Rig # 0.25 are common above and below the pycnocline. Most of the regions with an unstable shear above the pycnocline are

NOVEMBER 2008

TAYLOR AND SARKAR

2547

FIG. 13. (a) Mean gradient Richardson number, &Rig', and (b) probability of occurrences of the local Rig " 0.25. Vertical profiles have been averaged in terms of the distance from the maximum temperature gradient.

associated with a negative streamwise shear perturbation. It also appears that radiated internal waves (visible by phase lines of du!/dz that slope up and to the left) are sometimes associated with Rig " 0.25, but this only occurs near the pycnocline; in the outer layer, Rig is always large. Because changes in the local shear and the local stratification can cause variations in Rig, it is of interest

to examine the distribution of low Rig events. Figure 15 shows a scatterplot of the deviation of the local shear and buoyancy frequency from the background values, for events with 0 " Rig " 0.25. Only one height is shown for clarity, z/# $ 0.195, which corresponds to the secondary peak in Rig above the pynocline, as shown in Fig. 13b. At this location, the mean buoyancy frequency is 73f, and the mean gradient Richardson number is

FIG. 14. Instantaneous streamwise shear, du!/dz, with overlaid isopycnals for N%/f $ 75. Circles show where Rig " 0.25 locally.

2548

JOURNAL OF PHYSICAL OCEANOGRAPHY

FIG. 15. Perturbation streamwise shear and buoyancy frequency for events with 0 $ Rig $ 0.25 at z/) ! 0.195. The mean gradient Richardson number, "Rig# ! 3.15.

Rig ! "N2#/"S#2 ! 3.15. Because the mean shear is dominated by the streamwise component at this location, only this component is considered. It is clear from Fig. 15 that most of the occurrences of Rig $ 0.25 are when du%/dz $ 0 and nearly all occurrences of Rig $ 0.25 are when the stratification is less than its mean value. Few exceptions occur when the shear is very large and negative, with Rig $ 0.25 even though the local stratification is high.

6. Turbulence-generated internal waves In the outer layer above the pycnocline, vertically propagating internal waves can be observed. These waves, which are generated as turbulent eddies interact with the stratified ambient, were also observed at a low Reynolds number and are described in detail by Taylor and Sarkar (2007a). The importance of the energy radiated by turbulence-generated internal waves to the mixed layer growth has been discussed in several previous studies (Linden 1975; E and Hopfinger 1986). An

A&kh, !, z' ! A0&kh, !, z0'

!

examination of the steady-state turbulent kinetic energy equation, integrated through the boundary layer, reveals that the turbulent production can be balanced by three sink terms: the integrated dissipation rate, the integrated buoyancy flux, and the vertical energy flux associated with the internal wave field. To compare the relative sizes of the three sink terms, Fig. 16 shows the vertical energy flux, " p%w%#, normalized by the integrated dissipation and the integrated buoyancy flux. To leading order, the turbulent production is balanced by the dissipation, indicating that the bulk mixing efficiency is very small. This is expected given that most of the dissipation occurs in the unstratified region near the seafloor. Of the two remaining sink terms, the vertical energy flux is comparable to the integrated buoyancy flux. Taylor and Sarkar (2007a) made a similar analysis of the boundary layer energetics for a lower Reynolds number boundary layer and also found that the vertical energy flux was much smaller than the integrated dissipation but of the same order as the integrated buoyancy flux. Specifically, they found that the ratio of the vertical energy flux to the integrated dissipation at the top of the pynocline was between 0.01 and 0.03, while the ratio of the vertical energy flux to the integrated buoyancy flux at the same location was between 0.5 and 1. The present results are generally consistent with those of Taylor and Sarkar (2007a), but the vertical energy flux is about a factor of 2 smaller in the present study. Taylor and Sarkar (2007a) found that the internal wave energy spectrum in the outer layer could be explained by viscous damping of the waves based on the wavenumber and vertical propagation speed for a particular frequency. Given a turbulent spectrum that was assumed to be characteristic of the waves upon generation, they showed that the viscous decay term in the fully nonlinear numerical simulations resulted in outer layer waves in a relatively narrow frequency range. Viscous ray tracing was used to show that the vertical velocity amplitude for a specific frequency and wavenumber can be written as

("! | k 0| exp | k| kh&!2 ( f 2'1#2

where k0 and A0 are the wavenumber and wave amplitude at z ! z0 in the generation region. Because the molecular viscosity and diffusivity used by Taylor and Sarkar were two orders of magnitude smaller than the

VOLUME 38

"

z

z0

#

| k|4&N2 ( !2'(1#2 dz$ ,

&21'

present values, it is of interest to evaluate their viscous internal wave model using the present high Reynolds number simulations. The viscous decay model is compared to the observed frequency spectra of the outer

NOVEMBER 2008

2549

TAYLOR AND SARKAR

FIG. 16. Vertical energy flux normalized by (a) the integrated turbulent dissipation and (b) the integrated buoyancy flux.

layer waves in Fig. 17. To compare the predicted wave amplitude to the simulations, it is convenient to consider the spectral distribution of dw/dz, which is estimated by multiplying A(kh, !, z) by the vertical wavenumber m " #kH

!

N2 # !2

!2 # f 2

"

1"2

.

$22%

Because the subgrid-scale eddy viscosity that is used as part of the LES model is not negligible compared to the molecular viscosity in the outer layer, we have used & ' &sgs in Eq. (21). Figure 17 shows the spectral amplitudes of (w)/(z normalized by * and u . To show the com* bined contributions of all values of kH, the square root of the sum of the squared amplitudes of (w)/(z is shown as a function of !/f. The left-hand side shows the observed spectra at z " z0, corresponding to a location just above the pycnocline (solid line). The spectrum at z " z0 was then smoothed and used as input (dashed line) to the viscous decay model. The viscous decay model predicts the amplitude of each frequency and wavenumber component after propagating a distance z # z0. The predicted amplitude from the model is compared to the observed amplitudes from the numerical simulation on the right-hand side, corresponding to a location at the top of the computational domain. The

qualitative agreement between the observed and predicted wave amplitudes is good; in particular, the decrease in amplitude of the low-frequency waves is captured well. At high frequencies, as ! → N, the observed wave amplitudes are significantly lower than the model prediction. It is possible that this is the result of nonlinear wave–wave and wave–mean flow interactions that are neglected in the linear viscous model.

7. Evaluating methods for estimating the wall stress Because we have high-resolution velocity and density profiles through a steady Ekman layer, we are able to evaluate the performance of several methods for estimating the friction velocity from observational data under idealized conditions. The profile method estimates the friction velocity by assuming that the mean shear follows the unstratified law-of-the-wall, to give u

*,p

" #z

d+u, . dz

$23%

When the turbulent dissipation is measured, the friction velocity can be estimated by assuming a balance between the turbulent production and dissipation. As shown in the inset of Fig. 11, the production and dissi-

2550

JOURNAL OF PHYSICAL OCEANOGRAPHY

VOLUME 38

FIG. 17. Spectra of *w+/*z from simulated results (solid line) and the viscous decay model (dashed line) for (a),(b) N, f ! 31.6 and (c),(d) N, f ! 75. Vertical lines show N/f.

pation are in approximate balance throughout the mixed layer. The so-called balance method assumes a production and dissipation balance using the observed mean shear and disspation rate, so that the friction velocity is u ,b ! *

!" !

d"u# . dz

$24%

A third method called the dissipation method can be formed through a combination of the profile and balance methods by using the production and dissipation balance and assuming that the mean shear follows the unstratified law-of-the-wall. The friction velocity predicted by this method takes the form u

*,!

! "!"z#1#3.

$25%

Estimates from these models are compared to the observed friction velocity in Fig. 18. To illustrate the temporal variability in the friction velocity, &1' (the standard deviation) is also shown. As has been found by previous studies (Perlin et al. 2005; Johnson et al. 1994; Lien and Sanford 2004), the performance of the friction velocity estimates depends on the location of the observed shear and on dissipation. Therefore, estimates of the friction velocity are shown using the mean shear and/or the dissipation rate at two heights, z/( ! 0.06 and z/( ! 0.12, which correspond to about 1.98 and 3.96 mab, respectively. When the outer flow is unstratified,

as shown in the left column of Fig. 18, all of the above methods provide a reasonable estimate of the friction velocity. However, when the outer layer is stratified, the profile method and to a lesser degree the balance method do not provide accurate estimates for the friction velocity. The dissipation method appears to be the most accurate of the three methods in the mixed layer. However, because direct measurements of the dissipation rate are not always available from observations, it is desirable to have a method for estimating the friction velocity using more commonly measured quantities such as velocity and density profiles. The error in the profile method is the result of the increase in the mean shear with stratification as was seen in Fig. 9, which is not accounted for by the traditional law-of-the-wall. A modified law-of-the-wall that accounts for the increase in mean shear at the top of a stratified boundary layer was proposed by Perlin et al. (2005). They proposed that the mean velocity gradient could be modeled as d"u# u* , ! dz lp

$26%

lp ! "z$1 ) z#hd%,

$27%

where

and hd is a measure of the boundary layer depth, which is limited by stratification. To use Eq. (26) to predict the friction velocity for a given velocity profile, an es-

NOVEMBER 2008

2551

TAYLOR AND SARKAR

FIG. 18. Estimates of the friction velocity using several different methods at two locations in the mixed layer. Horizontal lines show the friction velocity observed in the simulations and .1/ of the time series. (left to right) Models are the balance method Eq. (24), the dissipation method Eq. (25), the profile method Eq. (23), the modified law-of-the-wall Eq. (26), and the modified profile method Eq. (30). Note that when the flow is unstratified, the modified law-of-the-wall and the modified profile method are identical to the profile method.

timate for hd must first be obtained. Perlin et al. (2005) proposed using the Ozmidov scale, lOz ! "#/N3, at the top of the mixed layer to set hd. We have used a similar criterion to evaluate this method in Fig. 18. Specifically the mixed layer height d is defined as the location where $%/d ! 0.01d%/dz&. Then hd is set so that lp ! lOz at z ! d. It is evident from Fig. 18 that the modified law-ofthe-wall provides a significant improvement over the profile method. However, as shown in Fig. 19, the Ozmidov scale varies very rapidly between the mixed layer and the pycnocline, so the estimate of hd depends strongly on the definition of the mixed layer depth. In the case when N&/f ! 31.6, the rapid change in the Ozmidov scale leads to hd ! 0.67', which is significantly larger than the boundary layer height, h ! 0.215'. As a result, lp is significantly larger than the observed shear length scale in the upper half of the mixed layer, as shown in Fig. 20. An alternative method to account for the decrease in the length scale with stratification was proposed by Brost and Wyngaard (1978) and can be written as 1 1 1 ! ( , l !z lb

length scale is shown in Fig. 20 and compares favorably to the observed shear length scale below the center of the pycnocline. Practically, it is difficult to measure the vertical velocity, especially in the boundary layer where it cannot be deduced from isopycnal displacements. We have found that most of the decrease in the length scale with height is due to an increase in the local stratification rather than to a change in the turbulent velocity. An alternative length scale can then be formed by re-

)28*

where lb ! Cb+w ,w ,-1/2/N is a buoyancy length scale. Nieuwstadt (1984) suggested the value of Cb ! 1.69, which was consistent with his local scaling theory. This

FIG. 19. Ozmidov scale.

2552

JOURNAL OF PHYSICAL OCEANOGRAPHY

VOLUME 38

FIG. 20. Length-scale profile derived from the mean shear from the LES (line with filled circles) compared to several model profiles. Shaded regions show where d&p'/dz ) d&p'/dz(.

placing the vertical turbulent velocity with the friction velocity: 1 1 N#z$ . ! " l !z Cbu *

#29$

This simplified form still provides a reasonable estimate for the shear length scale, as shown in Fig. 20. The friction velocity can then be recovered from the mean velocity and density profile without the need for direct turbulence measurements, specifically u

*,m%p

! !z

!" # $ # $ % d&u' dz

2

"

d&"' dz

2

1#2

&

% N#z$ ,

&30'

where we have taken Cb ! 1. In the limit of an unstratified boundary layer, this method becomes equivalent to the profile method, so we will refer to this as the modified profile (m-p) method. Estimates of the friction velocity based on Eq. (30) are shown as triangles in Fig. 18. While the estimated friction velocity is somewhat large, the modified profile method provides a significant improvement over the profile method and is comparable to the modified law-of-the-wall. It is worth noting that because stratification effects enter into the modified law-of-the-wall through hd in Eq. (27), which is independent of z, stratification effects are nonlocal in this model. By comparison, the modified profile

method is the result of a local balance between shear and stratification. The length scale given in Eq. (29) can also be used to form a model velocity profile by integrating dU u* ! , dz l to obtain 1 U#z$ 1 ! log#z#z0$ " u ! Cbu * *

#31$

'

z

N#z$$ dz$.

#32$

0

If a constant buoyancy frequency, say, N(, were used in place of N(z) in Eq. (32), the form of the velocity profile would become log linear. The model velocity in Eq. (32) can, therefore, be viewed as an analog to the loglinear profile from Monin–Obukhov theory that is commonly used in the atmospheric literature, with the Obukhov length replaced by u /N. Figure 21 shows * profiles of the observed mean velocity compared to Eq. (32) for both a constant and a nonconstant N. To smooth fluctuations in the instantaneous profile of N(z), the profiles of N(z) have been averaged over a time window of t ! 1/f. When a depth-dependent buoyancy frequency, N(z), is used in Eq. (32), the parameter is set to Cb ! 1 to be consistent with Eq. (30). When N(z) is replaced by the constant free-stream buoyancy frequency, N(, it is found that Cb ! 0.2 provides a good

NOVEMBER 2008

TAYLOR AND SARKAR

2553

FIG. 21. Simulated velocity profiles compared to predictions by the model profiles of Eq. (32) for a constant and nonconstant N.

fit to the observed profiles. Using either a constant or a nonconstant buoyancy frequency in Eq. (32), the deviations from the logarithmic law owing to stratification are well represented. Note that in the outer layer, the velocity profiles are affected by the boundary layer height and Eq. (32) is not expected to hold.

8. Conclusions We have examined a benthic Ekman layer formed when a uniformly stratified, steady geostrophic flow encounters a flat, adiabatic seafloor. The thermal field rapidly develops a three-layer structure with a wellmixed region near the wall separated from the uniformly stratified outer layer by a pycnocline. The outer layer is populated by upward-propagating internal waves that are generated by the boundary layer turbulence. After the initial spinup, a quasi-steady state is reached, characterized by a slow mixed layer growth and a nearly constant density gradient in the pycnocline. When the strength of the outer layer stratification is increased, the wall stress increases slightly, but the boundary layer thickness decreases significantly. The structure of the boundary layer is clearly confined by stratification as evidenced by the Reynolds stress, turbulent heat flux, and turning angle, which all nearly vanish above the pycnocline. Because the Ekman transport balances the wall stress in the integrated streamwise momentum equation, and the wall stress remains

relatively constant, the increase in outer layer stratification is accompanied by an increase in the magnitude of the cross-stream velocity. The increased cross-stream velocity in the mixed layer leads to a broadening of the Ekman spiral. The rate of veering in the mixed layer is a function of z but does not depend strongly on the external stratification. When the stratification is increased and the boundary layer is thinner and the surface turning angle is larger, the amount of veering that occurs in the pycnocline increases. This finding is consistent with the results of Weatherly and Martin (1978), who found that most of the veering occurred in the pycnocline using onedimensional simulations. We have also found that when the outer layer stratification is large, the rapid rate of turning in the pycnocline causes the mean velocity and the mean shear to be maximum at the same location near the center of the pycnocline. An interesting feature of the observed Ekman layer structure is the appearance of local shear instabilities above the pycnocline, despite the fact that the mean shear is stable with respect to the local stratification. Between 10% and 25% of the vertical profiles exhibit a local gradient Richardson number less than 0.25. A similar observation has been made in terms of the occurrence of overturning events. Mixing during these events is significant and appears to cause a local minimum in the density gradient above the pycnocline. Analysis of events with Rig ! 0.25 at a height above the

2554

JOURNAL OF PHYSICAL OCEANOGRAPHY

pycnocline indicates that these events are more likely to be associated with a below average density gradient than an above average shear. In the outer layer, turbulence-generated internal waves are observed radiating away from the boundary layer. The vertical energy flux associated with these waves is negligible compared to the turbulent dissipation integrated through the boundary layer. This is not surprising because most of the dissipation occurs near the seafloor where the flow is unstratified. However, the vertical energy flux at the top of the boundary layer is nearly half of the integrated buoyancy flux. This finding is consistent with a study at a lower Reynolds number by Taylor and Sarkar (2007a) and implies that the turbulence-generated internal waves may remove enough energy from the boundary layer to affect the growth rate of the mixed layer. The viscous internal wave model of Taylor and Sarkar (2007a) has been applied using the combined molecular and turbulent viscosity to estimate the decay rate of the turbulencegenerated internal waves. The model qualitatively captures the decay of low-frequency waves but overestimates the amplitude of high-frequency waves. As was seen in previous studies (e.g., Perlin et al. 2005; Johnson et al. 1994), an increase in the mean shear has been observed at the top of the mixed layer. The increase in mean shear with respect to an unstratified boundary layer can lead to significant errors in the friction velocity estimated from observed velocity profiles using the profile method. It has been uncertain whether this increase in the mean shear could be explained in terms of the local stratification. Because the unstratified logarithmic law appears to hold very close to the wall, the profile method is adequate, in principle, if the mean velocity very near the wall can be obtained. However, this is often difficult or impossible in practice. We have evaluated the performance of a variety of techniques for estimating the friction velocity given mean quantities at various heights in the mixed layer. Because the turbulent production and dissipation are the dominant terms in the turbulent kinetic energy equation throughout the mixed layer, the dissipation method agrees very well with the observed friction velocity. The modified law-of-the-wall, proposed by Perlin et al. (2005), shows considerable improvement over the standard profile method, especially near the top of the mixed layer. However, like the balance and dissipation methods, the modified law-of-the-wall requires knowledge of the turbulent dissipation rate. Direct observation of the dissipation rate is difficult, especially in active turbulent regions because it involves small-scale velocity gradients. When such information is not available, it is desirable to have an alternative method for

VOLUME 38

estimating the friction velocity u . We have introduced * a buoyancy length scale, u /N(z), which leads to a * modification of the well-known log-law by a buoyancyrelated augmentation of the mean velocity. Use of this modified mean velocity profile instead of the unstratified log law leads to a significant improvement in deducing the friction velocity. Acknowledgments. We are grateful for the support provided by Grant N00014-05-1-0334 from ONR Physical Oceanography, program manager Scott Harper. REFERENCES Armenio, V., and S. Sarkar, 2002: An investigation of stably stratified turbulent channel flow using large-eddy simulation. J. Fluid Mech., 459, 1–42. Armi, L., 1978: Some evidence for boundary mixing in the deep ocean. J. Geophys. Res., 83, 1971–1979. Bewely, T., 2008: Numerical Renaissance: Simulation, Optimization, and Control. Renaissance Press, 547 pp. Brost, R., and J. Wyngaard, 1978: A model study of the stably stratified planetary boundary layer. J. Atmos. Sci., 35, 1427– 1440. Caldwell, D., and T. Chriss, 1979: The viscous sublayer at the sea floor. Science, 205, 1131–1132. Coleman, G., 1999: Similarity statistics from a direct numerical simulation of the neutrally stratified planetary boundary layer. J. Atmos. Sci., 56, 891–900. ——, J. Ferziger, and P. Spalart, 1990: A numerical study of the turbulent Ekman layer. J. Fluid Mech., 213, 313–348. ——, ——, and ——, 1992: Direct simulation of the stably stratified turbulent Ekman layer. J. Fluid Mech., 244, 677–712. Dewey, R., P. Leblond, and W. Crawford, 1988: The turbulent bottom boundary layer and its influence on local dynamics over the continental shelf. Dyn. Atmos. Oceans, 12, 143–172. Durran, D., 1999: Numerical Methods for Wave Equations in Geophysical Fluid Dynamics. Springer-Verlag, 465 pp. Garrett, C., P. MacCready, and P. Rhines, 1993: Boundary mixing and arrested Ekman layers: Rotating stratified flow near a sloping boundary. Annu. Rev. Fluid Mech., 25, 291–323. Germano, M., U. Piomelli, P. Moin, and W. Cabot, 1991: A dynamic subgrid-scale eddy viscosity model. Phys. Fluids A, 3, 1760–1765. Haidvogel, D., and A. Beckmann, 1999: Numerical Ocean Circulation Modeling. Imperial College Press, 320 pp. Howard, L., 1961: Note on a paper of John W. Miles. J. Fluid Mech., 10, 509–512. Johnson, G., T. Sanford, and M. Baringer, 1994: Stress on the Mediterranean outflow plume. Part I: Velocity and water property measurements. J. Phys. Oceanogr., 24, 2072–2083. Lele, S., 1992: Compact finite difference schemes with spectrallike resolution. J. Comput. Phys., 103, 16–42. Lien, R.-C., and T. Sanford, 2004: Turbulence spectra and local similarity scaling in a strongly stratified oceanic bottom boundary layer. Cont. Shelf Res., 24, 375–392. Linden, P., 1975: The deepening of a mixed layer in a stratified fluid. J. Fluid Mech., 71, 385–405. Marusic, I., G. Kunkel, and F. Porte-Agel, 2001: Experimental study of wall boundary conditions for large-eddy simulation. J. Fluid Mech., 446, 309–320.

NOVEMBER 2008

TAYLOR AND SARKAR

Miles, J., 1961: On the stability of heterogeneous shear flows. J. Fluid Mech., 10, 496–508. Munk, W., 1966: Abyssal recipes. Deep-Sea Res., 13, 707–730. Nieuwstadt, F., 1984: The turbulent structure of the stable, nocturnal boundary layer. J. Atmos. Sci., 41, 2202–2216. Orszag, S., 1971: Numerical simulation of incompressible flow within simple boundaries. I. galerikin (spectral) representation. Stud. Appl. Math., 50, 293. Perlin, A., J. N. Moum, J. M. Klymak, M. D. Levine, T. Boyd, and P. M. Kosro, 2005: A modified law-of-the-wall applied to oceanic bottom boundary layers. J. Geophys. Res., 110, C10S10, doi:10.1029/2004JC002310. ——, ——, ——, ——, ——, and ——, 2007: Organization of stratification, turbulence, and veering in bottom Ekman layer s . J . Geophys . Re s., 1 1 2 , C 05 S9 0, d o i : 1 0 . 1 0 2 9 / 2004JC002641. Pope, S. B., 2000: Turbulent Flows. Cambridge University Press, 810 pp. Shingai, K., and H. Kawamura, 2002: Direct numerical simulation of turbulent heat transfer in the stably stratified Ekman layer. Thermal Sci. Eng., 10, 25. Stahr, F., and T. Sanford, 1999: Transport and bottom boundary layer observations of the North Atlantic Deep Western Boundary Current at the Blake Outer Ridge. Deep-Sea Res. II, 46, 205–243. Stoll, R., and F. Porte-Agel, 2006: Effect of roughness on surface boundary conditions for large-eddy simulation. Bound.Layer Meteor., 118, 169–187.

2555

Taylor, J., and S. Sarkar, 2007a: Internal gravity waves generated by a turbulent bottom Ekman layer. J. Fluid Mech., 590, 331– 354. ——, and ——, 2007b: Near-wall modeling for large-eddy simulation of an oceanic bottom boundary layer. Proc. Fifth Int. Symp. on Environmental Hydraulics, Tempe, AZ, Arizona State University, 6 pp. ——, and ——, 2008: Direct and large eddy simulations of a bottom Ekman layer under an external stratification. Int. J. Heat Fluid Flow, 29, 721–732. ——, ——, and V. Armenio, 2005: Large eddy simulation of stably stratified open channel flow. Phys. Fluids, 17, 116602. Thorpe, S., 2005: The Turbulent Ocean. Cambridge University Press, 439 pp. Weatherly, G., and P. Martin, 1978: On the structure and dynamics of the oceanic bottom boundary layer. J. Phys. Oceanogr., 8, 557–570. Xuequan, E., and E. J. Hopfinger, 1986: On mixing across an interface in stably stratified fluid. J. Fluid Mech., 166, 227– 244. Zilitinkevich, S., and A. Baklanov, 2002: Calculation of the height of the stable boundary layer in practical applications. Bound.Layer Meteor., 105, 389–409. ——, and I. Esau, 2002: On integral measures of the neutral barotropic planetary boundary layer. Bound.-Layer Meteor., 104, 371–379.