NOAA Technical Memorandum OAR PSD-313

IMPLEMENTATION OF A TWO-EQUATION VERTICAL TURBULENT MIXING 

SCHEME IN A MESOSCALE ATMOSPHERIC MODEL

J.-W. Bao

S. A. Michelson

L. H. Kantha

J. M. Brown

Earth System Research Laboratory 

Physical Sciences Division 

Boulder, Colorado

August 2008

NATIONAL OCEANIC AND noaa ATMOSPHERIC ADMINISTRATION

Office of Oceanic and

Atmospheric Research /

NOAA Technical Memorandum OAR PSD-313

IMPLEMENTATION OF A TWO-EQUATION VERTICAL TURBULENT MIXING 

SCHEME IN A MESOSCALE ATMOSPHERIC MODEL

J.-W. Bao

NOAA Earth System Research Laboratory, Boulder, CO

S. A. Michelson

NOAA Earth System Research Laboratory, Boulder, CO

CIRES, University of Colorado, Boulder, CO

L. H. Kantha

PAOS and Department of Aerospace Engineering Sciences, 

University of Colorado, Boulder, CO

J. M. Brown

NOAA Earth System Research Laboratory, Boulder, CO

Earth System Research Laboratory 

Physical Sciences Division

Boulder, Colorado 

August 2008

UNITED STATES

DEPARTMENT OF COMMERCE

Office of Oceanic and

Atmospheric Research

Dr. Richard Spinrad

Assistant Administrator

NATIONAL OCEANIC AND

ATMOSPHERIC ADMINISTRATION

VADM Conrad C. Lautenbacher, Jr.

Under Secretary for Oceans

and Atmosphere/Administrator

Carlos M. Gutierrez

Secretary
ii

NOTICE

Mention of a commercial company or product does not constitute an endorsement by 

NOAA/Earth System Research Laboratory. Use of information from this publication concerning 

proprietary products or the test of such products for publicity or advertising purposes is not 

authorized.

______________________________________________________

For sale by the National Technical Information Service, 5285 Port Royal Road

Springfield, VA 22061
iii

CONTENTS

ABSTRACT ....................................................................................................................................v

1. Introduction ............................................................................................................................... 1

2. The Vertical Turbulent Mixing Scheme ................................................................................... 3

3. The Closure Constants of the q

2

l Equation .............................................................................. 7

4. Constraints on the Length Scale .............................................................................................. 11

5. Numerical Integration ............................................................................................................ 15

6. Numerical Experiments .......................................................................................................... 17

7. Summary and Discussion ........................................................................................................ 27

8. References ............................................................................................................................... 29
iv
v

ABSTRACT

A two-equation vertical turbulent mixing scheme is implemented in the National 

Center for Atmospheric Research and the Pennsylvania State University mesoscale model 

(MM5-V3). This scheme is based on the Mellor-Yamada level 2.5 second-moment closure 

(MY closure), and consists of two prognostic equations: one for the turbulent kinetic energy 

(TKE) and the other for the length scale multiplied by twice the TKE. In this first effort to apply 

the scheme to simulating vertical turbulent mixing in the atmospheric boundary layer (ABL), it is 

found that the physical meaning of the length scale is of fundamental importance in the 

determination of the closure constants in the length-scale equation for simulating vertical 

turbulent mixing in the ABL. If the length scale is defined as the characteristic length scale of 

the largest energy-containing eddies and is related to the distance that these eddies travel in the 

vertical direction before losing their initial TKE due to turbulent mixing and buoyancy effects, it 

is concluded that the closure constants in the length-scale equation should be different than those 

proposed previously for oceanic applications. To ensure physically sensible performance of the 

scheme, the technique developed by Janjić (2002) is applied for deriving necessary constraints 

on the length scale and for numerically integrating the scheme. The constraints on the length 

scale are derived by requiring that the TKE equation be nonsingular under different stability 

regimes in terms of the gradient Richardson number. The numerical integration is performed 

using an innovative splitting algorithm to control the computational modes encountered when 

using conventional numerical schemes. The results from a series of numerical experiments 

indicate that when properly choosing these constants, the evolution of the ABL structure 

simulated by the scheme is similar to the original scheme of the MY closure in MM5-V3 where 

the length scale is diagnosed. Despite the theoretical appeal and the successful implementation 

of this two-equation scheme, future effort will be required to use comprehensive turbulence 

observations to further evaluate its performance.

1

1. Introduction

One of the popular methods to simulate vertical turbulent mixing in the atmospheric 

boundary layer (ABL) is to use the Mellor-Yamada level 2.5 second-moment closure 

(MY closure) turbulence model (Mellor and Yamada 1974, 1982), which is centered on the 

prognostic equation for the turbulent kinetic energy (TKE, denoted as q

2

/2). An important 

closure parameter in the turbulence model is the so-called length scale (denoted as l). Although 

the length scale is essential to the level 2.5 second-moment closure, there is no unique approach 

to its specification. The simplest approach to specifying the length scale is to assume that l is 

proportional to the distance from the surface within the ABL (see, e.g., Blackadar 1962, Mellor 

and Yamada 1982, and Janjić 1994). While this method works well in many applications, it is 

diagnostic and its treatment of the length scale far above the surface has a rather empirical nature 

(Holtslag 2003). Therefore, it is more desirable to use a prognostic equation with minimum 

empiricism to predict the length scale. A vertical turbulent mixing scheme of the MY closure 

uses a prognostic equation for specifying l is often called a two-equation scheme, in contrast with 

a one-equation scheme in which l is diagnosed.

Various approaches to deriving the prognostic equation for l have been proposed. Unlike 

the TKE equation that is derived from the Reynolds-averaged Navier-Stokes equations, the 

equation for l cannot be derived systematically from the governing laws of fluid motion. 

Essentially, it is cast as a transport equation for a quantity involving l that is analogous to the 

TKE equation with a different set of closure constants. In general, the equation for l can be 

written as a transport equation for a quantity involving q and l: c

p

qml

n

, where c is a constant, and

p, m and n are integers (see, e.g., Burchard 2002 and Kantha 2004). It has been suggested that 

for a homogeneous shear layer all the approaches to deriving the equation for l are equivalent 

(Burchard 2002, Umlauf and Burchard 2003, and Kantha 2004), and thus no one approach is 

fundamentally superior to another. It should be noted, however, that the ABL is not a 

homogeneous shear layer, and the non-uniform thermal stratification in the ABL plays an 

important role in controlling vertical turbulent mixing. Therefore, various length-scale equations 

may perform differently for the ABL.

In this paper, we present the implementation of a two-equation vertical turbulent mixing 

