Abstract

The first part of this paper is generic; it demonstrates a problem associated with one-dimensional, ocean surface layer model comparisons with ocean observations. Unlike three-dimensional simulations or the real ocean, kinetic energy can inexorably build up in one-dimensional simulations, which artificially enhances mixing. Adding a sink term to the momentum equations counteracts this behavior. The sink term is a surrogate for energy divergence available to three-dimensional models but not to one-dimensional models.

The remainder of the paper deals with the Mellor–Yamada boundary layer model. There exists prior evidence that the model’s summertime surface temperatures are too warm due to overly shallow mixed layer depths. If one adds a sink term to approximate three-dimensional model behavior, the warming problem is exacerbated, creating added incentive to seek an appropriate model change. Guided by laboratory data, a Richardson-number-dependent dissipation is introduced and this simple modification yields a favorable improvement in the comparison of model calculations with data even with the momentum sink term in place.

1. Introduction

This paper is in two parts. The first part, sections 2 and 3, is generic. We show that numerical one-dimensional simulations of mixed layer behavior should include a sink parameterization to account for energy flux divergence that is generally present in the ocean and in three-dimensional, numerical ocean simulations. Without an energy sink, wind-driven, surface layer velocities generated by a one-dimensional model will inexorably increase and, in the course of a multimonth simulation, erroneously enhance mixed layer deepening and surface temperature cooling. One presumes that this has introduced error in past comparisons of one-dimensional model results with observational data. Introducing a momentum sink term is a completely empirical means of simulating three-dimensional energy divergence, but is deemed better than excluding the sink. And Pollard and Millard (1970) did incorporate a sink term in a simple model in order to match the model with current meter observations. Since their paper, the emphasis has been predominantly on the prediction of mixed layer temperatures and salinities wherein currents were ignored.

Whereas sections 2 and 3 should apply to most models, the second part of the paper focuses on the Mellor-Yamada (henceforth M–Y) model, which is often used in atmospheric and ocean general circulation models. It is not a planetary boundary layer model per se; it is a more general turbulence closure model that happened to ape quite a few properties of stratified planetary boundary layers (Mellor 1973; Mellor and Yamada 1974, 1982). Other second moment models with application to geophysical flows include those by Lewellen and Teske (1973), Lumley et al. (1978), Andre and Lacarrere (1985), Therry and Lacarrere (1983), and Rodi (1987). The model of Andre and Lacarrere is labeled a third moment model by the authors, since closures are constructed for the fourth moment terms in an equation for the diffusional third moments. However, the more important pressure–velocity gradient terms and the dissipation related to the gradient of the two-point triple correlation are still closed at a lower order so that the model is predominantly a second moment model. For a more general review of mixed layer models including eddy viscosity and bulk models, see Kantha and Clayson (1994).

By incorporating model hypotheses and nondimensional constants derived from neutral laboratory data, the M–Y model very nearly reproduced the Monin–Obukhov similarity relations for near-surface stratified flows, both stable and unstable. However, in application to the entire surface boundary layer, the model has been criticized for overly shallow oceanic surface layer depths and overly warm, sea surface temperatures during summertime stable conditions. In particular, Martin (1985) showed that summertime temperatures at Ocean Weather Station Papa exceeded observations by about 2°C relative to an annual range of 8°C. For this reason, Kantha and Clayson (1994) deemed it necessary to amend the M–Y model and directly add eddy viscosity and diffusivity (thus abandoning a modeling principal that empirical constants should be nondimensional) to the base of the mixed layer.

The findings in section 2 and 3 exacerbate the problem since, with the addition of a momentum sink, summertime surface temperatures are increased. Therefore a modification is made to the M–Y model and that is to replace the existing submodel for dissipation with a Richardson number (based on turbulence variables) dependent, model dissipation. The need for such a modification was first revealed by a laboratory experiment of decaying, density-stratified turbulence (Dickey and Mellor 1980, henceforth DM) wherein the turbulence at first decays in approximately classical fashion (turbulence kinetic energy ∝ time−1) but then, rather abruptly, transitions to a random ensemble of internal waves with very little decay, an equilibrium state similar to that occurring in the ocean (Garrett and Munk 1979). It is somewhat of a jump in parameter space from that experiment to the shear-dominated regimes of planetary boundary layers. Still, one would hope that the model would include the DM case. As we will demonstrate, inclusion of a Richardson-number-dependent, energy dissipation increases the turbulence kinetic energy in the outer portion of the layer, increases boundary layer thickness, and improves comparison with observations using wind stress and heat flux compiled by Martin (1985).

There are, of course, other sources of disageement between model calculations and one-dimensional data. Aside from oversimplication of model physics and data errors, they include model finite difference errors, unknown vertical and horizontal advection of observed mean properties, interaction with unknown remotely generated internal and surface waves or with baroclinic velocity shears, and errors or lack of resolution in the surface forcing data that drive a model.

2. Artificial velocity buildup in one-dimensional ocean surface layer models

Pollard and Millard (1970) adopted a simple ocean mixed layer model:

 
formula

wherein the vertical divergence of the wind stress, (∂τzx/∂z,τzy/∂z) = (τox, τoy)/h, was distributed uniformly over a layer of fixed depth, h, and the velocity components, u and υ, are independent of depth. According to Pollard and Millard,

The linear damping term models the (three-dimensional) dispersion effect by introducing a decay factor of the form, exp(−ct); c−1 is the e-folding time. It should be noted that the damping term has the effect of shifting the dominant frequency response of the model to a frequency slightly less than f when the frequency spectrum of the forcing is flat. This is not a good property when modeling the real ocean where inertio–gravity waves have frequencies between f and the Brunt–Väisälä frequency N. Since N is usually greater than f, the dominant frequency observed is usually slightly greater than f.

Figure 1 is a comparison of the simple model for c−1 = 8 inertial days and h = 45 m, with observational data. The point to be made is that this empirical damping term was considered important toward improved comparison with current meter measurements. Other cases were considered in the paper and various values of c were chosen to minimize model error.

Fig. 1.

Wind record and response of a simple slab model (heavy line) by Pollard and Millard (1970) and compared with current meter measurements (light line) at 39°20.5′N, 69°59′W and a depth of 7 m during Oct 1965

Fig. 1.

Wind record and response of a simple slab model (heavy line) by Pollard and Millard (1970) and compared with current meter measurements (light line) at 39°20.5′N, 69°59′W and a depth of 7 m during Oct 1965

