The equations and their solution for the surface layer of the boundary layer have been changed in several ways. The most extensive change has been the elimination of the iterative solutions for the Obukhov length, L. Additional changes include the recognition of a "rough sub-layer", Z*, and also to changes in the formulations that apply to strongly stable and unstable bulk Richardson number, RB. For the stable case, the restriction that the bulk Richardson not exceed a "critical" Richardson number has been removed. For the unstable case, a restriction has been placed upon \Zo/L\.. That is, \L\ is not permitted to approach Zo. Ignoring this restriction can result in unrealistic upward sensible heat fluxes that exceed by several times the solar content. We shall consider the stable case first. [N. B.: the notation used in this appendix conforms to the notation used in the main body of the text].

(a) Monin-Obukhov functions FM,H (z = Z/L) that are linear for either all z > 0 (or for all sufficiently large z) yield critical Richardson numbers. That is, R --> RC for z --> + . As the critical number is approached, FM,H,Q vanish. For example if we take FM,H to be

where M, H and bM, H are constants, then the critical gradient Richardson number is given by

The critical Richardson number for the profile relations in the MRF is 1.32. This is about 6-7 times higher than the RC for the Dyer (1974) or for the commonly used Businger et al. (1971) profile relations. It is likely, however, that even RC = 1.32 is too low. When R equals or exceeds, RC, there is no downward heat flux. The main consequence of a low RC is that too much of the upward daytime heat flux is trapped in the boundary layer and cannot return during stable nocturnal conditions. A relaxation of the RC is needed to enhance the downward flux.

The lack of a precise knowledge of the stable boundary layer's profile laws permits some constructive tampering, even to the extent of eliminating (moving to +) the critical Richardson number. This should not be interpreted as a denial of the existence of RC (probably ~0.20) at a single point, but rather the recognition that an area average of the nocturnal boundary layer includes portions of small, patchy turbulent eddies mixed with gravity waves. Flux leakage is realistic and reasonable.

(b) In this section we examine a simple formulation for FM, H that exhibits neither a critical gradient nor a critical bulk Richardson number. In addition, the formulation does not require dividing z into separate weakly, mildly, and strongly stable regimes.

The formulation is motivated by the turbulent diffusion equations for R o that the MRF uses above the surface layer,

To make this relation compatible with the surface layer, we note the following surface layer relations,

As usual, FH = FQ. For mildly stable cases, experimental evidence suggests FM = FH. For increasing stability, it is possible that FM and FH begin to differ, but we will not deal with this possibility at this time. We shall assume that FM = FH for all positive z.

Eqs. (4.4) - (4.6) can also be written in 'K-theory' form:


follow by comparison of (4.121) - (4.123) with (4.124) - (4.126). These diffusion coefficients can be related to (3) provided that in the surface layer we have: the (neutral) mixing length: 1 = kz; the wind shear: ; the stability function:

since the connecting relation between z and R is

the 'universal' functions FM, H are simply

this means that



For small values of z (weak stability), we have the standard linear results

but for large , z we have

Eqs. 4.21 - 4.22 mean that for strongly stable conditions the fluxes of momentum and heat are nonzero and are given by

As is increased, the fluxes decrease.

(c) In this section we show how z can be estimated from the bulk Richardson number, RB, without iteration. We first use the equation that relates z to RB,

in which

From (4.132), we have,

since FM, H is conventionally written in the form

we must change the form of (4.142). The result is

in which

Another standard form for (4.143) is

in which

Now, some manipulation shows that FM,H can be approximated for= 5 by,

In addition, we can roughly approximate z by

which leads to the quadratic equation

in which

The physical correct solution to (4.152) is

The solution for z can be divided into two portions. Numerical computations show that z can be adequately approximated by

Substitution of the approximate value of into (4.144) leads to momentum and heat transfer, coefficients CM, H that differ only slightly from the 'exact' values. For the approximation should not be used and should be replaced instead by where is given by (4.153) - (4.154). As with should be substituted (4.144) in order to approximate FM, H which then leads to the transfer coefficients,

  1. The Unstable Regime
(a) As noted in the main body of the text, most of the relations for the unstable case (z<0) that have been used in numerical forecast or simulation models are of the form

where aM, H, M, H and pM, H are positive constants. Some examples are: (A) Businger et al. (1971), {aM, M and pM} = {1, 15, 1/4}; {aH, H and pH} = {0.74, 9, 1/2}; (B) Dyer (1967) aM=aH=1, M,H=15, pM=0.275, pH=0.55; (C) Carl et al. (1973), {aM, M and pM} = {1, 16, 1/3}; {aH, H and pH} = {.74, 16, 1/4}. The MRF physics formulation uses the results of Dyer (1974) and Hicks (1976):

(b) In this section we show how (4.157) can be used to compute FM, H = ln (Z/Zo) - YM, H (z,zo) = ln(Z/Zo) - YM, H + YM, H ( zo). When - zo is small, then (4.160) can be approximated by

as it commonly done. Now, for 0 - z 0.5, YM, H (z ) can be approximated by

To determine a0, a1, b1, ( 4.161) is expanded:

Also, we have,

By equating a0=-4, a1-a0 b1 =-20 and requiring YM, (-z=0.5) = 0.7934, we get

for -0.5. Additional accuracy can be achieved by forcing collocation at -z=0, 0.05, 0.10, 0.25, and 0.50:

Similar calculations for FH give,