scheme in the National Center for Atmospheric Research and the Pennsylvania State University 

mesoscale model (MM5-V3). The prognostic equation for l in this scheme is formulated in 
2

terms of the turbulence transport of twice the TKE. This two-equation scheme was first 

proposed by Mellor and Yamada (1974, 1982) for general applications in geophysical flows, and 

later revised by Kantha and Clayson (1994) and Burchard (2001) to simulate the vertical 

turbulent mixing in the ocean. Although this scheme has been applied extensively in simulating 

oceanic mixed layer, this is the first effort to implement it in an atmospheric mesoscale model. 

Three major issues are encountered during the implementation. The first issue is related 

to the closure constants in the length-scale equation. The values of the closure constants in the 

length-scale equation are determined a priori according to turbulence measurements or direct 

turbulence simulations. In oceanic applications, measurements of simple turbulent flows 

generated in laboratories are commonly used to determine these closure constants (see, e.g., 

Mellor and Yamada 1982, Umlauf and Burchard 2003). It has been recently recognized that 

there exist constraining relationships among the closure constants that can be utilized to 

determine the closure constants according to comprehensive turbulence measurements, (see, e.g., 

Burchard 2002). Although in many oceanic applications the closure constants determined from 

the laboratory measurements have appeared to work well, uncertainties remain in the constants’ 

values and continuing efforts are still being made to refine them. It is also not clear whether the 

closure constants that have worked well in oceanic applications are well suitable for simulating 

vertical turbulent mixing in the ABL. Therefore, how to properly determine the closure 

constants that are suitable for applications in the ABL needs to be investigated.

The second issue is related how to ensure physically sensible performance of the scheme. 

Vertical turbulent mixing schemes based on the MY closure have been known to require the use 

of the so-called realizability constraints on the shear and stability parameters as well as on the 

length scale to perform reasonably and robustly (see, e.g., Mellor and Yamada 1982, Galperin 

et al. 1988, Helfand and Labraga 1988, and Burchard 2002). The equation for the length scale 

itself does not provide a guarantee that the predicted l is always in agreement with the 

realizability constraints such that the modeled turbulent mixing is physically meaningful for a 

given flow situation. Therefore, properly established constraints are required to limit the 

solution of the length-scale equation. Although various approaches to imposing the constraints 

on the length scale have been proposed for oceanic applications, they are not as theoretically 

appealing as that proposed by Janjić (2002). Janjić’s approach is innovatively different than the 

others in that the constraints on the length scale are derived based on requiring that the analytical 
3

solution of the TKE equation be nonsingular, and thus there is no need to constrain the shear and 

stability parameters like the other approaches. However, this approach has not been applied to a 

two-equation scheme.

The third issue is concerned with numerical integration of the scheme. In oceanic 

applications, the scheme is discretized implicitly in time with centered in space differencing for 

numerical integration. It is found that the numerical solution of the equations for both the TKE 

and the length scale are apparently contaminated with computational modes. To make the results 

from the two-equation scheme physically meaningful, the computational modes need to be 

eliminated or significantly reduced. The numerical scheme proposed by Janjić (2002) is found to 

be effective to that end. 

The rest of the paper is organized as follows. The two-equation vertical turbulent 

mixing scheme is described in section 2. Section 3 presents results from the investigation on 

how to determine the closure constants in the length-scale equation. How to apply Janjić’s 

approach to establishing the constraints on the length scale and to numerical integrating the twoequation scheme is presented respectively in sections 4 and 5. Results from a series of numerical 

experiments are shown in section 6 to illustrate the sensitivity of the scheme to different values 

of the closure constants within their admissible ranges, which is followed by summary and final 

remarks (section 7).

2. The Vertical Turbulent Mixing Scheme

The vertical turbulent mixing scheme consists of two equations: the equation for the 

TKE and the equation for l. The TKE equation is written as 

 = + − ε



 





































∂

∂

















∂

∂ − ∂

















∂

q Ps Pb

q

z

l q S

t z

q

2

2 2

2

, (1)

where

















∂

∂ −















∂

∂ = − zV

wv

z

U P wu s , Pb = βgwθv , (2)
4

B l

q

1

3

ε = , (3)

z

U

wu Km ∂

∂ − = , z

V

wv Km ∂

∂ − = , z

w K v

v H ∂

∂Θ − θ = , (4)

K m = lqSM , H H K = lqS , (5)

( ) ( ) 1 2 2 2 1 2 2 S M 6 A A GM + S H 1− 3 A B GH −12 A A GH = A (6a)

( ) ( ) ( ) 1 2 1 1

2

1 2 1

2

1 S M 1+ 6 A GM − 9A A GH − S H 12A GH + 9A A GH = A 1− 3C (6b)

































∂

∂

+ 















∂

∂















 =

2 2

2

2

z

V

z

U

q

l GM , z

g

q

l G v H ∂

∂Θ















 = − β 2

2

, (7)

where U , V and Θv are the velocity components and virtual potential temperature of the mean 

flow, and the rest of symbols are of the same meaning as in Mellor and Yamada (1982). 

The success of (1)-(7) in the simulation of vertical turbulent mixing is critically 

dependent on well-calibrated closure constants( , , , , , ) A1 A2 B1 B2 C1 Sq and an appropriate 

specification of the length scale. The values of the constants in (1)-(6) are chosen as those by 

Janjić (2002), i.e.,

(A1 , A2 , B1 , B2 , C1 , Sq , β ) =

(0.65988838, 0.65742096, 11.877992, 7.226971, 0.00083092297, 0.2, 1./270).

Note that they are different than those proposed in Mellor and Yamada (1982) that are also used 

in Kantha and Clayson (1994). They are chosen here because they are used in the original 

vertical turbulent mixing scheme (hereafter referred to as the ETA scheme) in the mesoscale 

model in which this two-equation scheme is implemented. This original scheme, developed by 

Janjić (1990, 2002) for the National Centers for Environmental Prediction’s (NCEP) operational 

limited area model, is the same as this two-equation scheme except that the length scale is 

diagnosed. Using these values will make it convenient and fair to compare the results from this 
5

two-equation scheme and those from the original one-equation scheme (see the detail of the 

comparison later in section 6). 

Following Mellor and Yamada (1982) and Kantha and Clayson (1994), the equation 

for l is expressed as a transport equation of q

2

l, i.e., 

( ) ( ) l Psl Pbl l q l

z

qlS

z

q l

t = + − ε 



 





∂

∂

∂

∂ − ∂∂ 2 2 , (8)

where

sl Ps P = lE1 , bl Pb P = lE2 , B l

q

l F