We wish now to derive somewhat more general results related to one-dimensional, ocean, surface layer models. To simplify the matter, consider the vertical integrals of (1a,b), which, for zero stress at depth, are

 
formula

where (Sx, Sy) ≡ (0hu dz, 0hυ dz) and h, a constant, is large enough so that (τhx, τhy) = (0, 0). The conventional Ekman transport is obtained when c = 0 and the tendency terms are nil. Equations (2a,b) are an exact, zero-dimensional model, which conveniently excludes consideration of detailed profile structure.

We have available the 3-hourly wind data used by Martin (1985) at OWS Papa (50°N, 145°W) and November (30°N, 140°W). These winds are used to numerically compute (S2x, + S2y)1/2 from (2a,b) and the results are plotted as solid lines in Figs. 2a,b for c = 0 and c = (10 d)−1.

Fig. 2.

The growth of Ekman transport for (a) OWS Papa and (b) OWS November. The solid lines are numerical solutions of Eqs. (2a,b) where the dashed lines are obtained from (6) and (7). In each diagram the upper curves are for c = 0 whereas the lower curves are for c = (10 d)−1

Fig. 2.

The growth of Ekman transport for (a) OWS Papa and (b) OWS November. The solid lines are numerical solutions of Eqs. (2a,b) where the dashed lines are obtained from (6) and (7). In each diagram the upper curves are for c = 0 whereas the lower curves are for c = (10 d)−1

Signals from the c = 0 case increase with time and this result can also be obtained from integrals of velocity profiles calculated by the full z-dependent model (see section 4) if we set h ≥ 160 m, the depth of maximum penetration of turbulence kinetic energy, Reynolds stress, and currents. On the other hand, the signal from the c ≠ case flattens, on average, to constant values. We now show that this behavior is to be expected for systems such as (2a,b).

3. Statistical expectation

For the present purpose, the winds at Papa and November can be approximated as a random time signal after the means of (Sx, Sy) = f−1(τox, −τoy) have been subtracted. Without change in nomenclature and after subtracting the means, (2a,b) may be written

 
formula

where στox + oy and SSx + iSy and i ≡ (−1)1/2. A solution to (3) for resting initial conditions is

 
formula

a. Expected behavior with no decay

In Eq. (4), let c = 0. If now we write the modulus squared of S, we obtain

 
formula

where the overbars denote ensemble means and an asterisk denotes a complex conjugate. Now let t2t + τ/2, t1tτ/2, and σ(t1)σ∗(t2) = Rσ(τ) = Rσ(−τ). After some manipulation (Taylor 1921), we obtain

 
formula

For small T compared to the decorrelation time, SS = T2Rσ(0), whereas for large T

 
formula

Wind stress correlations have been obtained from the yearlong records for stations Papa and November (courtesy P. Martin) and are plotted in Fig. 3. An approximation to Rσ(τ) is Rσ(0)ea|τ| after which the quadrature in (6) may be evaluated such that 0Rσ(τ) cos = Rσ(0)a/(a2 + f2). This result plus the values, Rσ(0) and a, obtained from Fig. 3, provide the dashed lines in Figs. 2a and 2b wherein (SS)1/2T1/2. Note that one should ideally compare the average of an ensemble of data with the statistical model. Nevertheless, there is little doubt that the statistical model explains the behavior of the two experiments in that a one-dimensional Ekman layer velocities, forced by a wind, approximated by a random process, will inexorably increase.

Fig. 3.

The wind stress correlation functions for OWS Papa and November. The units of Rσ are m4 s−4

Fig. 3.

The wind stress correlation functions for OWS Papa and November. The units of Rσ are m4 s−4

b. Expected behavior with decay

We now repeat the analysis for c ≠ 0 and obtain

 
formula

for large T instead of (6). With the exponential approximation for Rσ(τ) noted above, 0Rσ(τ)e cos = Rσ(0)(a + c)/[(a + c)2 + f2]. Since, according to Pollard and Millard, ca, the ect term in (7) can be replaced by unity. Thus one can find the expected constant response for large time by replacing 2T with c−1 in (6); the dashed constant values of (SS)1/2 in Figs. 2a,b may be easily obtained for a specific value of c. The illustrations are for c = (10 d)−1.

c. One-dimensional behavior as a function of decay rate

In the preceding analysis, we have shown that, statistically, a wind-driven, one-dimensional, ocean surface layer model will accumulate increasingly larger Ekman transport and, therefore, velocity and kinetic energy. This should lead to enhanced mixing and mixed layer deepening relative to the same model imbedded in a three-dimensional, ocean model where horizontal inhomogeneity is enabled.

In Fig. 4, we show the sensitivity of monthly averaged, SST on the decay rate. We use the full M–Y model (described below). It would therefore seem that, in one-dimensional simulations, a momentum decay term as in (1a,b)—or, perhaps, a better substitute yet to be devised—is necessary, especially for long-term model runs. Whereas decay rates have been previously used, they have not been employed in most one-dimensional simulations where the focus has been on simulating mixed layer deepening and temperature and salinity distributions.

Fig. 4.

Sensitivity of the behavior of monthly averaged, OWS Papa, model SST to the decay rate, c. The time units are inertial days. The circles are data

Fig. 4.

Sensitivity of the behavior of monthly averaged, OWS Papa, model SST to the decay rate, c. The time units are inertial days. The circles are data

4. The Mellor–Yamada model

Anticipating the need to correct the M–Y model, a short summary of that model is now provided. A complete derivation, justification of model constants, and model data comparisons are to be found in Mellor and Yamada (1982) and some numerical details are in Mellor (1996).

We introduced the model some years ago (Mellor 1973; Mellor and Yamada 1974; Mellor and Durbin 1975). The model demonstrated skill in solving diverse turbulence problems; these include homogeneous decaying isotropic and anisotropic turbulence, neutral boundary layer, channel flow, pipe flow, and other neutral flows such as separating flows (Celenligil and Mellor 1985) and flows over curved walls (which can suppress or enhance turbulence much as density stratification does; Mellor 1975) and stratifed flows such as the Monin-Obukhov similarity flow variables are predicted.