We now consider the more difficult general problem of determining YM, H from YM, H = (1-z ) -p. The task is tedious, but the asymptotic result, Eq. 4.76 is fairly simple. As -z increases, the accuracy of YM, H increases.

We first express Y as

in which x = - Z/L. We also express (4.168) as Y(x) = I1 (x0)-I2(x)

in which


Integral I1 can be broken up into Ja and Jb:

The second term in J a is the exponential integral defined by (Arfken, 1985)

and given by the series expansion

for which is the Eulaer-Masheroni constant,

Substitution of (4.174) - (4.175) into (4.172) yields

We are left with a logarithmic divergence in Ja, but this will be shown to cancel with another term.

We now examine Jb by noting that an integral representation of the logarithm of the gamma function G(t) can be written as

The integral can be differentiated to get

which, of course, is the same as [1/G(u)] [dG(u)/du], the definition of the Euler 'digamma' or 'psi' function (Gradshteyn and Ryzhik, 1965),

We can change variable in the second term of Jb:

to get

and, therefore,

Values of the digamma function have been tabulated. Some values that are of interest to us are:

It remains to evaluate I2. We write

in which S (s;p) is the series

For FM (Z/L) with p= 1/4 and FH (Z/L) with p=1/2, we get


The presence of ln r in cancels the divergence in Ja. The final result is obligingly compact:

When (4.190) is applied to FM and FH, we get

We need to retain only the (-z) -1/4 and (-z) -1/2 terms for -z>0.5 to achieve an adequate accuracy of 1.8% or greater. As expected, as -z increases, the accuracy increases.

(c) In this section, we delve into the physical and computational difficulties that arise when h is not many times larger than Z0 and when |L|->Z0. These conditions are not normally encountered in field conditions whose goal is to measure FM, H (Z/L), but large Z0 and small |L| must be expected to occur fairly frequently in global modeling and must be dealt with.

We deal with the problem of Zo~O(h) first. It is not expected that the logarithmic wind law is valid for Z ~Zo. Tennekes (1973) has suggested that the region for validity is Z> 10-100 Zo. Following Garratt (1980), we denote by Z* the lowest Z for which the profile laws are approximately valid for neutral and unstable lapse rates. The region Z<Z* is Garratt's 'roughness sublayer'.

At heights lower than Z* the individual roughness elements can be 'sensed' as the result of turbulent wakes created by flow around individual elements. This violates the similarity conditions by introducing an additional length scale, Z*, into the conventional profiles laws, FM, H. Garratt's analysis of data from very rough terrain suggests that Z*35 Zo for momentum and that Z*100 Zo for heat. This means that if the largest roughness lengths in the MRF are say, Zo5m, then the standard profile laws do not apply to heights below 175m for momentum and 500m for heat. These depths are as great as or greater than the heights of the surface layers of most NWP models.

In our calculations, we shall take a more liberal approach and merely require h>10 Zo. For points where Zo>h/10 for the FM, H calculations.

We now deal with the second question. The Businger-Hicks profile laws are regarded as reasonably accurate for -z<3. Model computations, however, can significantly exceed this moderate limitation.

Numerical calculations have uncovered another difficulty. Consider the following situation: and h are held constant as U(h) increases -RB and -z. This decreases FH and FM, which, in turn, increases CH. Under normal conditions (that is, |L| >>Zo), the decrease in U(h) overcomes the increase in CH, and the surface heat flux decreases. For very large Zo and small |L|, there is an anomaly: As U(h) decreases, the heat flux decreases, reaches a plausible minimum, and then increases. Eventually, heat flux that equal and exceed by several times the solar constant can be (mathematically) reached. This is counter-intuitive.

These 'runaway' heat fluxes can be eliminated by noting the following pattern that begins to emerge after observing the behavior of FH for many examples: Although the examples do not seem to indicate a minimum -u* Q* for any particular value of U(h), RB, h, or a minimum occurs when

for m20-50. To help eliminate the problem of runaway heat and humidity fluxes, we choose

This choice has the added advantage of allowing us to neglect Y(Zo/L) compared to Y(Z/L). This means that YM,H (Z/L; Zo/L) can be accurately approximated by YM, H.

To achieve the restriction imposed by (4.194), U(h) is increased in

That is, we require

for the CM,H,Q and flux calculations. Once these quantities are computed, U(h) is returned to its original value.


Arfken, G., 1985: Mathematical Methods for Physicists, Third Edition. Academic Press, Inc., Orlando. 985pp.

Carl, D. M., T. C. Tarbell, and H. A. Panofsky, 1973: Profiles of wind and temperatures from towers over homogeneous terrain. J. Atmos. Sci., 30, 788-794.

Dyer, A. J., 1967: The turbulent transport of heat and water vapor in an unstable atmosphere. Quart. J. Roy. Meteor. Soc., 87, 401-412.

Garratt, J. R., 1980: Surface influence upon vertical profiles in the atmospheric near-surface layer. Quart. J. Roy. Meteor. Soc., 106, 803-819.

Gradshteyn, I. S. and I. M. Ryzhik, 1965: Table of Integrals, Series and Products, Fourth Edition. Academic Press, Inc., New York. 1086pp.

Hicks, B. B., 1976: Wind profile relationships from the 'Wangara' experiment. Quart. J. Roy. Meteor. Soc., 102, 535-551.

Tennekes, H., 1973: A model for the dynamics of the inversion above a convective boundary layer. J. Atmos. Sci., 30, 558-567.