1

3

ε = . 

In the above equation, E1 , E2 , and F are the closure constants. Note that in both Mellor and 

Yamada (1982) and Kantha and Clayson (1994), F is not a constant, but a wall function that 

approaches a constant far from the surface. It seems that the inclusion of the wall function is not 

absolutely necessary if E1 is chosen properly. As discussed later in section 4, if E1 and E2 are 

chosen to be constants, F should be a constant according to the validity of (8) in a neutral 

surface layer. Rotta (1951) was first to propose an equation of transport form like (8) to 

describe the evolution of q

2

l according to the two-point correlation equations derived from the 

Navier-Stokes equations (see also Kantha 2004).

The determination of the closure constants in (8), like that for the TKE equation, is 

empirical. We recommend that the following set of values be used for the results that are in the 

best agreement with those from the ETA scheme:

( , , , ) (2.75, 3, 3.0, 0.2) E1 E2 F Sl = .

As discussed later, they are determined based on the close analogy between the TKE equation 

and the length-scale equation and the properties of the turbulent flow in the neutral surface layer. 

This analogy is regarded as physically appropriate if the length scale is defined as the 

characteristic length scale of the largest energy-containing eddies and is related to the distance 

that these eddies travel in the vertical direction before losing their initial TKE due to turbulent

mixing and buoyancy effects. 
6

It should be pointed out that the values of the constants in both Mellor and Yamada 

(1982) and Kantha and Clayson (1994) were determined based on reliable measurements of 

decaying grid-generated homogeneous turbulence or turbulence near a wall without thermal 

stratification. The values that they recommended are 

( , , ) (1.8, 1.8, 0.2) E1 E2 Sl = and 































 = +

2

3 4 1

L

l F E E

κ ,

where ( E3 , E4 ) = (1.0, 1.33), L is the distance from the boundary and κ is the von Karman 

constant. It is assumed that the applicability of the constants obtained in such a way can be 

extrapolated to other flow situations of more complexity with thermal stratification such as in the 

ABL where reliable and comprehensive measurements are relatively difficult to make. However, 

the validity of this assumption has not been theoretically and experimentally well established, 

leading to the uncertainties in the determination of the closure constants. Recently, Umlauf and 

Burchard (2003) proposed an analytical approach in which the properties of the two equations 

for the TKE and the length scale are used as extra constraints to calibrate the constants. To 

include the buoyancy effects for stably-stratified flows, Freedman and Jacobson (2003) proposed 

how to modify the constants in the length scale equation that are calibrated with measurements 

of flows without thermal stratification. Despite the extensive effort to optimize the values of 

these closure constants for flow situations of different complexity, consensus has not been 

reached on the possibility of seeking a set of universal values suitable for the simulation of 

vertical turbulent mixing in geophysical flows, reflecting the fact that it is still an ongoing 

research how to properly specify the length scale as an essential part of parameterized turbulence 

modeling.

Boundary conditions are required to solve (1) and (8). In our implementation of the 

scheme in MM5-V3, the lower boundary condition for the TKE equation is provided using the 

surface flux algorithm used in the ETA scheme (Janjić 1990, 1994). This surface flux algorithm 

is based on the well accepted Monin-Obukhov similarity theory with the correction for the free 

convective limit proposed by Beljaars (1994). The lower boundary condition for q

2

l is 0.
7

3. The Closure Constants of the q2l Equation

The physical meaning of the length scale is of fundamental importance to the 

simulation of vertical turbulent mixing in the ABL. In the literature on the simulation of vertical 

turbulent mixing in the ABL, the length scale l is most often interpreted as the characteristic 

length scale of the largest energy-containing eddies (see, e.g., Moeng and Wyngaard 1989 and 

Wyngaard 1992). It is directly related to the distance that these eddies travel in the vertical 

direction before they loss their initial TKE due to turbulent mixing and buoyancy effects. In 

some references, the length scale is also linked to the so-called integral length scale (see, e.g., 

Kantha and Clayson 2000, p. 16). Strictly speaking, the integral length scale is associated with 

unbounded isotropic turbulence and is three-dimensional. The length scale l discussed here is 

associated turbulence in the atmospheric boundary layer, which is much more complicated than 

isotropic turbulence due the existence of thermal stratification in the vertical direction. The 

nature of the closure as described in (1)-(5) requires that the direction of l be along the vertical 

wind shear and thermal gradient. 

The closure constants in the length-scale equation must be determined empirically. 

Baumert and Peters (2000) indicate that the original values proposed by Mellor and Yamada 

(1982), and later adapted by Kantha and Clayson (1994), are not able to simulate the most basic 

properties of weakly stratified shear flows. They showed that for the well-understood near-wall 

flow situations without thermal stratification, there exist constraints on the determination of the 

values of the closure constants. Here we use a different approach to address the issues related to 

how to determine these constants.

The equation for the length scale, (8), can be rewritten as

( )



 





∂

∂

∂

∂ = − − +

∂

∂

+

∂

∂

q l

z

qlS B z

q E l K S E l K N F

t

l

q

t

q l m H l

2

1

3

2

2

2

1

2

2

 , (9)

where 

2 2

2

















∂

∂

+ 















∂

∂ = zV

z

U S , z

N g v

∂

∂Θ = β 2 . (10)
8

Substituting (1) into (9) and using the so-called Kolmogorov-Prandtl relationships (i.e., qlSq = 

K m σ m / and qlSl = K m σ l / where both σ m and σ l are empirically determined constants) yield

( ) ( ) ( )

( ) ( )



 





∂

∂

∂

∂ −



 





∂

∂

∂

∂

+

= − − − − − ∂

∂

2 2

1

3

2

2

2

1

2

/ /

2 2 2

q

z

K

z

q l l

z

K

z

B

q E l K S E l K N F

t

l

q

m l m m

m H

σ σ

 . (11)

Now let us now apply (11) into a steady neutral surface layer in which ∂l/∂t = 0, 2 N = 0, 

S u /κ z = * , 2

*

2 1/ 2

q 2 c u − = μ and l =κ z (where * u is the surface friction velocity, κ is the von 

Karman constant and μ c is an empirical constant). It is straightforward to verify that if σ l = 1, 

the following relationship holds for the neutral surface layer:

( ) 1

2 1/ 2

2

1 = c F − E − κ μ . (12)

This equation imposes an important constraint on the relationship between E1 and F . It also 

implies that only the difference between E1 and F controls the balance among the production, 

dissipation and diffusion of the length scale, l, and therefore q

2

l in the neutral surface layer. This 

implication is important for the determination of both E1 and F as shown later because no 