The M–Y model consists of the following:

  1. A hypothesis by Rotta (1951a,b) for the pressure, velocity gradient covariance term and extended to the pressure, temperature (density) gradient covariance terms. The hypothesis relates these terms to linear functions of the Reynolds stress and heat (density) flux and are tensorially unambiguous.

  2. The Kolmogorov (1941) hypothesis for the dissipation extended to other dissipation terms involving the fluctuating temperature (density) gradients. The choice here is unambiguous.

  3. Fickian type gradient terms for turbulence diffusion terms. The choice here is somewhat ambiguous but these terms play a relatively minor role (they are excluded from the level 2 model).

  4. Steps 1 and 2 lead to the definition of four length scales and a nondimensional constant. It is then assumed that all are proportional to one another or to a “master length scale.” There are then five constants to be determined. These constants are unambiguously related to measured laboratory turbulence data for neutral flows.

  5. A process of formula simplification (Mellor and Yamada 1974) based on assumed small (but not zero) departures from isotropy. This led to a hierarchy of model versions, labeled levels 1, 2, 3, and 4. In the level 4 version, one would need to prognostically solve for all components of the Reynolds stress and heat flux tensors and other variances, thus the need for simplication.

There is left the specification of the master length scale, l. But assuming that l = kz = 0.4z is valid near a solid surface (the five constants were scaled as a group so that l became asymptotically equal to Prandtl’s mixing length as a solid surface is approached), the Monin–Obukhov similarity variables were derived and compared with near-surface wind data obtained from a meteorological tower in Kansas (Businger et al. 1971) as shown in Fig. 5. It was this result (Mellor 1973) that encouraged further development and application of the model. Apparently, not all second moment models have reproduced this result.

Fig. 5.

Monin–Obukhov similarity variables where, appropriate to the atmosphere, ϕMκz(∂U/∂z)/uτ, ϕMκzuτ(∂ϒ/∂z)/H, and ζκzgβH/u3τ. The distance from the surface is z, κ is the von Kármán constant, uτ is the friction velocity, H is the kinematic surface heat flux, β is the thermal expansion coefficient, g is the gravity constant, ∂U/∂z is the velocity gradient, ∂ϒ/∂z is the potential temperature gradient. These quantities can all be shown to be functions of GH using (9a,b,c) and the level 2 approximation (Mellor and Yamada 1982); the functions provide the solid curves in the figures. The data are from Businger et al. (1971) 

Fig. 5.

Monin–Obukhov similarity variables where, appropriate to the atmosphere, ϕMκz(∂U/∂z)/uτ, ϕMκzuτ(∂ϒ/∂z)/H, and ζκzgβH/u3τ. The distance from the surface is z, κ is the von Kármán constant, uτ is the friction velocity, H is the kinematic surface heat flux, β is the thermal expansion coefficient, g is the gravity constant, ∂U/∂z is the velocity gradient, ∂ϒ/∂z is the potential temperature gradient. These quantities can all be shown to be functions of GH using (9a,b,c) and the level 2 approximation (Mellor and Yamada 1982); the functions provide the solid curves in the figures. The data are from Businger et al. (1971) 

The so-called level 2½ models, arising out of step 5 and the boundary layer approximation, lead to expressions for the vertical diffusivities, KM and KH, such that

 
formula

where q2/2 is the turbulence kinetic energy and l is the turbulence master scale. The coefficients, SM and SH, are functions of a Richardson number given by

 
formula

where

 
formula

is a Richardson number. The factor, ρ̃/∂z, is the vertical density gradient minus the adiabatic lapse rate. The five constants in (9a,b) were evaluated from near-surface turbulence data (law-of-the-wall region) and decaying homogeneous turbulence data; they were found (Mellor and Yamada 1982) to be (A1, B1, A2, B2, C1) = (0.92, 16.6, 0.74, 10.1, 0.08). The stability functions, SM(GH) and SH(GH), are plotted in Fig. 6. The stability functions limit to infinity as GH approaches the value 0.0288. Absent discretization error, model flows cannot exceed this value since stratification would have been destroyed by indefinitely large KM and KH.

Fig. 6.

The stability factors SM and SH as obtained from Eqs. (2a,b)

Fig. 6.

The stability factors SM and SH as obtained from Eqs. (2a,b)

Since the modeling hypotheses and nondimensional constants were all based on neutral data, the extension to stratified flows is almost entirely a derived result.

In Mellor and Yamada (1982), SM and SH were also listed as additional functions of either the parameter GMl2[(∂U/∂z)2 + (∂V/∂z)2]/q2 or the parameter (Ps + Pb)/ɛ, where Ps, Pb, and ɛ are shear production, buoyancy production, and dissipation respectively. Subsequently, Yamada (1983) found that (Ps + Pb)/ɛ could be set to unity without much change in calculated results. Then Galperin et al. (1988) justified this step by recognizing that (Ps + Pb)/ɛ = 1 + O(a2), where a2 is a nondimensional measure of the departure from isotropy (Mellor and Yamada 1974). Within the rules of the level 2½ model, the O(a2) term could be neglected. In practice, we have also found negligible differences in performance between the two versions and the newer version does avoid some numerical complexities. However, there is a disadvantage: when one examines the resulting turbulence components, they do not add up to q2; the difference is q2O(a2).

In the level 2½ version of the model, a prognostic equation is solved for q2 (twice the turbulent kinetic energy equation). If we define the material derivative, Df/Dt = ∂f/∂t + ∂(Uf)/∂x + ∂(Vf)/∂y + ∂(Wf)/∂z, then, the equation is

 
formula

The terms on the right are the vertical turbulence diffusion, the shear production, the buoyancy production, the dissipation, and horizontal diffusion. The dissipation is speciifically labeled for later amendment.

Large eddy simulations by Moeng and Wyngard (1986, 1989) of convectively driven boundary layers have demonstrated the role of buoyancy terms in addition to that which has entered into (9a,b) and (10) so that we may add such terms to the model in the future as has been done by Sun and Ogura (1980) in their applications of the M–Y level 2½ model; Lumley et al. (1978) and Therry and Lacarrere (1986) have also included such terms in their models. However, modeling the mean variables in a convective layer is particularly easy: we do not expect that these terms would have a large net effect on the mean variables and, in this paper, we would rather deal solely with that which we deem to be a more important model feature.

Left for last is a prescription for the master length scale, l, which entails a greater degree of empiricism than what has preceeded. We have used an equation for q2l (the master length scale equation, which is related to integrals of the two-point, correlation equations). Thus,

 
formula

in which ≡ 1 + E2(l/κL), where L−1 ≡ |z|−1 + |D + z|−1 is a “wall proximity” function and 0 < z < −D;Kq is the vertical turbulence diffusivity; and (E1, E2, E3) = (1.8, 1.33, 1.0) are additional nondimensional constants (see Mellor and Yamada 1982). Boundary conditions for (10) and (11) are that ∂q2/∂z = l = 0 on bounding surfaces. It can be shown, from (10) and the law of the wall, that ∂q2/∂z = 0 is formally equivalent to the alternate condition, q2 = B2/31u2τ; here, uτ is the friction velocity, defined as the square root of the surface stress divided by density, and B1 is one of the model constants. For alternative boundary conditions, which are meant to account for breaking waves, see Craig and Banner (1994) and Stacey and Pond (1997).

Since Eq. (11) is the most empirical part of the model, it has often been replaced by an algebraic equation for l, namely, l/lo = kz/(kz + lo) where loα ∫ zq dz/∫ q dz. The early value of α = 0.1 (Mellor and Yamada 1974) has been revised to α = 0.2 (Mofjeld and Lavelle 1984; Martin 1985). The later values also agree more nearly with calculated results using (11). We generally use (11) instead of the algebraic relation because it seems to yield credible results for multiple turbulent regions for, say, an ocean with surface and bottom boundary layers separated by an inviscid region or for merged layers. It also expresses the fact that eddy scales are transportable quantities; that is, they have memory, an application of which will be discussed in the next section.

Solutions to (10) and (11) together with (8) and (9) yield values of KM and KH, that can be used to solve for the horizontal velocity components, potential temperature, and salinity. Thus,

 
formula

where f is the Coriolis parameter and FU, FV, Fϒ, FS are horizontal diffusion terms. The solar radiative flux divergence is given by ∂R/∂z. The understanding is that, in three-dimensional calculations, c = 0, whereas in one-dimensional calculations, c ≠ 0 and Df/Dt reduces to ∂f/∂t.

Notice that, if one forms a kinetic energy equation from (12a,b), energy can be lost (or perhaps gained) in three-dimensional flow by the divergence of the rate of work term,  · (pU), and tranfer to potential energy via the term, Wρg, (Pollard 1970); whereas, in the one-dimensional world, it is accomplished by the term, c(U2 + V2).

5. A Richardson-number-dependent dissipation submodel

In this section, we set the stage for a correction to the Mellor–Yamada model. The correction is suggested by an experiment by Dickey and Mellor (1980, henceforth DM). In the experiment, a square grid of rods of mesh size 2.54 cm was towed vertically upward in a tank (0.66 × 0.66 m cross section and 2.44 m high) of, first, fresh unstratified water and, second, water in which a vertical salinity stratification had been created such that the Brunt–Väisälä frequency N = 0.378 rad s−1. The result is shown in Fig. 7. In the unstratified experiment, the grid turbulence decayed like the classical wind tunnel experiments of Batchelor and Townsend (1948) and others (see additional references in DM), such that q2t−1 [data depart somewhat from the t−1 behavior for large t as explained by Domaradski and Mellor (1984)]. In the stratified experiment, the turbulence decayed as in the unstratified case until a critical Richardson number was reached; whence the decay process nearly ceased and did so rather abruptly. What had been nonlinearly interacting, turbulent eddies became nearly linear internal waves and q2 was approximately constant. According to DM, the small decay that remains is due to moloecular viscosity and would decrease as the Reynolds number increased.

Fig 7. Decay of turbulence for neutral, ○, and stratified, ‖ cases; in the latter, ∂ρ/∂z = 1.46 × 10−4 g cm−4. Curve 1 is the solution to Eqs. (14a), (15), and (16a,b,c). Curve 2 is the best fit to the data (the variance of the data around the straight line fit is explained in DM). Curve 3 is the initial period decay law. For remarks on curve 4, see Domaradzki and Mellor (1984). The grid speed Wg = 95.0 cm s−1, the mesh M = 5.08 cm, and the grid Reynolds number = 48 260

Fig 7. Decay of turbulence for neutral, ○, and stratified, ‖ cases; in the latter, ∂ρ/∂z = 1.46 × 10−4 g cm−4. Curve 1 is the solution to Eqs. (14a), (15), and (16a,b,c). Curve 2 is the best fit to the data (the variance of the data around the straight line fit is explained in DM). Curve 3 is the initial period decay law. For remarks on curve 4, see Domaradzki and Mellor (1984). The grid speed Wg = 95.0 cm s−1, the mesh M = 5.08 cm, and the grid Reynolds number = 48 260

The governing equation for decaying, homogeneous turbulence is

 
formula

The dissipation is modeled such that

 
formula

to which we add a model equation for the macro-dissipation scale

 
formula

Note that (15) derives from (11) for this simple flow case since Λ = B1l. Also DM showed that Λ is, not surprisingly, of the order of the integral macroscale (Batchelor 1956). From (14a,b) and (15), we recover the initial period decay behavior such that q2t−1 and Λ ∝ t1/2; this simulated the DM unstratified flow case and other classical data.

To simulate the stratified case and the cutoff in the decay rate, the following model was proposed for ɛ:

 
formula

where

 
formula

The power 3/2 behavior was suggested in the DM paper, but it is not critical to the main purpose of Eq. (16), which is to provide a cutoff for ɛ around RqRqc; the empirical need for the cutoff was unambiguous, and to fit the DM turbulence decay data, it was determined that Rqc ≅ 100 or ΛN/q ≅ 10. The result is shown in Fig. 7. Not shown is the fact that Λ first increases as t1/2 and then becomes constant.

The buoyancy production term in the DM experiment was previously assumed to be negligibly small. This assumption has been questioned by a reviewer and upon review of DM and Dickey (1977) we find that the reviewer is correct; the buoyancy production is generaly small and not important to the main results in DM, but it is not negligible. Therefore we provide appendix A wherein the issue is discussed and the buoyancy production is estimated by the M–Y model. It should be noted that the level 2½ or 3 model does not apply to the DM experiment since dissipation and production are very much out of balance in the turbulence energy equation, one must appeal to the level 4 model.

6. Station Papa and November data and a M–Y model correction