matter how F is determined, E1 must be determined through (12) for it to be dynamically 

consistent with F . A consequence of applying (12) to the determination of E1 and F is that if 

F involves a wall function such that it depends on z, E1 should also be dependent on z such 

that the difference between them is independent of z. On the other hand, F must be a constant 

according to (12) if E1 is chosen as a constant.

Baumert and Peters (2000) proposed E1 = 2 according to the hypothesis introduced by 

Tennekes (1989) that the integral length scale does not, based on dimensional grounds, depend 

on the homogenous shear. It should be noted that setting E1 = 2 implies that the evolution of 

the length scale is controlled by the buoyancy and dissipation terms. Given the fact that 
9

turbulence grows in a stratified unstable flow as well as in a weakly stratified stable flow, one 

has to use a positive E2 for the unstable flow and a negative E2 for the stable flow (see, e.g., 

Baumert and Peters 2000). Although this may not pose great difficulty in the simulation of 

vertical turbulent mixing in the oceanic mixed layer where the stratification is often weakly 

stable, it is definitely not desirable for the ABL because with a diurnal cycle, it is common for 

the stratification within the ABL to vary between stable stratification and unstable stratification. 

Moreover, although Tennekes’ hypothesis may be acceptable for idealized unbounded turbulence 

with homogenous shear, it is perhaps too simple to be applicable to turbulence in the ABL where 

the shear is hardly homogenous and the buoyancy effect on turbulence is much stronger than 

most of the engineering flows in laboratories (more discussion on this point is provided in 

Kantha et al. 2005).

We wish to argue that in the simulation of vertical turbulent mixing in the ABL, if the 

length scale is defined as the characteristic length scale of the largest energy-containing eddies 

and is regarded as being related to the distance that these eddies travel in the vertical direction 

before they loss their initial kinetic energy due to turbulent mixing and buoyancy effects, its 

evolution as described by (11) should be closely analogous to that of the TKE in that the 

deviation from the balance between the turbulence production (by both shear and buoyancy) and 

the turbulence dissipation controls the change of the length scale. During the initial development 

of turbulence, the characteristic length scale of the largest energetic turbulent eddies increases 

with the intensification of turbulence. The greater the characteristic length scale is, the farther 

the largest energy-containing eddies are able to move in the vertical direction before they 

completely mix with the surrounding environment. The intensity of turbulence and therefore l

cannot increase indefinitely because turbulent mixing tends to smooth out the vertical velocity 

shear and unstable stratification that are favorable for turbulence production. Therefore, the 

evolution of the characteristic length scale of the largest energy-containing eddies is governed by 

the same processes that control the evolution of the TKE. In order for the length-scale equation 

(11) to be valid for both unstable and weakly stable stratifications with the same set of closure 

constants, it is required that all the parameters before the turbulence production terms (i.e., the 

shear production and the buoyancy production) and the dissipation term be positive. With this 

requirement, (11) indicates that the lower bounds of the closure constants should be 2. 
10

Requiring F > 2 is not in agreement with the conventional value of F (which is close 

to 1.0 corresponding to the decay exponent of ~1.0; see, e.g., Durbin and Pettersson-Reif 2001 

and Umlauf and Burchard 2003) that is derived from the experimental data of decay exponent for 

homogenous, isotropic turbulence. If the behavior of the vertical turbulent mixing in the ABL is 

the same as in the turbulence in thermally unstratified flows, the values of all the closure 

constants can be determined uniquely from simple measurements in laboratories that isolate each 

term. However, there is no observational evidence to indicate whether or not this is indeed the 

case. 

The simulation of vertical turbulent mixing in weather prediction models perhaps should 

be different than the simulation of engineering turbulence because the former emphasizes the 

effect of the largest energy-containing (i.e., flux carrying) eddies on vertical turbulent mixing in 

a thermally stratified environment. Because the two-equation turbulence model is not exact and 

the thermal stratification strongly affects vertical turbulent mixing in geophysical flows, it may 

not be practically possible to make the closure constants universally applicable to geophysical 

and engineering flows. Take the night-time atmospheric boundary layer for example. The 

largest eddies travel in the vertical direction against the restoring force caused by the stable 

stratification. According to the scaling argument made by Wyngaard (1992), this implies that in 

an equilibrium, stably stratified ABL, the distance (i.e., l) that the largest energy-containing 

eddies travel is proportional to the square root of the eddies’ initial TKE. Thus, if the thermal 

stratification is held constant, l increases (or decreases) as q increases (or decreases). When 

F > 2 is chosen, the behavior of our length-scale equation appears to be in agreement with this 

argument, but apparently contradictory to the prediction based on the data of homogeneous, 

isotropic turbulence. However, if F is chosen to be smaller than 2 as done conventionally for 

engineering-turbulence modeling, the prediction of (11) is that the decrease of q due to the 

turbulence dissipation will lead to the increase of l, which is contrary to the above argument.

Currently, there is no practical way to observationally/experimentally determine the 

decay rate of the length scale as a function of the turbulence dissipation rate in the ABL. This 

study suggests that it is not unreasonable to set F = 3.0. Because typically μ c = 0.09, κ = 0.4 

(see, e.g., Durbin and Pettersson-Reif 2001) and E1 should be slightly smaller than F , this 

leads to E1 ≈ 2.9. In accordance with the consistent analogy between the length-scale equation 
11

and the TKE equation, E2 should be close to E1 , assuming that the Prandtl number is very 

close to unity. Results from our numerical experiments indicate that E2 = 3 is acceptable 

when E1 varies between (2.5, 3.0). Hence, the set of values recommended in section 2, 

i.e., ( , , , ) (2.75, 3, 3.0, 0.2) E1 E2 F Sl = , are obtained for the results that are in the best

agreement with those from the ETA scheme.

Summarizing the above discussions, the close analogy between the evolutions of the 

length-scale and the TKE should be used to determine the constraints on the values of the closure 

constants in the length-scale equation. This is based on the interpretation that l is the 

characteristic length scale of the largest energy-containing eddies, which is directly related to the 

turbulent mixing of these eddies in the vertical direction. It is required that the values of E1 , 

E2 and F be greater than 2. E1 should be slightly smaller than F , and E1 and E2 should 

be close to each other. The actual values of E1 and E2 are determined according to F . The 

ultimate validation of the value of F , and those of the closure constants, must rely on the 

comparison of the predictions of the scheme and the observations of the turbulence in the ABL.

4. Constraints on the Length Scale 

In general, the two-equation vertical turbulent mixing scheme does not ensure that the 

quantities of the TKE and l are both always within the physically meaningful range for a given 