The work by Martin (1985) indicated that the M–Y model produced overly shallow, stable summertime surface layer depths and consequently overly warm SSTs. With the decay terms in (12a,b) in place, the M–Y model produces even warmer summertime SSTs as shown in Fig. 4 so that there is need for a model change.

We follow the prescription in section 5 with a change in nomenclature to conform to the M–Y model. Thus, Rq = −B21GH such that

 
formula

where GHc is a constant yet to be specified. The only difference in (16a–c) and (17a,b), other than the nomenclature change, is the introduction of the value 0.1, instead of 0.0, for GHGHc. This is cosmetically motivated. During model execution, where negative GH is large, q2 is very small, and l is physically ill-defined; then the addition of the small value, 0.1, creates smooth distribution of l(z) that, otherwise, is noisy.

Henceforth, Eqs. (17a,b) is substituted for (10b). This is equivalent to a modification of B1 in (10b) to B1 = B1;im(GH)−1. Consistant with the development of the level 2½ model, this B1 modification should also be applied to (9a,b) and also to B2 = B2;im(GH)−1 so that it is assumed that the dissipation of density variance (see appendix A) is similarly modified.

Problematic is the actual choice of the decay constant in the momentum equation when dealing with one-dimensional simulations. Based on the paper of Pollard and Millard (1970), we hereby choose the value c = (8 inertial days)−1 for OWS Papa and November and for the data discussed in sectons 7 and 8. With c fixed, we will determine the best value of GHc to fit the data previously invoked by Martin (1985). If and when a better estimate of c is available, GHc may have to be amended. The parameters of the calculation are displayed in Table 1. As generally required by the leapfrog temporal discretization, a temporal smoother (Asselin 1972) is used in the model. However, the value, 0.02, can be shown to correspond to a decay constant of ∼(100 d)−1 and is small enough as to not compete with the explicit decay constant, c = (8 inertial days)−1, introduced into (12a,b). However, it does have an effect relative to cases where c = 0; for example, the buildup of transport in Fig. 2 is decreased somewhat if we apply an Asselin smoother to that calculation.

Table 1.

Characteristics of the 3-hourly surface forcing (Martin 1985) and the numerical parameters

Characteristics of the 3-hourly surface forcing (Martin 1985) and the numerical parameters
Characteristics of the 3-hourly surface forcing (Martin 1985) and the numerical parameters

For c = 0, Fig. 4 includes the case for GHc = −∞ such that ;im(GH) = 1 for OWS Papa. Figures 8 and 9 are the results for c = (8 inertial days)−1 for Papa and November, respectively. From these results we conclude that GHc ≅ − 6.0. This translates to the value ΛN/q ≅ 40 and is larger than the value appropriate to the DM experiment. This introduces the possibility that ;im = ;im(GH, GM).

Fig. 8.

Monthly averaged OWS Papa data (○) and model calculation for c = (8 inertial days)−1 and various values of GHc. The value, GHc = −∞ represents the uncorrected M–Y result

Fig. 8.

Monthly averaged OWS Papa data (○) and model calculation for c = (8 inertial days)−1 and various values of GHc. The value, GHc = −∞ represents the uncorrected M–Y result

Fig. 9.

As in Fig. 9 but for OWS November

Fig. 9.

As in Fig. 9 but for OWS November

Figure 10 shows sample two-week plots of the model velocity components at a depth of 8 m for the three cases: (a) no momentum decay and the unrevised M–Y model, (b) momentum decay and the unrevised M–Y model, and (c) momentum decay and the revised M–Y model. In Fig. 10a the currents are quite large.

Fig. 10.

Sample time series of the model velocity components for OWS Papa and at a depth of 8 m for (a) c = 0, GHc = −∞; (b) c = (8 d)−1, GHc = −∞; and (c) c = (8 d)−1, GHc = −6.0. The units are centimeters per second

Fig. 10.

Sample time series of the model velocity components for OWS Papa and at a depth of 8 m for (a) c = 0, GHc = −∞; (b) c = (8 d)−1, GHc = −∞; and (c) c = (8 d)−1, GHc = −6.0. The units are centimeters per second

To cover the entire year, we next plot time–depth contours of current speed in Fig. 11 for the different conditions cited in Fig. 10; current speed unto itself filters out the mainly rotary inertial oscillations. Figure 12 shows plots of the square root of twice the turbulence energy, q, and Fig. 13 shows plots of temperature. Figure 13c compares favorably with OWS Papa data, which we reproduce here in Fig. 14 from Martin’s paper. Notice evidence of upwelling; this could be added to the model with some guesswork as to the vertical velocity distribution. The calculated data have been averaged over 5 days before plotting. It is clear that, in comparing the (b) plots with the (a) plots, the turbulence and mixing is reduced due to reduced velocity shear, resulting in a shallower and warmer mixed layer in the summer and fall. The revised model, shown in the (c) plots, exhibits increased turbulence in the outer portion of the mixed layer; this is, of course, due to the reduced dissipation. The increased turbulence and mixing is also responsible for reduction in velocity shear.

Fig. 11.

Current speed contours for the same parameters as in Fig. 10. The contour interval is 0.05 m s−1

Fig. 11.

Current speed contours for the same parameters as in Fig. 10. The contour interval is 0.05 m s−1

Fig. 12.

Contours of q, the square root of twice the turbulence energy, for the same parameters as in Fig. 10. The contour interval is 0.01 m s−1

Fig. 12.

Contours of q, the square root of twice the turbulence energy, for the same parameters as in Fig. 10. The contour interval is 0.01 m s−1

Fig. 13.

Temperature contours for the same parameters as in Fig. 10. The contour interval is 1°C

Fig. 13.

Temperature contours for the same parameters as in Fig. 10. The contour interval is 1°C

Fig. 14.

Temperature contours for 1961 from OWS Papa for the year 1961. Compare with the model contours of Fig. 13c (after Martin 1985)

Fig. 14.

Temperature contours for 1961 from OWS Papa for the year 1961. Compare with the model contours of Fig. 13c (after Martin 1985)

It should be clear that GHc will henceforth be a model closure constant; it is rooted in laboratory data but, at present, it must be regarded as a “tuning constant” somewhat reluctantly admitted to the M–Y model. On the other hand, c is not a closure constant; it is empirically introduced to compensate for missing processes in a one-dimensional model formulation necessary for comparison with one-dimensional data. When the closure model is incorporated into three-dimensional ocean models, one should set c = 0; our expectation is that three-dimensional performance will be improved using (17a,b).

7. The LOTUS dataset

The Long-Term Upper Ocean Study (LOTUS) was conducted during 1982 and 1983 at a site in the Sargasso Sea (34°N, 70°W) as described in papers by Stramma et al. (1986) and Price et al. (1987). The data of July 1982 has been compared with models by Stramma et al., Gaspar et al. (1990), Large et al. (1994), and Kantha and Clayson (1994).

The surface forcing data are plotted in Fig. 15a. The solar insolation was directly measured during the LOTUS experiment together with wind velocity and air and sea surface temperatures. Longwave radiation, sensible and latent heat flux, and the “remainder heat flux” in Fig. 15a, were calculated using standard formulas. The ocean surface albedo was set at 0.05. The cloud cover was fixed at 0.3, but tests showed negligible sensitivity of the longwave radiation to cloud cover. Following Stramma et al. (1986), the relative humidity was fixed at 75%. There was a sensor deployed on the hull of the attendant meteorological buoy, which measured temperature at a depth of 0.6 m (and taken equal to the SST mentioned above); the sensor recorded a strong diurnal signal during times of light wind. The period, starting at midday of 12 July 1999 and extending for the next 7–10 days was selected variously by the previous authors as being particularly free of advective effects. We have chosen to extend the study period to 28 days in order to illustrate the degree to which advective processes can affect model data comparisons.

Fig. 15.

The LOTUS data: (a) Surface heat flux (W m−2); (b) Temperature depth contours (CI = 1°C); (c) Heat storage (°C m) from surface heat flux (dashed curve) and from depth integrals of temperature profiles (solid curve)

Fig. 15.

The LOTUS data: (a) Surface heat flux (W m−2); (b) Temperature depth contours (CI = 1°C); (c) Heat storage (°C m) from surface heat flux (dashed curve) and from depth integrals of temperature profiles (solid curve)

The measured temperature contours, smoothed with a 24-h boxcar filter are shown in Fig. 15b. The heat storages, one obtained from the meteorological forcing and the other by integrating the temperature data down to 50 m, are displayed in Fig. 15c. Although previous authors have cited unimportant imbalances of the two heat storages in the range of Julian day 195–205, the imbalances of Fig. 15c in this range, presumably related to advection, do seem significant. After day 205, the imbalance changes sign and is quite large.

The model was initialized at day 194.5 with measured temperatures and currents, thus avoiding a rapidly developing imbalance in the preceeding day or two. The model was then executed for 28 days. Figures 16a and 16b show the 0.6-m temperatures—exhibiting a strong diurnal signal—and the water column temperature contours, respectively. The near-surface diurnal excursions of the model are sensitive to the radiative penetrative algorithm as noted by Stramma et al. (1986), which motivated them to adopt the two component radiation formulas of Paulson and Simpson (1977) as derived from Jerlov (1976); the present model incorporates these formulas.

Fig. 16.

(a) Surface temperatures from LOTUS (solid curve) and from the model (dashed curve). (b) The full depth–time temperature contours: CI = 1°C

Fig. 16.

(a) Surface temperatures from LOTUS (solid curve) and from the model (dashed curve). (b) The full depth–time temperature contours: CI = 1°C

The peak temperatures are higher than obtained by Large et al. (1994) and Kantha and Clayson (1994), but are about the same as those obtained by Stamma et al. (1986). The present model peak temperatures are also higher than the data in the range 195 to 204 days but lower in the range 208 to 222 days. However, this mismatch is well correlated with the imbalance of the heat storage obtained from the surface heat flux and that obtained from the measured temperature profiles. (Note that the two heat storage determinations are identical for the model.) The model and observed temperatures are nearly coincident when the imbalance is small in the range 204–208 days.

One might enlist the known imbalances to approximately correct the model, but we are not motivated to do so here.

8. The FLEX dataset

In the spring of 1976 the Fladen Ground Experiment (FLEX) was carried out in the North Sea at 58°55′N, 0°32′E (Soetje and Huber 1980) in a water depth of about 145 m. Hourly meteorological data was obtained including solar insolation and surface properties from which Burchard and Baumert (1995), using standard formulas, estimated wind stress and the remainder heat flux that comprised sensible and latent heat flux and downward longwave radiation. The heat flux record is shown in Fig. 17a. CTD data was analyzed into six hourly temperature profiles displayed in Fig. 17b; storms occurred on days 112, 133, and 153.

Fig. 17.

The FLEX data: (a) Surface heat flux (W m−2). (b) Temperature depth contours: CI = 0.1°C. (c) Heat storage (°C m) from surface heat flux (dashed curve) and from depth integrals of temperature profiles (solid curve)

Fig. 17.

The FLEX data: (a) Surface heat flux (W m−2). (b) Temperature depth contours: CI = 0.1°C. (c) Heat storage (°C m) from surface heat flux (dashed curve) and from depth integrals of temperature profiles (solid curve)

The model calculations are displayed in Fig. 18. Figure 18b compares well with Fig. 17b. However, the model surface temperatures are high by a fraction of a degree compared to the measurements. However, as in the LOTUS case, the difference is well correlated with the imbalance of the two heat storage determinations in Fig. 17c. Also measured was radiation penetration from which variations of the deeper pentration component of the radiation were incorporated into the model. The model results were not too sensitive to this detail compared to the use of Jerlov optical Type IIA constants.

Fig. 18.

(a) Surface (1.25 m) temperatures from FLEX (solid curve) and from the model (dashed curve). (b) The full depth–time temperature contours: CI = 0.1°C

Fig. 18.

(a) Surface (1.25 m) temperatures from FLEX (solid curve) and from the model (dashed curve). (b) The full depth–time temperature contours: CI = 0.1°C

9. Summary and discussion

In this paper, we first point out an error in comparing one-dimensional model simulations of ocean surface boundary layers with ocean data due to an unrealistically large accumulation of model kinetic energy. There are of course other sources of error but this represents an omnipresent bias. Given the rather large kinetic energy increase in one-dimensional calculations, it is surprising that the error in SST and concommitant mixed layer deepening after six months is not larger. Still, it is of the same order as the perceived summertime error in the M–Y model results obtained by Martin (1985) and in this paper and it is additive. Following Pollard and Millard (1970), we correct for the lack of three-dimensional energy dispersion through the use of Rayleigh damping in the momentum equations.