flow situation. It does not even guarantee that both quantities are positive definite. Thus, 

constraints on the two quantities should be established to prevent nonphysical results. While it is 

straightforward to set the lower bounds for the two quantities (they should not be negative), 

setting the upper bounds is not straightforward.

It has been recognized for a long time that because the closure hypothesis and the 

specification of the length scale involve empirically determined closure constants, the 

MY closure has performance problems (see, e.g., Mellor and Yamada 1982, Galperin et al. 1988, 

and Burchard 2002). To alleviate the problems requires that the shear and stability parameters as 

well as the length scale be constrained properly. A proper specification of the limiting 

constraints on the shear and stability parameters as well as the length scale is very critical for the 

successful application of the MY closure. Various approaches to imposing the constraints on the 

length scale have been proposed for oceanic applications to constrain the shear and stability 
12

parameters along with the length scale. However, these approaches are not as theoretically 

appealing and straight forward as that proposed by Janjić (2002). Janjić’s approach is 

innovatively different than the others in that it only applies constraints on the length-scale by 

requiring that the solution of the TKE equation be nonsingular. The effectiveness of this 

approach was illustrated in the implementation of the ETA scheme. The essence of the approach 

can be summarized as the following.

When the thermal stratification is unstable corresponding to a negative gradient 

Richardson number (Ri), the TKE production and dissipation, i.e., the right-hand side of (1) 

should be nonnegative and bounded. Following Janjić (2002) and using the definitions in (4) and 