The core of the M–Y model is contained in the constants (A1, B1, A2, B2, C1) that connect the governing turbulence length scales and relate to known turbulence structure properties. In this paper, to correct for the warm bias in summertime SSTs in the model, we introduce a Richardson number cutoff to the turbulence dissipation. This is not an arbitrary decision; it is virtually dictated—at least qualitatively—by the DM experiment and, indeed, by observations of random internal waves in the ocean (Garrett and Munk 1979). While the principal of the cutoff is well founded, the constant, GHc, unlike the length scale constants, must be regarded as a tuning constant so that model calculations agree more closely with ocean surface data.

A separate finding, in the course of this research, was that B1 and B2 in (9a,b) may be kept as constants, in which case GHc = 2.5, and the effect on the model results is nearly identical to that obtained with Richardson-number-dependent B1 and B2 in (9a,b) and where GHc = −6.0. In either case, in the range of the data, the model results in Fig. 5 are not affected.

Although this paper has dealt with surface boundary layers, the model can and has been applied to bottom boundary layers.

Fig. A1. The terms in the energy budget from the model, whereby tendency = buoyancy production + dissipation. The dashed curve is the neutral case, where buoyancy production is nil. The abscissa has been offset by the value, 35, in accordance with the virtual origin of the data in Fig. 7 

Fig. A1. The terms in the energy budget from the model, whereby tendency = buoyancy production + dissipation. The dashed curve is the neutral case, where buoyancy production is nil. The abscissa has been offset by the value, 35, in accordance with the virtual origin of the data in Fig. 7 

Acknowledgments

This paper owes much to Paul Martin who supplied me with the OWS Papa data and with his code to calculate the surface fluxes. And his 1985 paper obviously triggered this paper. Jim Price helped me access the LOTUS dataset. Lakshmi Kantha and Boris Galperin contributed some valuable comments. This research was primarily funded by the Office of Naval Research under Grant. N00014-93-1-0037 and partially funded by the NOPP Coastal Marine Demonstration Project.

REFERENCES

REFERENCES
Andre, J. C., and P. Lacarrere, 1985: Mean and turbulent structure of the oceanic surface layer as determined from a one-dimensional third order simulation. J. Phys. Oceanogr., 15, 121–132.
Asselin, R., 1972: Frequency filters for time integration. Mon. Wea. Rev., 100, 487–490.
Batchelor, G. K., 1956: The Theory of Homogeneous Turbulence. Cambridge Press, 103 pp.
——, and A. A. Townsend: 1948: Decay of isotropic turbulence in the initial period. Proc. Roy. Soc., 193A, 539–566.
Burchard, H., and H. Baumert: 1995: On the performance of a mixed-layer model based on the k-e(greek epsilon) turbulence closure. J. Geophys. Res., 100, 8523–8540.
Businger, J. A., J. C. Wyngaard, Y. Izumi, and E. F. Bradley, 1971: Flux profile relationships in the atmospheric boundary layer. J. Atmos. Sci., 28, 181–189.
Celenligil, M. C., and G. L. Mellor, 1985: Numerical solutions of two-dimensional turbulent separated flows using a Reynolds stress closure model. J. Fluids Eng., 107, 467–476.
Craig, P. D., and M. L. Banner, 1994: Modeling wave-enhanced turbulence in the ocean surface layer. J. Phys. Oceanogr., 24, 2546–2559.
Dickey, T. D., 1977: An experimental study of decaying and diffusing turbulence in neutral and stratified fluids. Ph.D. dissertation, Princeton University, 133 pp. [Available from Xerox University Microfilms, 300 North Zeeb Road, Ann Arbor, MI 48106.]
.
——, and G. L. Mellor, 1980: Decaying turbulence in neutral and stratified fluids. J. Fluid Mech., 99, 13–31.
Domaradzki, J. A., and G. L. Mellor, 1984: A simple turbulence closure hypothesis for the triple-velocity correlation functions in homogeneous isotropic turbulence. J. Fluid Mech., 140, 45–61.
Galperin, B., L. H. Kantha, S. Hassid, and A. Rosati: 1988: A quasi-equilibrium turbulent energy model for geophysical flows. J. Atmos. Sci., 45, 55–62.
Garrett, C., and W. Munk, 1979: Internal waves in the ocean. Annu. Rev. Fluid Mech., 11, 339–369.
Gaspar, P., Y. Gregoris, and J.-M. Lefevre, 1990: A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station Papa and Long-Term Upper Ocean Study site. J. Geophys. Res., 95, 16 179–16 193.
Jerlov, N. G., 1976: Marine Optics 14. Elsevier Science, 231 pp.
Kantha, L. H., and C. A. Clayson, 1994: An improved mixed layer model for geophysical applications. J. Geophys. Res., 99, 25 235–25 266.
Kolmogorov, A. N., 1941: The local structure of turbulence in incompressible viscous fluid for very large Reynolds number (in Russian). Dokl. Akad. Nauk. SSSR, 30, 301–310.
Large, W. G., J. C. McWilliams, and S. C. Doney, 1994: Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. Rev. Geophys., 32, 363–403.
Lewellen, W. S., and M. E. Teske, 1973: Prediction of the Monin–Obukhov similarity functions from an invariant model of turbulence. J. Atmos. Sci., 30, 1340–1345.
Lumley, J. L., O. Zeman, and J. Seiss, 1978: The influence of buoyancy on turbulent transport. J. Fluid Mech., 84, 581–597.
Martin, P. J., 1985: Simulation of the mixed layer at OWS November and Papa with several models. J. Geophys. Res., 90, 903–916.
Mellor, G. L., 1973: Analytic prediction of the properties of stratified planetary boundary layers. J. Atmos. Sci., 30, 1061–1069.
——, 1975: A comparative study of curved flow and density stratified flow. J. Atmos. Sci., 32, 1278–1282.
——, 1996: User’s guide for a three-dimensional, primitive equation, numerical ocean model. Program in Atmospheric and Oceanic Sciences, Princeton University, 41 pp. [Available online at http://www.aos.princeton.edu/WWWPUBLIC/htdocs.pom/.]
.
——, and T. Yamada, 1974: An hierarchy of turbulence closure models for planetary boundary layers. J. Atmos. Sci., 13, 1791–1806.
——, and P. A. Durbin, 1975: The structure and dynamics of the ocean surface mixed layer. J. Phys. Oceanogr., 5, 718–728.
——, and T. Yamada, 1982: Development of a turbulent closure model for geophysical fluid problems. Rev. Geophys., 20, 851–875.
Moeng, C.-H., and J. C. Wyngaard, 1986: An analysis of closures for pressure–scalar covariances in the convective boundary layer. J. Atmos. Sci., 43, 2499–2513.
——, and ——, 1989: Evaluation of turbulent transport and dissipation closures in second-order modeling. J. Atmos. Sci., 46, 2311–2330.
Mofjeld, H. O., and J. W. Lavelle, 1984: Setting the length scale in a second order closure model of the unstratified bottom boundary layer. J. Phys. Oceanogr., 14, 833–839.
Paulson, C. A., and J. J. Simpson, 1977: Irradiance measurements in the upper ocean. J. Phys. Oceanogr., 7, 952–956.
Pollard, R. T., 1970: On the generation of inertial waves in the ocean. Deep-Sea Res., 17, 795–812.
——, and R. C. Millard, 1970: Comparison between observed and simulated wind-generated inertial oscillations. Deep-Sea Res., 17, 813–821.
Price, J. F., R. A. Weller, C. M. Bowers, and M. G. Briscoe, 1987: Diurnal response of sea surface temperature observed at the Long-Term Upper Ocean Study (34°N, 70°W) in the Sargasso Sea. J. Geophys. Res., 92, 14 480–14 490.
Rodi, W., 1987: Examples of calculation methods for flow and mixing in stratified flows. J. Geophys. Res., 92, 5305–5328.
Rotta, J. C., 1951a: Statistische Theorie nichthomogener Turbulenz. Z. Phys., 129, 547–572.
——, 1951b: Statistische Theorie nichthomogener Turbulenz. Z. Phys., 131, 51–77.
Stacey, M. W., and S. Pond, 1997: On the Mellor–Yamada turbulence closure scheme: The surface boundary conditions for q2. J. Phys. Oceanogr., 27, 2081–2086.
Stramma, L. P., P. Cornillon, R. A. Weller, J. F. Price, and M. G. Briscoe, 1986: Large diurnal sea surface temperature variability: Satellite and in situ measurements. J. Phys. Oceanogr., 16, 827–837.
Sun, W.-Y., and Y. Ogura, 1980: Modeling the evolution of the convective planetary boundary layer. J. Atmos. Sci., 37, 1558–1572.
Taylor, G. I., 1921: Diffusion by continuous movement. Proc. London Math. Soc., 20, 196–211.
Therry, G., and P. Lacarrere, 1986: Improving the eddy-kinetic energy model for boundary layer description. Bound.-Layer Meteor., 37, 129–148.
Yamada, T., 1983: Simulations of nocturnal drainage flows by a q2−1 turbulence closure model. J. Atmos. Sci., 40, 90–106.

APPENDIX

The Dickey–Mellor Experiment

A reviewer questioned neglect of the buoyancy production term, gwρ/ρ0, in the energy equation (14) and thus this appendix is meant to address his concern, but in so doing will add a new exercise of the M–Y model and some modeling information to the DM experiment.

Experimental data

In DM, it is simply stated that “it was possible to determine from mean salinity measurements that the buoyancy flux production term was negligible.” A review of DM and Dickey (1977) now reveals that this statement is erroneous; the reviewer was correct.

The density flux was not measured directly but an upper bound was estimated. In the experiment, the grid could be towed upward many times without affecting the salinity in the middle of the tank. This is due to the spatial homogeneity, that results in zero salinity flux divergence. The effect of vertical flux could only be seen at the ends of the tank where the flux itself must be nil. If we integrate the salinity equation from the middle of the tank, z = −H/2, to the top, z = 0, we obtain

 
formula

Note that ∂S/∂t was nonzero only for a short distance inward from the ends of the tank. Therefore the flux can be obtained from measurements of mean salinity. But salinity profiles could be measured only before and after each run so that only a time integral of ws(−H/2) was obtained. Furthermore, the mean salinity changes were not large, so there was error in measurement. Nevertheless, noting that the flux was only active over the first 3 s of the experiment—estimated from salinity fluctuation measurements in the center of the tank—Dickey (1977) estimated an upper bound that /ρ0 < 5 × 10−4 cm s−1. After multiplication by the gravity constant, one now finds that this value is not negligible. However, the buoyancy production is a turbulence energy sink term in a stably stratified environment and would tend to reduce the turbulence and not arrest the decay process as seen in the experiment. Also, the stratified data initially tracked the neutral data (Fig. 7) so that the buoyancy production could not have greatly affected the turbulence.

A model of the DM experiment

In view of the uncertainties concerning buoyancy production, we next seek model results. Because the tendency terms and dissipation presumably dominate the flow, one must turn to the full level 4 M–Y model. If we let θρ′/ρo be the turbulence density petrurbation (or the product of the salinity perturbation and the salinity coeficient of expansion) and w the vertical velocity peturbation, the relevant equation set is

 
formula

since spatial derivatives of mean quantifies are nil. The cutoff variable, ;im⁠, is provided by (16c) where Rqc = 100. Since q2t−1 and lt1/2 as t → 0 we can obtain the constants of proportionality from the DM experimental data for use as initial conditions. From the asymptotic behavior of (A1)–(A5) as t → 0, we find that

 
formula

an interesting result unto itself.

The equations were solved with simple one-step differencing; variables that appear as first and zero order derivatives in the same equation were treated explicitly;the other terms were explicit. Time steps were chosen successively smaller until convergence was assured at a time step of 0.01 s.

The energy budget is displayed in Fig. A1, where it is seen that the buoyancy production, while not negligible, is nonetheless not important to the decay process. Subsequent evalution of q−2 yields values close to the analytical results presented in Fig. 7. Although data for direct comparison with model buoyancy production is lacking, the model-derived role of the buoyancy production is consistant with the decay of the turbulence kinetic energy in the stratified case and the fact that the kinetic energy is close to that of the neutral case in the range, 0 < Wgt/M < 300.

For large Wgt/M, the model buoyancy production does exhibit a decaying oscillatory behavior (barely visible in Fig. A1) and this is probably unrealistic.

Footnotes

Corresponding author address: Dr. George L. Mellor, Program in Atmospheric and Oceanic Sciences, Princeton University, Box CN710, Sayre Hall, Princeton, NJ 08544-0710.