(7), the contribution of the TKE production and dissipation can be written as

) [ 1/ ]( ) 2 ( 3

/ 1

2

S G S G B q l

t

q

prod dissp = M M + H H − ∂

∂ . (13)

Substituting the definition of GH and GM along with SH and SM obtained by solving (6a) and (6b), 

(13) can be written in the following form





































−

















+ 















+ 















































+ 















 = − 

 





∂

∂

1

4 2

4 2

1

1

1

B

q

l

q

l

q

l

q

l

q t

l

γ δ

α β

 , (14)

where

[ ( ) ( ) ( ) ] H M H α A A β g g A A A B C A C B g β g g 1 2 2 2 1 1 1 2

2 2 2

1 2 = − 9 + 3 3 + 3 +18 −

[ ( ) ( ) ] M H β A C g A β g g 1 1 2 = 1− 3 −

[ ( )( ) ( ) ( ) ] A A A B β g g H A A B A g M β g g H γ 2 2 2

2

1

2 2

1 2

2 = 9 1 2 12 + 3 + 2 − 3

[ ( )( ) ] δ A gM A2 A1 B2 β g gH

2 = 3 2 1 + 7 +
13

z

g v h ∂

∂ = θ

2 2





 





∂

∂  + 

 





∂

∂ = z

V

z

U

gM

In order for the right-hand side of (13) or (14) to be bounded and nonnegative when the 

stratification is unstable or neutral, it requires that 

( )

1/ 2

q P1 l < (for Ri < 0) , (15)

where P1 is the nonnegative root of the following equation:

1 0

4 2

+ = 















+ 















q

l

q

l

γ δ .

It is known that an equilibrium solution may exist when the thermal stratification is 

neutral or stable (i.e., GH ≤ 0), but the gradient Richardson number (Ri) is smaller than a critical 

value (Ric) (see, e.g., Gerrity et al. 1994). This critical value of the gradient Richardson number 

can be obtained by letting the right-hand side of (13) or (14) be zero. It is equivalent to the 

following:

0 1

1

1

4 2

4 2

− =

















+ 















+ 















































+ 















B

q

l

q

l

q

l

q

l

γ δ

α β

, (16)

which can be rewritten as

0

4 2

 + = 

 



  +  





χ κ

l

q

l

q (17)
14

with 

[ ( ) ]( ) [ ( )] A A B A B g g H A A 1B 1 C 1 g M

2 χ = 3 2 7 1+ 2 + 2 1 β + 6 1 − 1− 3

[ ( )]( )

[ ( ) ( )] ( ) M H

H

A A B A B C A C B A A B A g g g

A A B A A A B g g

β

κ β

2 2 2

2

1 2 1 2 2 1 1 1 2 1

2 2

1 2

2

1 1 2

2

1 2

3 3 3 18 18 3

9 9 12 3

+ + − + −

= + + .

By definition, when Ri > Ric, q / l vanishes, and thus κ = 0 . This implies 

[ ( )]

( ) ( ) 2 2 2

2

1 2 2 2 1 1 1 2 1

1 2

2

1 1 2

2

1 2

3 3 18 6 3

3 12 3

A A A B C A C B A A B A

A A B A A A B

Ric

+ + − + −

+ + = . (18)

For the empirical constants listed in section 2, it is obtained that Ric = 0.505.

Since q and l may not vanish when GH ≤ 0 and the gradient Richardson number Ri < Ric, a 

constraint on l needs to be established accordingly. Again, following Janjić (2002), one can 

derive such a constraint by analyzing the diagnostic equation for 2 2

w q in the MY closure. 

Appealing to the expression for 2 2

w q in the Mellor and Yamada level 2.5 turbulence model 

(Mellor and Yamada 1982), Janjić showed that q/l depends on 2 2

w q through the following 

equation:

0

4 2

 + = 

 



  +  



 ξ σ ω

l

q

l

q , (19)

where

( ) ξ = 1− 3R s

σ {[18 A C 1 g m ( 9 A 1A 2 3 A 2 B 2 )( β g)g H ] 3R s δ } 2 = 1 + + −
15

ω {[27 A A B (β g) g H 54 A A 2 B 2C 1 g M (β g)g H ] 3R s γ} 2

1

2 2

2

2 = 1 2 + −

2 2 Rs = w q

It can be shown then that the smallest possible value of 2 2

w q (corresponding to q/l = 0) 

for the above given closure constants is 0.144. It can also be verified that the greater 2 2

w q is, 

the greater q/l. Therefore, the constraint on l is

( )

1/ 2

q P2 l < (for Ri ≥ 0) , (20)

where P2 is the nonnegative solution of (19) with 2 2 R w q s = = 0.144.

It is important to emphasize that the constraints derived above are based on the TKE 

equation and thus irrelevant to how the length scale is specified. In other words, no matter what 

prognostic or diagnostic approach to obtaining the length scale, it should be limited by the same 

constraints as expressed by (15) and (20), so long as (1)-(7) are used for the TKE prediction.

5. Numerical Integration 

The implicit forward in time and centered in space differencing scheme is commonly 

used to discretize the two-equation model, (1) and (8), for numerical integration, in which the 

production and dissipation terms are specified at the current time. This discretization scheme 

results in a set of linear algebraic equations that can be easily solved with the Thomas algorithm 

(see, e.g., Fletcher 1991). However, numerical experiments with this scheme often produce 

results that are contaminated with computational mode (see an example later in section 6). To 

circumvent this problem, the equations for the TKE and the length scale are solved using the 

numerical scheme developed by Janjić (2002), in which the TKE equation is solved first, and 

then the updated TKE is used in the integration of the length-scale equation. This numerical 

scheme solves either the TKE equation or the length-scale equation in two steps. Each equation 

is first integrated to obtain the intermediate solution with respect to the production and 
16

dissipation terms using an iterative scheme derived by linearizing the production and dissipation 

terms. Then, with the intermediate solution as input, the final solution is obtained using the 

Thomas algorithm to integrate the equation with respect to the diffusion. 

Because the split method is used to integrate the model, the length scale in the TKE 

equation can be moved into the differentiation sign on the left-hand side of (14), yielding





































−

















+ 















+ 















































+ 















 = − 

 





∂

∂

1

4 2

4 2

1

1

1

B

q

l

q

l

q

l

q

l

q t

l

γ δ

α β

 . (21)

By denoting the right hand-side as R, the solution of (21) for the next time step, i.e., the 

intermediate solution of the TKE equation with respect to the production and dissipation terms,

can be iteratively obtained using the following equation

( ) ( ) [( ) ] [( ) ]

{ [ ( ) i ] [( ) i ] [( ) ( ) i ]} { [( ) i ]} i i i i

R l q R l q l q l q t R l q

l q l q R l q R l q

'

0

'

'

1

+ + − exp ∆

+ = −

where i is the index of current time step, (/q)i is the value of /q around which the linearization 

is performed, (/q)0 is the initial value of /q at the beginning of the time step, and

( ) [( )( ) ( ) ( )] [ ( ) ( ) ] ' 5 3 4 2 2 R = d R d l q = − 2 α δ − β γ l q + 2α l q + β l q γ l q + δ l q +1 .

 

As proposed by Janjić, the iterations are started from the equilibrium solution of the 

equation (21) to ensure that the solution behaves well. After the solution converges, the 

diffusion operation is applied to obtain the final solution of the TKE equation.
17

For the length-scale equation, the counterpart of (21) is

( )





































− −

















+ 















+ 















































+ 















= − ∂

∂

1

4 2

4 2

2

1

( ' '

1/

B

F

q

l

q

l

q

l

q

l

t

l

q

γ δ

α β

 , (22)

according to (11), where

( )

3 ( 2)(3 3 18 ) ( ) ]

' [9 ( 2)

1 2 1 2 2 1 1 1 2

2 2

2

2

1 2

M H

H

A A E A B C A C B g g g

A A E g g

β

α β

+ − + + −

= − −

and

[ ( ) ( ) ] β A E C g M A E β g g H ' ( 2) 1 3 ( 2) = 1 1 − − 1 − 2 2 − . 

The procedure to solve this equation is the same as (21). As with in solving the TKE equation, 

after the converged solution is obtained for the intermediate l with respect to the production and 

dissipation terms, the corresponding diffusion operator is applied to the variable q

2

l, and thus the 

final solution of l is yielded.

6. Numerical Experiments

The mesoscale model in which the two-equation vertical turbulent mixing scheme is

implemented is the National Center for Atmospheric Research and the Pennsylvania State University 

mesoscale model (MM5-V3) (Grell et al. 1994 and http://www.mmm.ucar.edu/mm5/mm5v3.html). The 

mesoscale model is run on multiple 1-way nested meshes of 36 km, 12 km, and 4 km horizontal 

grid spacings (Fig. 1) with 50 vertical layers in total and 20 layers within the lowest 1 km. The 

initial and boundary conditions for the 36-km mesh are from the NCEP’s ETA Data Assimilation 

System (EDAS). All the simulations are started at 0000 UTC 30 July 2000. The configuration 

of MM5 physics for the 36-km and 12-km meshes includes the mixed-phase cloud physics, the 

Grell scheme for subgrid (convective) condensation (only for 36-km and 12-km meshes), the 

Dudhia’s 5-level simple soil model, and the MM5 simple short-wave and long-wave radiation 

parameterization schemes. The 4-km mesh is chosen to be small (with 19 x 19 horizontal grids) 
18

so that within the entire mesh the landuse characteristics are homogeneous. Since the synoptic 

forcing is weak in this case, we think that the extension of the domain is sufficient for the 

purpose of testing the two-equation vertical turbulent mixing scheme, and compare it with the 

ETA scheme that has been implemented in MM5-V3.

The differences in the length-scale distributions from the two-equation scheme and the 

ETA scheme are shown in Fig. 2 with an east-west cross-section cutting through the middle of 

the 4-km domain that is valid at 47 hours into the simulation (1700 local standard time). 

E1 = 2.875, E2 = 1.0 and F = 3.0 are used in the two–equation scheme. It is seen that the two 

schemes yield different magnitudes and gradient of the length scale for a given location. 

Associated with the differences in the length scale are the differences in q

2 (i.e., twice the TKE), 

the eddy exchange coefficients ( K m and K H ) and virtual potential temperature (Θv ) are 

depicted in Fig. 3. Significant differences between the two schemes are seen in the TKE and the 

eddy exchange coefficients fields. It can also be seen that the two-equation vertical turbulent 

mixing scheme produces more mixing of Θv near the surface than the ETA scheme, while far 

above the surface the former mixes slightly less that the latter. There is no significant difference 

in the ABL heights (as defined as the inversion of Θv at the top of ABL) simulated respectively 

by the two schemes, suggesting that observations of turbulence fields within ABL will be more 

discriminating than the ABL height for the validation of ABL models.

It is desirable to have some idea of how the performance of the two-equation scheme is 

sensitive to the variation of E1 and E2 . Therefore, Figs. 4-7 present a sample of the results from

the simulations using the two-equation scheme with different sets of values for E1 and E2 . If 

E2 is set to 1 and E1 varies, one can see that although the vertical distribution of the TKE with 

various E1 is quite similar, the maximum TKE increases only slightly as E1 decreases. On the 

other hand, if E1 is set to 1 and E2 varies, the vertical distribution of the TKE does not change 

very much with E2 , but the maximum of the TKE increases as E2 decreases, indicating that the 

TKE is more sensitive to variation in E2 than in E1 .
19

FIGURE 1. Location of the 3 nested MM5 domains.

 

a b

FIGURE 2. East-west cross-sections of the length scale cutting through the middle of the 4-km 

domain valid at 47 hours into the simulation (1700 local standard time). (a) is for the twoequation scheme with E1 = 2.875, E2 = 1.0 and F = 3.0, and (b) is for the ETA scheme. The 

contour interval is 10 m.
20

 a b

 

c d

 

e f

FIGURE 3. Continued on next page.
21

 g h

FIGURE 3. Cross-sections of 2

q (a and b, contour interval of 0.25 m2

s

-2

), m k (c and d, contour 

interval of 10 m2

s

-1

), H k (e and f, contour interval of 10 m2s

-1

), and Θv (g and h, , contour 

interval of 0.1 K). The left column (i.e., a, c, e and g) are for the two-equation scheme with E1 = 

2.875, E2 = 1.0 and F = 3.0, while the right column (b, d, f, and h) are for the ETA scheme. The 

cross-sections are at the same location as in Fig. 2 and are valid at 47 hours into the simulation 

(1700 local standard time).

 

It is mentioned in section 5 that a commonly used finite-difference scheme to discretize 

the two-equation model, (1) and (8), often leads to a solution that is apparently contaminated by 

computational modes. In this scheme, (1) and (8) are discretized using the implicit forward in 

time and centered in space differencing scheme except that the production and dissipation terms 

are specified at the current time. To illustrate the contamination of computational modes, the 

experiment E1 = 2.875, E2 = 1.0 and F = 3.0 is rerun with the implicit scheme and the results are 

shown in Fig. 8. It is seen that there are multiple extremes in the distributions of the TKE and 

the eddy exchange coefficients, but not in Θv . These multiple extremes cannot be explained 

physically given the smooth mean state, leading to the conclusion that they result from the 

contamination of computational modes. It is for the purpose of eliminating the computational 

contamination of the solution that Janjić developed his innovative numerical scheme to 

integrating the TKE equation (2003, personal communication).
22

 a b 

 c d

 e f 

FIGURE 4. Cross-sections of 2

q (contour interval of 10 m2

s

-2

) from the simulation of the 4-km 

domain using the two-equation scheme with (a) E1 = 1.0 and E2 = 0.5, (b) E1 = 1.0 

and E2 = 0.75, (c) E1 = 1.0 and E2 = 0.1, (d) E1 = 0.5 and E2 = 1.0, (e) E1 = 0.75 and E2 = 1.0, 

and (f) E1 = 0.875 and E2 = 1.0. The cross-sections are at the same location as in Fig. 2 and are 

valid at 47 hours into the simulation (1700 local standard time).
23

 a b

 c d

 e f

 

FIGURE 5. Cross-sections of m k (contour interval of 10 m2s

-1

) from the simulation of the 4-km domain using 

the two-equation scheme with (a) E1 = 1.0 and E2 = 0.5, (b) E1 = 1.0 and E2 = 0.75, (c) E1 = 1.0 and E2 = 

0.1, (d) E1 = 0.5 and E2 = 1.0, (e) E1 = 0.75 and E2 = 1.0, and (f) E1 = 0.875 and E2 = 1.0. The crosssections are at the same location as in Fig. 2 and are valid at 47 hours into the simulation (1700 local standard time).
24

 a b

 c d

 

 e f

FIGURE 6. Cross-sections of H k (contour interval of 10 m2s

-1

) from the simulation of the 4-km domain using 

the two-equation scheme with (a) E1 = 1.0 and E2 = 0.5, (b) E1 = 1.0 and E2 = 0.75, (c) E1 = 1.0 and E2 = 

0.1, (d) E1 = 0.5 and E2 = 1.0, (e) E1 = 0.75 and E2 = 1.0, and (f) E1 = 0.875 and E2 = 1.0. The crosssections are at the same location as in Fig. 2 and are valid at 47 hours into the simulation (1700 local standard time).
25

 a b

 c d

 e f

FIGURE 7. Cross-sections of Θv (contour interval of 0.1 K) from the simulation of the 4-km domain 

using the two-equation scheme with (a) E1 = 1.0 and E2 = 0.5, (b) E1 = 1.0 and E2 = 0.75, (c) E1 = 1.0 

and E2 = 0.1, (d) E1 = 0.5 and E2 = 1.0, (e) E1 = 0.75 and E2 = 1.0, and (f) E1 = 0.875 and E2 = 1.0. The 

cross-sections are at the same location as in Fig. 2 and are valid at 47 hours into the simulation (1700 

local standard time). 
26

 a b

 c d

 e

FIGURE 8. Cross-sections of (a) the length scale (contour interval of 10 m), (b) 2

q (contour interval of 

0.25 m

2

s

-2

), (c) m k (contour interval of 10 m2s

-1

), (d) H k (contour interval of 10 m

2

s

-1

), and (e) Θv (contour 

interval of 0.1 K) for the simulation using the two-equation scheme with the conventional implicit numerical 

scheme and E1 = 2.875, E2 = 1.0 and F = 3.0. The cross-sections are at the same location as in Fig. 2 and are 

valid at 47 hours into the simulation (1700 local standard time).
27

7. Summary and Discussion

A two-equation vertical turbulent mixing scheme is implemented in the National 

Center for Atmospheric Research and the Pennsylvania State University mesoscale model

(MM5-V3). This scheme is based on the MY closure with two prognostic equations: one for the 

TKE and the other for the length scale multiplied by twice the TKE. Although the scheme has 

been popularly used in oceanic modeling community, the work reported here is the first effort to 

implement the scheme in a mesoscale atmospheric model.

The physical meaning of the length scale is of fundamental importance in the 

simulation of vertical turbulent mixing in the ABL. In this work, the length scale is defined as 

the characteristic length scale of the largest energy-containing eddies, and is related to the 

distance that these eddies travel in the vertical direction before losing their initial kinetic energy 

due to turbulent mixing and buoyancy effects. With this definition, it is concluded that the 

closure constants in the length-scale equation should be different than those proposed previously 

for oceanic applications. 

In order to ensure physically sensible performance of the scheme in MM5-V3, the 

approach developed by Janjić (2002) is applied in the implementation of the two-equation 

scheme to impose the necessary constraints on the prognostic length scale. These constraints are 

derived by analyzing the TKE equation under different stability regimes in terms of the gradient

Richardson number. They are functions of the TKE as well as the stability functions in the TKE 

equation, and therefore are generally applicable to limiting the length scale predicted from 

various versions of the prognostic equations for the length scale. The use of these constraints 

circumvents the need for the realizability constraints on the non-dimensional vertical gradients of 

mean velocity and temperature. It is also shown that the numerical scheme proposed by Janjić 

(2002) for integrating the TKE equation should be used to solve the prognostic equations of the 

two-equation scheme to control the computational modes encountered when using conventional 

numerical schemes to solve the equations for the TKE and the length scale. 

The numerical experiments are performed using MM5-V3 to illustrate that with properly 

chosen closure constants in the length-scale equation, the two-equation scheme produces more 

(and less) mixing near the surface (and way above the surface layer within the ABL) during 

daytime than the original TKE based, ETA scheme in MM5-V3 where the length scale is 

prescribed using a diagnostic approach. The turbulence fields from the two-equation scheme are 
28

significantly different than those from the ETA scheme, while the difference in the ABL height 

is almost the same. A series of sensitivity experiments indicate that the two equation scheme is 

more sensitive to the closure constant associated with the wind shear than to that associated with 

the buoyancy effect in the length-scale equation. Although the use of the prognostic equation for 

the length scale is more theoretically appealing than the use of the diagnostic equation, the 

results presented here are not enough to generalize whether or not the gain from using the twoequation scheme is proportional to the extra computational cost for the simulation of the ABL 

turbulent mixing in weather prediction models.

It has been known in the ABL modeling community that a successful one-dimensional 

test and calibration of a vertical turbulent mixing scheme does not guarantee that the same 

scheme will be computational robust and physical meaningful when it is used in much more 

complicated three-dimensional models. We would like to emphasize that a well-tested oceanic 

vertical turbulent mixing scheme may not automatically work well when used in simulating the 

vertical turbulent mixing in the atmosphere because the thermal stratification of the atmosphere 

varies more quickly and widely than in the ocean, and winds in the atmosphere vary more 

rigorously than currents in the ocean. While we do not have adequate turbulence observations to 

indicate which one of the available vertical turbulent mixing schemes in thee-dimensional 

mesoscale model such as MM5-V3 is the best for a given flow situation, this two-equation 

scheme provides another choice in ABL parameterizations, which could contribute to the effort 

of effectively perturbing model physics for generating statistical ensemble forecasts. 
29

8. References

Baumert, H., and H. Peters, 2000: Second-moment closures and length scales for weakly 

stratified turbulent shear flows. J. Geophys. Res., 105, 6453-6468.

Beljaars, A.C.M., 1994: The parameterization of surface fluxes in large-scale models under free 

convection. Quart. J. Roy. Meteor. Soc., 121, 255-270.

Blackadar, A.K., 1962: The vertical distribution of wind and turbulent exchange in neutral 

atmosphere, J. Geophys. Res., 67, 3095-3102.

Burchard, H., 2001: On the q2l equation by Mellor and Yamada (1982). J. Phys. Oceonogr., 31, 

1377-1387.

Burchard, H., 2002: Applied Turbulence modeling in marine waters, Springer-Verlag, Berlin, 

215 pp. 

Durbin, P.A., and B.A. Pettersson-Reif, 2001: Statistical Theory and Modeling for Turbulent 

Flow, John Wiley and Sons, Chichester, U.K., 282 pp.

Fletcher, C.A.J., 1991: Computational Techniques for Fluid Dynamics, Vol.1: Fundamental and 

General Techniques, Springer-Verlag, Berlin/New York.

Freedman, F.R., and M.Z. Jacobson, 2003: Modification of the standard ε-equation for the stable 

ABL through enforced consistency with Monin-Obukhov similarity theory. BoundaryLayer Meteorol., 106, 383-410.

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.

Gerrity, J.P., T.L. Black, R.E. Treadon, 1994: The Numerical Solution of the Mellor-Yamada 

Level 2.5 turbulent kinetic energy equation in the Eta model. Mon. Wea. Rev., 122, 1640–

1646.

Grell, G.A., J. Dudhia, and D.R. Stauffer, 1994: A description of the fifth-generation Penn State/

NCAR mesoscale model (MM5). NCAR Technical Note, NCAR/TN-398+STR, 122 pp.

Helfand, H.M., and J.C. Labraga, 1988: Design of a nonsingular level 2.5 second-order closure 

model for the prediction of atmospheric turbulence. J. Atmos. Sci,. 45, 113-132.

Holtslag, A.A.M., 2003: Boundary layers/Modeling and Parameterization, Encyclopedia of 

Atmospheric Sciences, Edited by J.R. Holton, J.A. Curry, and J.A. Pyle, Vol.1, 253-261, 

Academic Press, New York.
30

Janjić, Z.I., 1990: The step-mountain coordinate: physical package. Mon. Wea. Rev., 118, 1429-

1443.

Janjić, Z.I., 1994: The step-mountain eta coordinate model: further developments of the 

convection, viscous sublayer and turbulence closure schemes. Mon. Wea. Rev., 122, 927-

945.

Janjić, Z.I., 2002: Nonsingular Implementation of the Mellor–Yamada Level 2.5 Scheme in the

NCEP Meso model, NCEP Office Note, No. 437, 61 pp.

Kantha, L.H., 2004: The length scale equation in turbulence models, Nonlin. Processes

Geophys., 11, 83-97. 

Kantha, L.H., and C.A. Clayson, 1994: An improved mixed layer model for geophysical 

applications, J. Geophys. Res., 99, 25, 235-25,266. 

Kantha, L.H., and C.A. Clayson, 2000: Small Scale Process in Geophysical Fluid Flows. 

Academic Press, 888 pp.

Kantha, L.H., J.-W. Bao, and S. Carniel, 2005: A note on Tennekes hypothesis and its impact on 

second moment closure models. Ocean Modelling, 9, 23-29.

Mellor and Yamada, 1974: A hierarchy of turbulence closure models for planetary. J. Phys., 

Oceanogr., 5, 718-128.

Mellor, G.L., and T. Yamada, 1982: Development of a turbulent closure model for geophysical 

fluid problems, Rev. Geophys. Space Phys., 20, 851-875. 

Moeng, C.-H., and J.C. Wyngaard, 1989: Evaluation of turbulent transport and dissipation 

closures in second-order modeling. Jour. Atmos. Sci., 46, 2311–2330.

Rotta, J.C., 1951: Statistische Theorie nichthomenger Turbulenze, 1, Z. Phys., 129, 547-572.

Tennekes, H., 1989: The decay of turbulence in plane homogeneous shear flow. In Lecture 

Notes on Turbulence, eds. J.R. Herring and J.C. McWilliams, World Scientific, 32-35

Umlauf, L., and H. Burchard, 2003: A generic length-scale equation for geophysical turbulence. 

J. Mar. Res., 61, 235-265. 

Wyngaard, J.C., 1992: Atmospheric turbulence, Annu. Rev. Fluid Mech., 24, 205-233.