Paul E. Long
January 1988

  1. Purpose of Subroutine PROGTN
PROGTN's main purpose is to predict the surface (air-ground interface) temperature and humidity and to estimate the fluxes of momentum, heat, and humidity in the surface layer of the atmospheric boundary layer; supply certain turbulent quantities required elsewhere in the MRF. To do this, PROGTN requires external calculations of long and short wave radiation; also, surface and subsurface quantities such as roughness length, snow cover, and soil thermal properties.

The subsections that follow explain the functioning, and the reasoning that underlies the functioning, of the major components of PROGTN.

We are unable, at the present time, to provide a complete documentation of the surface and subsurface processes in the model and will only provide the first part. The remaining sections are currently under preparation and w ill be furnished soon.

  1. Introduction to the Surface Layer
Motion within the atmospheric boundary layer is usually turbulent. Although the atmospheric boundary layer (ABL) is typically 1-2 km thick during vigorous daytime turbulent convection, it may shrink to only several decimeters at night during stable conditions. Under extremely stable conditions, turbulence can become so irregular that determining a unique height of the ABL may not be possible.

Subroutine PROGTN computes the turbulent fluxes of momentum, sensible heat, and latent heat () for both unstable and stable conditions within the MRF's surface layer. The surface layer of the atmosphere (less frequently used near-synonyms: constant flux layer; contact layer) occupies very approximately the lowest tenth of the boundary layer. That is, hs0.1 h, in which hs is the height of the surface layer (SL), and h is the height of the ABL.

The flux calculation serves two purposes. First, the sensible and latent heat fluxes are required for the equations that predict the air-ground interface temperature (T*). Second, can be used to specify the lower boundary conditions for the diffusion equations in VERDIF that simulate turbulent transport of momentum, heat, and specific humidity above the surface layer. At the present, however, the fluxes of PROGTN are not used directly in VERDIF's boundary conditions, but rather turbulent quantities from PROGTN are passed to MONIN which are then used to reconstruct the surface layer's turbulent fluxes.

The top of the ABL in the current MRF is not confined to a predetermined number of layers, in contrast with earlier versions in which all turbulent diffusion (and therefore the ABL) was restricted to 1.5 km above the ground. The surface boundary layer, however, is assumed to be the lowest model layer for the current and previous versions. Until recently, the MRF used a 56 mb thick surface layer. The current version's surface layer is, by contrast, only 10 mb deep. The earlier surface layer was much too deep. In the real atmosphere, 56 mb can contain much (for convective conditions) or all (for stable conditions) of the entire ABL. Flux calculations from the new MRF are, therefore, expected to be more realistic than those from the MRF that used the deeper SL.

A deep model surface layer creates several difficulties. First, the relations that are used to compute the surface fluxes may well be forced beyond the limits for which they are valid. The calculation of surface layer fluxes in the MRF depend upon the non-dimensional variable hs/L, where L, the Obukhov (or Monin-Obukhov) length, is a measure of surface layer stability. The Obukhov length is positive for stable cases, negative for unstable cases, and diverges to 1 as neutral conditions are approached. A small value of -L, coupled with a large hs, yields a large -hs/L, indicating a condition of extreme instability.

There are certain 'universal' functions used in most SL relations that depend upon the non-dimensional height Z/L. The experimentally verified range of validity of these SL relations is about -Z/L2.5 for unstable cases and Z/L2 for stable cases (e.g., Businger et al., 1971). Within this fairly narrow range of stability, the SL relations agree with field data reasonably well. Outside this range, it is uncertain whether the standard SL relations are valid. the few experimental field data that exist suggest that the standard SL relations are not valid (e.g., Carl et al., 1973). For strongly stable cases (Z/L1), the question of validity is particularly vexing, since there is no general agreement upon even the functional forms of the SL relations or if there is a 'critical' Richardson number (Ric) for which the SL fluxes vanish. Stated somewhat differently, it is not known if there is a Ri for which Z/L as RiRic.

For strongly unstable cases, on the other hand, it is not clear if the conventional SL relations used in PROGTN should be retained or if they should be replaced by some form of free convection relations. Existing empirical data are insufficient to definitively resolve the issue: lacking such a verdict, computational considerations become a prime factor. Free convection SL flux-profile relations have a computational advantage. With free convection relations, SL fluxes remain finite as the windshear is reduced to zero, in contrast to the standard relations for which heat and moisture fluxes become unbounded. Past experiments with the MRF have revealed isolated instances of inordinately strong heat fluxes. Unrealistically large surface fluxes in PROGTN can usually be suppressed by prohibiting the SL windspeed from falling below a fixed minimum value, which in PROGTN, is 1 ms-1. This tactic has eliminated most, but not all, unrealistic values of hs/L and ; however, several related problems remain.

First, strongly unstable potential (or virtual potential) temperature gradients can occasionally produce unrealistically large fluxes even if the SL windspeed remains above the minimum value. In addition, a compatibility problem, common to two-tiered ABL models, occurs at the interface of the SL and the second the layer of the MRF. The result: the diffusion coefficients are not continuous across the interface. The discontinuity intensifies with increasing -hs/L. To remove the discontinuity either the flux-profile relations in the SL or the diffusion coefficient relations above the surface layer must be reformulated. A simpler alternative is to substantially decrease hs. The considerable reduction of the height of the surface layer in the model is not a complete palliative. A very shallow surface layer, while reducing the chances of unrealistic simulated surface fluxes, creates fresh computational and parameterization problems.

One such parameterization conundrum occurs over regions with large roughness. Within these regions, surface flux calculations can be rendered invalid unless the (representative) Zo is very much smaller than the height of the surface layer. The largest roughness lengths in the current MRF are between one and ten meters. It is difficult to give a physically sound, compelling argument for these large values, but a plausible rationale is given in Baumgartner and Mayer (1977). The Zo field over land in the MRF is based upon their estimated values. We do not use their Zo fields over oceans. PROGTN computes its own Zo fields that depend upon the frictional stresses at the oceans' surface.

A strong caveat again the use of large Zo was given by Garratt (1977). His analysis of field data suggests that the standard flux-profile laws are invalid within the sublayers Z < 35Zo for momentum and Z100Zo for heat and (presumably) moisture. This disconcerting result implies that for Zo2m, the standard SL flux-profiles laws are highly suspect within the first 70 m above the surface for wind and 200 m for temperature and humidity.

Since hs is 10 mb in the MRF, Garratt's sublayer exclusion principle - if correct - indicates that the fluxes computed by PROGTN will be probably quite erroneous over the domains of maximum roughness in the MRF. We note that these domains cover a greater area in the ECMWF model (Tiedtke et al., 1981) than with the MRF. Arbitrarily reducing the values of Zo in the regions where Zo > 1 m is a crude but useful 'correction', as MRF experiments seem to show. The issue of the physical reality of large roughness lengths remains unsettled, as is the possible distinction between the Zo for momentum, heat, and moisture.

As noted, there can be computational problems created by the use of a very shallow surface layer. A thin surface layer implies that the adjacent layers are also relatively thin (layers 1, 2, and 3 in the MRF are 10, 17, and 25 mb deep). The adjacent layers in the MRF use turbulent quantities computed in the SL as lower boundary conditions for the diffusion equations ("K-theory" equations) for heat, momentum, and moisture. The diffusion equation for each of these variables is approximated by a fully (in a linear sense) implicit finite difference scheme. This implicit difference scheme for constant K and Z is known to be absolutely stable for all time steps (t). Computational experience, however, shows that, for layers of variable thickness and for diffusion coefficients that depend upon the local Richardson number, the accuracy, and even the computational stability schemes, can be severely degraded as the Fourier number ( ) is increased. Feed-back oscillations between the K profiles and the predicted , u, v, q profiles will cause the solutions to behave increasingly pathologically. While not necessarily unstable in the usual sense( the relations can become completely useless rather quickly. Our recent experiments show that much of the erratic behavior can be eliminated with an adaptation of a predictor-corrector method.

We see from the caveats discussed above that the choice of the depth of a model's surface layer is not arbitrary. The surface layer can be neither too deep not too shallow if unrealistic numerical numerical for flux parameterization computations are to be avoided. It is believed that the current value for hs used in PROGTN strikes a judicious balance.

  1. Surface Layer Relations: The Monin-Obukhov Flux Profile Laws
Two of the most frequently encountered problems in boundary layer meteorology are the determination, from minimal number of parameters, of the surface layer fluxes of momentum, heat, and moisture as well as the synthesis of the complete mean SL profiles of wind, temperature, and specific humidity. PROGTN computes the SL fluxes and profiles and employs results in the solution of the surface temperature over soil, snow, and sea ice. Turbulent quantities computed in PROGTN are also utilized as the lower boundary conditions for the computation of changes in wind, temperature, and humidity for the portion of the boundary layer that lies above the surface layer.

The focal point of most modern flux-profile relations is the similarity theory positioned by Obukhov (1946), and Monin and Obukhov (1953, 1954). The Obukhov similarity theory can be invoked to show that the turbulent fluxes in the SL are uniquely determined from six quantities: the local roughness length; the air-ground interface (surface) temperature and specific humidity ( Qs, qs); also, Q(Z), U(Z) (windspeed); and q(Z). The height Z can be any height within the SL, provided Z>>Zo. PROGTN uses this uniqueness property of Obukhov similarity theory to compute all the required SL turbulence variables required by the MRF.

We now delineate the basic Obukhov relations used in PROGTN. First, Obukhov similarity postulates the existence of 'universal' functions such that the nondimensional vertical gradients of wind, temperature, an d moisture are given by:

in which k is the von Karman constant ( k0.4), and b, c are nondimensional constants near unity; the constants b and c are the neutral turbulent Prandtl and Schmidt numbers. The turbulent scaling parameters ,, and have dimensions of wind speed, temperature, and humidity. These parameters are a measure of turbulent fluctuations and vanish when the SL becomes laminar.

Unfortunately, the original Obukhov similarity theory only assumes the existence of the FM,H,Q as functions of Z/L. No current theory, moreover, exists that accurately specifies the universal functions over the entire stable and unstable range found with the field data or the even more extreme range encountered in global forecast models. The so-called 'higher-order' closure models have been reasonably successful, however, in simulating measured flux-profile relations (see e.g., Mellor and Yamada, 1982).

The paucity of high-quality field data to definitively specify interpolation formulae for FM,H,Q over extreme stability ranges has created a proliferation of semi-empirical functions. Intuitively plausible but conflicting arguments are frequently marshaled to support one set of functions over another. Surveys by Dyer (1974) and Yaglom (1977) compare many such relations. The relations for the unstable and stable (Z/L > O) cases that have beenmost used in numerical models are given by Businger, et al. (1971),

in which k=0.35, b=0.74. Dyer (1974) and Hicks (1976) proffered comparable profiles:

The Businger relations are more often employed in numerical models than the Dyer-Hicks relations. Computational experiments, however, suggest that the Dyer-Hicks system provides a slightly more rapidly converging solution f or Z/L < O than the Businger relations. The MRF uses a rather contrived extension of the Dyer-Hicks relations in PROGTN and MONIN in order to circumvent an otherwise prohibitively small critical Richardson number. The calculation of L is central to subsequent calculations of surface stress, sensible and latent heat fluxes, and the eddy diffusion coefficients.

  1. Relation Between the Gradient Richardson Number and the Obukhov Length
There is a simple relationship between the gradient Richardson number and the scaled height Z/L, as can be seen by solving for , and from the Obukhov profile laws (4.1, 2) and substituting into one of the defining relations for L; that is,

In (4.8), is the standard buoyancy parameter, and Qo is the so-called kinematic heat flux (). It should be pointed out that these equations are 'dry' and should be modified when the latent heat flux assumes a significant role in the surface energy balance. The result of substituting (4.1,2) into (4.8a) is

in which Ri is the local gradient Richardson number, and we have assigned Z = hs in (4.9).

We see that L, , and the 'dry' gradient Richardson number must have the same sign; that is, Ri, L, and are negative when the sensible heat flux is positive (upward), and are positive when the sensible heat flux is negative (downward). The sign of should always be positive; a negative is always unphysical. A negative is a spurious solution of (4.9, 4.6, 4.7b), for example, in the strongly stable case when the critical Richardson number is exceeded. Negative are avoided in general and in PROGTN by constraining Ri to be 0 < Ri < Ric.

The formal (mathematical) question of whether a critical Richardson number exists hinges upon the functional forms of the stable flux-profile laws. That a critical Richardson number exists for the Dyer-Hicks relations follows immediately from the linear form of (4.6b) and (4.7b). Substitution into (4.9) yields

Taking the limit , we see that the Richardson number steadily increases to the limiting (critical) value,

For the general linear profile laws given by


the asymptotic value - the critical Richardson number - is

For the Businger profile laws the critical Richardson number is 0.21. The gradient Richardson number cannot exceed Ric without causing Z/L to become negative, a mathematically correct but physically prohibited solution. If , as is common in a model, Ri exceeds Ric, then either Ri ormust be decreased or must be increased to force Ri < Ric.

Such relatively small critical Richardson numbers which follow from the Hicks and the Businger profiler laws are too restrictive for use in PROGTN. It is a feature of PROGTN that the linear flux-profile laws (4.6,4.7b) are not considered valid for all positive Z/L but rather only over the mildly stable range 0 < Z/L < 0.5. For the stronger stability range 0.5 < Z/L < 10, and for the extremely stable range (in which turbulence is largely extinguished), PROGTN uses,


There is little empirical support for these relations that were first proposed by Carson and Richards (1978) as an extension to the linear Hicks flux-profile laws. Their major effects in PROGTN are to: shift Ric from 0.20 to 1.32; allow a small downward surface heat flux within the otherwise super-critical range 0.2 < Ri < 1.32; needlessly complicate PROGTN by imposing iterations with little evident practical gain. To simplify the stable flux-profile system used in PROGTN while maintaining the small downward super-critical heat flux requires new flux-profile relations in PROGTN. Although not particularly difficult, this is not an approach we will pursue further here. A change in the flux-profile relations will be made within PROGTN only after the influence of the current formulation upon the heat and vapor flux and also the diurnal temperature wave is better understood. Another reason for temporizing is the need to construct a suitable stable SL formulation that is also comparable with the existing K-relations above the SL. Suitability and compatibility are not trivial constraints.

  1. The Relation Between the Bulk Richardson Number and the Obukhov Length
The stable and unstable relations for Z/L given above depend upon local gradients of temperature, windspeed, and humidity. It is important to note that none of these quantities is known in PROGTN, however, where the SL is treated as a bulk layer with no internal refinement. Only bulk differences are given, for example, , between the ground and mid-way into the SL (or at the top of the SL)]; similarly, for q and U. Accordingly, the Richardson numbers used in PROGTN is not the gradient Richardson number as used above, but rather a bulk Richardson number (RiB). This bulk Richardson number is defined ('dry form') by

where it is assumed that the wind speed vanishes at Zo. Calculations with RiB and hs/L require that the flux-profile relations be integrated from Zo to hs where hs >> Zo (a physical rather than a mathematical requirement). The integrated flux-profile relations complicate the problem algebraically, although the relation between Z/L and RiB has the same simple form as for the gradient Richardson number relation,

in which

The quantities FM and FH are the integrated forms of the flux-profile laws and are given by, for unstable conditions (the Dyer-Hicks relations),

Approximate forms of these relations for -L>> Zo are,

For the Dyer-Hicks-Carson-Richards stable case, the relations corresponding to (4.20) are,

In the neutral case (), the relations for hs/L reduce to

It is difficult to see how (4.19) and (4.20) reduce to this simple relation in the near-neutral limit; accordingly, alternative, but equivalent, forms for FM, H are often used. This form is achieved by adding and subtracting Z-1 to the integrated profile relations, which, for the unstable case yields,

For the stable case, it is clear that FM,H reduces to (4.24).

To simplify the notation, this last system is usually written in the conventional form (as it is in PROGTN)

where the functionsmay be considered to be 'correction' terms that account for the departure from neutral. Bothand () are zero for and increase as |hs/L| increases. For conditions close to neutral, increase linearly; for strongly unstable conditions,increase logarithmically:

Eqns. (4.18 - 4.21) seem to show that hs/L is a strongly nonlinear function of RiB, at least in the unstable regime. Computational experience shows, however, that for RiB < 0 this is not the case and that a few iterations ar e sufficient to yield and adequate approximation to hs/L. In the extremely stable case near the critical Richardson number, hs/L is a very sensitive function of RiB. For such values the (downward) fluxes are quite small, and a close approximation to hs/L is unnecessary, especially since the relationship between hs/L and RiB is largely conjecture.

  1. Specific Forms of the Profile Laws and Choice for Selection of Certain Boundary Layer Constants
The exact functional forms for with their attendant constants is still a matter of active dispute among boundary layer meteorologists, as is the value of von Karman's constant, the most fundamental constant in turbulence theory (in fact, von Karman's constant may not be a universal constant at all, but may vary slightly with the slow regime). In any case, two differing values prevail: that of Businger, et.al. (1972), who give k = 0.35, and Dyer and Hicks (1974), who give k = 0.41. The reason for the differing values of the von Karman constants as well as the other free constants have been suggested by Wieringa (1980) to be systematic experimental error. The uncertainty in the current forms of , at least for the unstable case, is not likely to lead to very differing values of the surface fluxes.

The stable case is less settled. There is, first of all, the question of whether a critical Richardson number exists, and if it does, what its value is. Over the past several decades, the accepted value of Ric (assuming it exists) has fluctuated but has usually hovered about 0.2. All flux-profile laws, whatever their functional forms for small and moderate values of Z/L, are ultimately constrained by critical Richardson numbers provided become linear as . PROGTN provides a particular example of such functions, as proposed by Carson and Richards (1978). Although their proposed profiles are nonlinear for moderate values of Z/L, for larger values their return to linear. Their system therefore, falls prey to the fate of all linear profile laws: a critical Richardson number, albeit displaced. For our use in the MRF, their manipulation of the profile laws seems hardly worth the additional computational effort, since the final result is merely a shifted Ric, with no fluxes whatever for values exceeding Ric = 1.32.

It is planned that the MRF will have the Carson and Richards profiles replaced by a set for which . These new profile-flux relations do not require iteration to estimate hs/L, and have the property that the surface fluxes decrease smoothly as RiB increases.

  1. The Solution for the Obukhov Length from the Nonlinear Equations Using the Standard Newton Raphson Method
The determination of the Obukhov length (L) is basic to the calculation of turbulence quantities in the surface layer. The turbulent surface fluxes of momentum, sensible heat, and latent heat () are given by,

These relations, derived from (4.1-4.3), show that the fluxes depend upon the local gradients of U, Q, and q. As they stand, Equations 4.30-4.32 are not of immediate use in PROGTN, since PROGTN does not deal with local gradients but ratherDQ,U and q . The bulk difference relations (4.19-4.20) can be written as,

Substitution of (4.33-4.35) into (4.30-4.32) gives

In the above, r is the density of air, cp is the specific heat of air at constant pressure, and Lv is latent heat of vaporization. Variables CD, CH, and CQ are the transfer coefficients for momentum (the drag coefficient), heat, and latent heat, and are fundamental to PROGTN. Combining (4.33-4.35) with (4.36-4.38), we get the important result,

assuming FH = FQ; b = c = 1 (Dyer-Hicks 1970). The transfer coefficients CD and CH are non-negative and vary smoothly with hs/L. The coefficients increase with increasing instability and decrease with increasing stability. An order of magnitude of change in CD,H can be reasonably expected within a typical MRF forecast.

If L (or hs/L) is known, then FM,H, CD,H, and the surface fluxes are readily computed. The central problem is this: given DQ , U, Dq,, hs and Zo, use (4.18) - or a relation derived from (4.18) - to compute L efficiently. We have noted that (4.18) seems to be rather nonlinear because of FM,H, but computational experience proves the nonlinearity to be benign. We also noted that for the near-neutral case, (4.18) collapses to

This simple relation - unexpectedly - proves to be a fairly accurate approximation to (4.18) over a fairly wide range of unstable cases (unless Zo is very large).

To help substantiate this claim, we examine the behavior of a general fM,H,

for small values of - Z/L. We have


in which. By substituting (4.46-4.47) into (4.18) and making some approximations, we get

For the Dyer-Hicks relations (unlike the Businger relations) gM= gH=16. Thus, the RiB2 term vanishes, and (4.12), fortunately, proves to be accurate to (at least) order RiB2. Calculations show that (4.12) tends to consistently overestimate -hs/L somewhat for most reasonable values of RiB and Zo.

To get an accurate approximation to hs/L for cases of strong instability and large Zo, an iterative solution must be invoked. For this the Newton-Raphson method is used to get rapid, quadratic convergence. Iteration is also required for much of the stable regime.

We first form the function G(zs, z0) in which hs/L. We must solve

The Newton-Raphson iteration to F (x) = 0 is

Thus, the n+1 iterate for z is

Expressions for over the entire stable-unstable range are rather tedious and will not be reproduced here. Convergence occurs within several iterations, except for seriously anomalous cases.

  1. The Introduction of Virtual Potential Temperature into the Surface Layer Equations
Upon introducing the Bulk Richardson number in (4.17), we noted that this definition is strictly valid only for 'dry' conditions in which water vapor does not significantly alter the density of air. To generalize RiB, we account for the effect of moisture by selectively introducing virtual potential temperature, qv, in place of q.

The turbulent kinetic energy equation (see, for example, Brutsaert, 1982) requires that the kinematic heat flux, wq, be replaced by wq+ 0.61 q wq in the buoyant production term. This has the effect of introducing virtual potential temperature and forces us to expand the definition of the 'dry' Obukhov length:

The kinematic heat flux Qo must be replaced by

The introduction of moisture does not change the original simple form of the relation between hs/L and RiB provided: (1) nq in RiB are replaced by nqv and; (2) Q(Z/L) =H(Z/L); (3) the roughness lengths for heat and humidity are equal, although they need not equal the roughness length for momentum. The assumptions that Q=H and ZoQ=ZoH are standard for most numerical models, including the MRF.

We begin with the defining relation for qv,

from which we have

The higher-order term will be omitted, thus giving,

or, equivalently,

We recall that

which, when substituted into yield,

The potential temperature difference in (4.59) can be approximated by

The two e-terms are generally of very different magnitude: eqDqv can be neglected in comparison to . This leads to

This relation for , the moist analogue of

contains the perturbation term, (bFQ - c FH)/FHFQ. The SL relations used in PROGTN set b = c = 1; further, both the Dyer-Hicks and Businger (b = c = 0.74) relations assume that

are equal. This equality subsumes both and Z0Q = Z0H . The presumed equality of FH of FQ conveniently causes the moisture perturbation term in to vanish. This yields


The form of the basic relations therefore remains unchanged, except of the replacement of by . The calculation of L and FM,H,Q proceed as before [N.B.: FH =-rCpCHUDq, as before; is not replaced by U. PROGTN requires the sensible rather than the virtual heat flux in the surface energy balance.]

Given the lack of compelling evidence, we have not explored the possibility that ZoH  ZoQ or fH  fQ although there is fairly strong evidence that, in general, ZoM > ZoH,Q (Garratt, 1978, Garratt and Francey, 1978). Theory and experiment show that ZoM should be much greater than ZoH or ZoQ over rough areas. The physics of momentum transfer involves viscous shear as well as form drag generated by sharp local pressure gradients. This process is considerably different from that of heat and moisture transfer that involve molecular diffusion as well as turbulent diffusion. Garratt (1973) has suggested that ZoM  ZoH. On the other hand, there is little evidence to suggest that ZoH  ZoQ or foH  foQ. For this reason, as well as for computational efficiency, we assume the perturbation term vanishes and that (4.66) can be used in PROGTN.

The static stability of the surface layer, therefore, depends upon the vertical gradient of the virtual potential temperature rather than the (dry) potential temperature. This distinction can lead to the seemingly paradoxical situation in which the (kinematic) sensible heat flux (-) flows in the opposite direction to the virtual heat flux (-). No paradox exists, in fact, since the fluxes play differing roles. The sensible heat flux, along with latent, ground, and radiation fluxes, enters into the conservation of energy relation that determines the surface temperature. The virtual heat flux, on the other hand, enters into the conservation of turbulent energy equation that determines the rate of change of turbulent kinetic energy. Along with the virtual turbulent heat flux, there are also terms representing shear production, viscous dissipation, and redistribution of turbulent energy.

Physical intuition suggests that opposing model fluxes of sensible and virtual heat should arise only for near-neutral conditions involving weak surface fluxes. Our model forecasts what this is not always the case. For a given U, , and , the functions FM,H,Q decrease with increasing Zo; concomitantly, the transfer coefficients CM,H,Q increase. For a mildly positive , and with large Zo and -, the sensible heat flux (downward) can be readily exceeded by the - 0.61 uoqo term, thus producing an upward (unstable) virtual heat flux. MRF forecasts demonstrate that these antiparallel heat fluxes arise most commonly over tropical land areas. In these cases, the negative sensible heat fluxes tend to raise the surface temperatures, while the positive virtual heat fluxes tend to increase the turbulent kinetic energies.

  1. The Computation of Surface Temperature
The surface temperature in the MRF is predicted over all solid surfaces. These surfaces include: soil, snow or ice lying upon soil; sea ice; snow lying upon sea ice. Temperatures of water surfaces are subject to rather mo re complex physical processes than solid surfaces and are beyond the current capacity of the MRF to forecast. Temperatures of water surfaces, important as they are for the calculation of evaporation over much of the globe, are treated as passive elements and uncoupled in the MRF from the prediction of solid surface temperatures. Water area surface temperatures are externally supplied. See Chapter 7 for more detail.

The surface temperature Ts(t) (hereafter, 'surface temperature' denotes the air-ground interface temperature over a solid surface at about Z = Zo) is influenced by the energy fluxes of short and long wave radiation and latent heat, ground heat transfer, and liquid-solid phase changes. The two fundamental techniques for computing surface temperature in NWP models are the energy balance and the temperature "marching" (or tendency) methods.

Energy balance techniques require the vector sum of the energy fluxes to vanish at the earth's surface at each time step. This stipulation leads to a nonlinear equation of the general form f(Ts) = 0, where each term of f(Ts) is an energy flux term that, except for the short wave radiation flux, strongly depends upon the surface temperature. Variables not directly related to Ts in an energy balance are treated as quasi-constants for a given time step.

Solutions to energy balance equations are commonly computed by using the Newton-Raphson itertive root-finding technique which requires the calculation of . The derivative of f(Ts) usually cannot be expressed completely analytically, since an energy balance equation usually contains terms in Ts that are not in closed form. Accordingly, must be approximated by the semi-linearization of f(Ts) or by finite difference approximations. Approximations can slow convergence. A simplified energy balance method has been used in the MRF in the past to initialize surface temperatures, after which they were 'marched' forward in time by a predictive equation.

Energy balance method in NWP models differ chiefly in their handling of the surface heat flux term (hereafter, to be denoted as the "ground heat flux", including cases in which the surface is snow or ice). Although the ground heat flux is generally relatively small compared to the above-ground fluxes during the day, it si the ground heat flux which largely balances the infrared surface heat loss at night. The ground plays an important role in absorbing and storing energy during the day for release at night. If is neglected, then the computed nocturnal infrared loss results in an unrealistically rapid surface temperature drop.

In some models, is equated to a specified fraction of certain above-ground fluxes. The approach simplifies the mathematical formulation of f(Ts). It is partially justified on empirical grounds and also on the grounds t hat an accurate estimate of is a chimera, due to the lack of sufficient reliable subsurface data - or complete absence of such data - over much of the globe. For example, Kasahara and Washington (1971) proposed that

for bare soil; Fuchs and Hadas (1972) suggested instead

in which Rnet is the net long and short wave radiation at the surface. Their analysis indicates that Rnet is better satisfied for wet soil than for dry. For dry soil, a phase shift apparently exists between and Rnet. Nickerson and Smiley (1975) proposed two values for r as the result of their analysis of the O'Neill data: 0.19 for downward Rnet and 0.32 for upward Rnet. These simplifications for are unlikely to be generally applicable, however, since the 'constants' of proportionally probably depend upon location, soil moisture, and time of day.

The simplifying limiting case, , was used by Manabe et al. (1974) and more recently by Mahrt et al. (1984) [later recanted: Mahrt et al. 1988)]. Nulling is equivalent to assuming that the surface is a perfect insulator with zero thermal conductivity. Mahrt et al. (1984) argued that a perfectly insulated land surface is a good approximation over much of the globe. This arises because of the poor heat conducting properties of accumulated vegetative debris. Their model experiments imply, however, that perfect insulation is generally a poor assumption, as was foreshadowed by Deardorff (1978).

Snow, an exception, is an excellent insulator: the heat conductivity for snow is more than an order of magnitude smaller than for (say) clay.

Differing soil types exhibit great variation in and . Novak and Black (1985) presented simulations in which the daytime maximum for wet sand (a good conductor) was 325 Wm2 while for dry peat (a poor conductor) was 25 Wm2. The surface temperature maximum for the peat exceeded that for sand by about 505C. In addition, the peat reached its maximum temperature several hours earlier than the sand. Similarly, Deardorff's simulations in which the ground was treated as a perfect insulator showed that the maximum calculated temperatures over a clay pasture, dry quartz sand, still muddy water, and the O'Neill soil were higher and were phase-advanced compared to the simulations with proper thermal conductive properties. Deardorff's nocturnal simulations also showed the expected unrealistic rapid temperature drop at night for perfect insulators . Grossly over-simplifying the calculation of or by neglecting altogether achieves a computational economy that degrades both simulations and forecasts.

The most widely used methods for computing surface temperatures in NWP models are those which solve time dependent forecast or "marching" equations for Ts Predictive equations for Ts are the most commonly used equations for calculating Ts in NWP models. The appeal of predictive equations lies partially in their abandonment of cumbersome iterative processes required by energy balance methods and partially in the intuitive physical reasonableness of the role of the flux forcing functions. That is, the direction of the individual fluxes (positive or negative) conforms to the intuitive notion of a rising or falling surface temperature.

PROGTN uses the more straight-forward (yet more time-consuming) of the two most commonly used predictive techniques: a finite difference solution of the subsurface heat transfer equation. The vector sum of the above-ground fluxes supplies the upper boundary condition. The lower boundary condition is the deep below ground temperature that is sufficiently deep to be immune to diurnal influences. In the current MRF, the depth of is 5 meters. varies from gridpoint-to-gridpoint but is independent of time (see the chapter on 'fixed-fields' ). The choice {;} is not clear-cut and partially depends upon the nature of the NWP model, e.g., mesoscale or GCM. The selection of {;} also depends upon subsurface thermal properties and the number, depth and separation of the subsurface levels.

The more frequently used method for predicting Ts, especially in mesoscale modeling, was devised by Bhumralkar (1975) and independently by Blackadar 91976) [hereafter denoted the BB method]. The BB method uses only two level ls, the air-ground interface and Z , but suffers from a similar ambiguity in {;} as does PROGTN.

Standard equations for heat transfer are used in the PROGTN and BB methods. Heat transfer within the ground can be predicted by

In these subsurface heat transfer equations, z is, by convention, positive downward; = density of the soild subsurface; = subsurface specific heat; = volumetric heat capacity; Ks = subsurface thermal diffusivity; subsurface heat flux at depth z [nb: continues to denote , also, ]. The lower boundary condition for the heat transfer equations is . The upper boundary condition is computed from the energy balance equations,

Note that the distinction between iterative energy balance techniques and techniques urging predictive relations is more apparent than real - the techniques are, of course, all expressions of the conservation of energy. The choice of one over the other id dictated by simplicity, computational economy, and numerical sensitivity (e.g., convergence rate for the energy balance, 'stiffness' for predictive differential equations).

In particular, the BB technique originates with an energy balance relation in which is approximated by an expression that contains the difference Ts - and a term proportional to dTs/dt:

Here the angular velocity of the rotation of the earth , , and , which has dimensions of length, defined by

Eqs. (4.73-4.74) are derived by solving (4.70-4.71) with (4.72) as a lower boundary condition and a sinusoidal upper boundary condition (diurnal, with period = :

The solution of this system is, witj transients removed,

The relation for [(4.73)] then follows. Eq. 4.76 shows that is the "e-damping depth" of a diurnal wave. The damping depth also plays an important role in understanding the behavior of PROGTN's prediction equation.

The BB formula for dTs/dt results from eliminating from the energy balance equation, thus yielding,

The net (positive upward) above-ground flux is denoted by ,

Note that Rnet is positive upward, as are H,Q.

In the absence (or cancellation) of the above-ground fluxes ('forces'), the BB formula tends to 'restore' Ts to the subsurface with a (e-damping) time scale . Hence, the BB relation is often referred to as the 'force-restore method'.

Simulation experiments suggest that the force-restore method - although derived form the stringent assumption that (z,t) is sinusoidal with a period is surprisingly robust and versatile. Force-restore predictions under varying physical conditions often compare favorably with the more accurate, but more complicated and more time and memory consuming finite difference solutions that use many computational levels within the subsurface.

Numerical experience also shows, however, that the force-restore method can also produce 'solutions' that are as unrealistic as the finite difference solutions computed in PROGTN that uses three predictive layers. One likely reason for failure is the 'stiffness' of the BB and PROGTN predictive differential equations that contain multiple time-scales of substantial disparity (minutes versus hours) between the individual forcing and force-restore terms. Stiffness is most acute for situations with strong sensible and (especially) latent heat fluxes. The study of simple, low-order nonlinear finite difference equations also reveals solutions that are similar to the strange solutions for Ts produced by PROGTN and by the BB methods: 'transcritical', 'flip', and 'fold' bifurcatons (Holden, 1986). That is, the erratic behavior of the solutions is inherent in the nonlinearity of the difference equations. Ill-behaved forecasts can also occur with highly insulating surfaces, that is, when .

As noted, diagnostic energy balance and the predictive equations originate from the conservation of energy. The conservation of energy can be expressed as a continuity equation for energy density,

in which is the energy density, and E is the energy flux vector. E is directed along the Z-axis and is above ground and below ground. Accordingly, the operator reduces to . If we integrate downward from the air-ground interface to a depth we have,

in which 0+ denotes a point at a vanishing small distant above the interface. The first term in (4.80) is the rate of change of (heat) energy in the sublayer . Denoting as the mean in at t, then we have,


The sublayer thickness in PROGTN is 5 cm. Its representative (average) temperature is taken as is defined in PROGTN to be Ts. Thus, the so-called interfacial temperature Ts is, in fact, a subsurface temperature. PROGTN uses

for the "surface" temperature prediction equation. The approximation is valid provided [see (4.74)].

The diagnostic energy balance equation (4.72) is recovered if the limit is taken in (4.79). In taking the limit of

we assume that has, at most, a finite (jump) discontinuity at the interface; thus, S1=0. If, however, has a delta-function discontinuity, then we have an interfacial energy source or sink .

Photosynthesis is a physical example of a sink or storage term. Source and sink terms are not negligible, but, with the exception of water phase changes, are beyond the scope of our consideration. We therefore assume that, in the absence of phase change, S1=0 and, therefore,

which is the diagnostic energy balance equation given by (4.72). Phase changes in the prediction equation for Ts are not included explictly in (4.83) but are taken into account over snow, ice or water if the predicted value of Ts exceeds or falls below the melting or freezing point.

In addition to the predictive equation for "surface" temperature, , PROGTN carries two additional subsurface predictive equations for T-1 and T-2 at depths Z-1= 10 cm and Z-2= 50 cm. The lowest level temperature at T-3(Z-3= 500 cm) is intended to be the lower boundary condition for the heat diffusion solution of {Ts, T-1, T-2, T-3}. [nb.: we use the notation {T-3, Z-3} and almost interchangeably; we also return to Z positive downward.] The geometry of the subsurface levels is shown in Fig. (4.1).

The finite-difference approximation for is motivated by (4.71) and Fig. (4.1),

Note in Fig. (4.1) that is located at Z = = 5 cm, but the midpoint between Ts and T-1 is Z=6.25 cm hence, is not an exact second-order difference approximation.

The semi-discrete approximations to are derived by letting T1 be the representative temperature in the layer Z1A < Z < Z1B ; T-2 in layer Z2A < Z < Z2B; T-3 the lower, unperturbed boundary condition: T-3 = = constant for Z > Z-3. The planes Z1A, Z1B, Z2A, Z2B are located at levels,

The fluxes at these levels are given by

The rates of change of T-1, T-2, T-3 are

These relations may also be written as

As in the case with the approximations for , in (4.92-4.94) are not exact centered differences.

To help keep the time-marched solutions for {Ts,T-1,T-2} numerically stable and well-behaved, the system (4.83), (4.92-4.97) is integrated implicitely; therefore, a nonlinear term is semi-linearized by

For example, the Boltzman term, , becomes

The CH,Q transfer coefficients are troublesome since there is no closed-form relationship involving Ts. Computations show that transfer coefficients can vary sharply from tie step to time step. This can cause Ts to become, while not numerically unstable, unphysically oscillatory: we call such solutions 'fibrillations'. Semilinearization, where practicable, is useful for holding fibrillations in check. FQ terms particularly need semilinearization since they contain the surface saturation humidity, q*(Ts). For high Ts over the tropics, q*(Ts) can change strongly from . Thus, we have approximately,

in which es is the saturation vapor pressure, mv, the 'molecular weight' of vapor
(18.016 g mol-1), and R is the gas constant (8.314 erg mol-1 k-1.

The time-difference relations for T1, T2 are

As noted, the spatial differences are offset by small distances and are not true centered differences. The difference scheme for Ts is

Solar (shortwave) and downward longwave radiation are denoted by S and .

Although the finite difference relations above are cast in as implicit a form as possible, there are a few inconsistences. The sensible heat flux term contains , since it is not yet possible to know . In essence, is assumed to remain constant from time step n to n+1. This means that atmospheric and subsurface temperatures are not marched forward in time as a unit, but rather as two independent atmospheric and soil subsystems. This compromise doubtless helps create patchy areas prone to numerical irregularity, but such numerical problems can be largely reined in by time filtering. A holistic solution does not appear to be worth the computational effort at this time.

The choice of merely three levels within the ground, although physically overly parsimonious, is computationally a prudent choice; the soil subsystem reduces to the simple matrix system,

The elements Aij over snow-free land are:

A12 =

A13 =

A21 =

A22 =

A23 =

A31 =

A32 =

A33 =

D1 =

D2 =

D3 =

The use of only three active soil levels allows a simple application of Cramer's rule for (Ts,T1,T2)n+1:




Despite the intuitive appeal and simplicity of this approach, we have found that it has several problems. The first is the apparent irrelevance of the deep-soil boundary condition () imposed at a depth of five meters. The e-1 penetration depth of the diurnal wave (8.5 cm) into the "universal" soil is a very small fraction of the 500 cm depth of the lower boundary condition. Numerical experiments show that Z-2 = 50 cm becomes, in fact, the effective lower (diurnal) boundary condition. This leaves two coupled ordinary differential equations for Ts and T-1. It can be shown by numerical experiments and analysis that this 2x2 system is not capable of faithfully reproducing even a simple diurnal sinuosidal surface temperature wave, given a sinusoidal net ground flux.

The second difficulty, which requires time smoothing to control, is the sporadically labile surface temperature forecast. When the unsmoothed temperature forecast become ill-behaved, Ts(t) does not exhibit the explosive instability that is the signature of an unstable difference scheme. Rather, the Ts forecast begins smoothly and reasonably, and then suddenly begins to oscillate - almost chaotically - producing unreasonable sensible and latent heat fluxes. These unrealistic fluxes are a once a cause and effect of the prediction scheme's poor behavior and are produced, in part, by the highly nonlinear dependence of the heat and moisture transfer coefficients upon the static stability of the surface boundary layer. The energy transported by evaporation in the model can easily (and falsely) overwhelm the incoming solar energy.

The first of these basic difficulties is a direct consequence of the number and placement of the soil levels, and also the soil thermal properties which are assumed to be everywhere constant (: 'universal' soil). The soil properties are, in cgs units : soil density == s = 1.5 gcm-3; specific heat = cs = 0.32 cal g-1 deg-1; heat diffusivity =Ks=3.0x10-3 cm2s-1; volumetric heat capacity == 0.48 cal deg-1cm-3. We can anticipate the lack of influence of upon by several means. First, we examine (4.96).

We see that the time scale associated with heat transfer from level Z-3 to Z-2 is roughly given by time scale(3,2)3.4x107sec=391 dy. Now let us assume {Z-2, T-2} is made the lower (diurnal) boundary condition with T-2 = 0. The time scale for the heat transfer from level Z2 to Z1 is is very approximately time scale(2,1) 3x106 or 35 dy. Thus, the influence time scale for level Z-3 from level Z-2 is about a year, while about a mobth is required for level Z-2 to influence Z-1. The choice of T-3 is irrelevant.

Another approach that suggests the irrelevance of T-3 is to examine the analytical solution for the surface temperature,

We see that the surface temperature depends upon the history of the above-ground flux, , and the initial soil temperature profile. We wish to determine the depth of the initial temperature that fails to influence the surface temperature within a specified time. More concretely, if we ask: for what depth does the exponential in (4.106) fall below, say, 0.01 for t = 30 dy, then find that Z = 378 cm. By contrast, we see by using the same conservative criterion that a temperature at a depth Z = Z-2 = 50 cm is felt in only 12.6 hr. The selection of at = 500 cm would appear to have only a modest effect upon a 30 day surface temperature forecast; setting =50 cm would seem to be a more prudent choice. The irrelevance of T-3 is, therefore, not a simple numerical fortuity created by the form of the finite difference scheme.

The second difficulty is the sporadically irregular behavior of the surface temperature in time. These irregularities ('fibrillations') can also occur with the BB scheme and evidently arise when and the time scales associated with the radiation, ground , and turbulent heat fluxes become widely disparate and range between minutes to hours. When the flux time scales are all significantly longer than the time step used in marching Ts, and strong fibrillationss generally do not appear and the integration proceeds smoothly.

The wide range of time scale of the turbulent fluxes of sensible and latent heat result from the nonlinearity of CH(= CQ) and q*. The heat transfer coefficient CH varies sharply for relatively small deviations from neutral; moreover, CH is strongly assymetric with respect to . For unstable , CH increases rapidly compared to the newtral value and approaches in the free convection limit. (NB: the curent version of PROGTN does not quite folow this power law; the diffusion process above the surface layer in the MRF does). For stable , CH decreases rapidly to zero for . Over hot, unstable, moist surface CH, and q*(Ts), can become quite large and can cause an acute drop in Ts over the time interval . This extreme daytime temperature drop forces the MRF's SL to become strongly stable and begins to create cycles of overshooting: stable, unstable, stable, etc. Periods of fibrillations are intermingled with period of calm.

A considerably smaller time step would allow Ts to 'catch up' with the air and ground temperature but t << 10 min is not practicable. Use of a time step in which is the shortest of the time scales associated with the sensible, latent, and ground heat fluxes) can greatly burden the finite difference scheme. Recall that the finite difference scheme for Ts is given by (4.106). The time scale for the ground heat flux between Zs and Z-1 is

The for the ground heat flux exceeds one hour and, accordingly, should be easily handled with the MRF's . The time scale IR of the long wave term is approximately

or about 10.3 hr. As expected, the long wave radiation poses no numerical problem. The turbulent terms are quite a different matter. For the sensible heat flux term, we have,

for CH = 10-2, while for the latent heat flux q 1.9 min. For CQ = CH= 10-2, Ts = 40oC and qs = 15 gkg-1. The transfer coefficient is typical of reasonably unstable conditions with Z0=0.1m. We see that the time scales are both short and that we should expect computational irregularities. Even if there are no evident fibrillations, forecasts over hot, moist regions are regarded with suspicion. Computational experience shows that for smooth spurious 'solutions' to the marching equation can be generated which bear little relationship to the physical solution. Such spurious solutions may well be examples of the bifurction phenomenon that commonly arises in the study of nonlinear difference equations. We cite two examples:



In both examples we shall assume that the point T = 0 is the 'physical' equilibrium point (equilibrium point fixed point critical point Tn+1 = Tn). For (4.122), T = 0 is a stable solution provided -2 < A < 0. However, there is a second equilibrium solution, Tn = Tn+1 = -A, that is a stable solution for 2 > A > 0. Thus, when the physical solution becomes unstable, the iterations are attracted to the spurious solution. In the second example, the physical solution becomes unstable for A > 0. When A > 0, a curious behavior begins: a 'flip' bifurcation which creates '2-cycle' equilibrium points,



which are stable if (1 - 2A)2 < 1/ An initial To near zero is repelled by T = 0 and attracted to T = + A, at which points successive iterates flip sign [see Holden, 1986 for a fuller discussion of (4.124) and numerous of her examples].

These examples show that the behavior of even relatively simple nonlinear differences schemes cannot be satisfactorily explained in terms of linearlized schemes. For these reasons, the use of a carefully tested smoothing technique is indispensable for the simulation of surface temperatures on a global model with extensively varying surface conditions.

  1. The Prediction of Surface Temperature in the Presence of Snow or Ice
Snow or sea ice that is melting substantially changes the prediction of surface temperature. The modeling of melting is one of the oldest portions of the MRF's surface physics; accordingly, melting is one of the most simplistic portions of PROGTN. Some of the simplifying assumptions are:

1) Melting takes place only in an infinitesimally thin layer at the air-ground interface regardless of the subsurface temperature; (2)Melting occurs only when the predicted Ts exceeds 273.16K for snow or 271.2K for sea ice; (3)Meltwater is not allowed to percolate into snow and refreeze at a lower, colder depth with the accompanying heat of fusion transfer and changed temperature at that level; (4)Meltwater percolates immediately into the ground (without refreezing) where it increases the volumetric soil moisture content, provided the 'field capacity' is not exceeded; (5)If, during a given time step, all the snow cover is melted, then any remaining heat goes into increasing Ts. Other subsurface temperatures are unchanged; (6)The water depth of snow is 1/10 of the snow depth; (7)Sea ice meltwater puddles inertly on the surface of the sea ice; (8)The subsurface heat flux through sea ice is inversely proportional to an assumed constant sea ice thickness of 2 meters; and (9) Compaction of snow is not allowed.

These rather primitive and restrictive assumptions can perhaps be justified by appealing to the equally primitive knowledge we have of snow cover, snow depth and sea ice cover, salinity and depth. Without improved knowledge , there is little reason to upgrade the simulation of melting.

Consider now the following conditions: If the Ts prediction equation yields as predicted Ts that is lower than the melting point of snow, then the computation continues uninterrupted. If, however, Ts >273.16K, then the computation


is replaced by

in which is the vector sum of the above-ground fluxes, is the subground heat flux just below the interface, and is the heat of fusion of ice: 79.9 cal g-1. The quantity

is the energy per unit area per unit time that is devoted to melting snow, an isothermal process at 273.16K. If is the heat of fusion, the is the 'material flux', so to speak, of the melted snow: mass per unit area per unit time. The mass is, of course, the water mass which can be approximately converted to physical snow depth by multiplying by the conversion factor of 10. The water-equivalent depth that is melted in time

t is t, assuming that all the snow is not melted during the time

t. If the snow is not completely melted, then the remaining snow depth (that is, the water equivalent depth) is,


The new Ts 273.15K. If new is negative, that is, there is more heat energy available than is needed to melt the snow, then the computation is again interrupted. We first write,


in which s is the mass flux required to melt the remaining snow during the time t. The associated energy flux is Sold/ t. The remaining energy flux available is:


This available energy is used to raise the temperature Ts of the interface slab or thickness ,


Ts is the new (elevated) temperature. In the case of sea ice, we assume that sea ice is two meters thick everywhere and that = 2 m.

  1. Evapotranspiration Over Land
    1. Introduction
Since the recent studies of Walker and Rowntree (1977), Warrilow (1986), Shukla and Mintz (1982), Yeh et al. (1984) and others, it is becoming apparent that numerical weather prediction (NWP) models are very sensitive to the parameterization of the surface exchange processes at the atmosphere - land interface. The development of the daytime planetary boundary layer (PBL) is strongly dependent upon the parameterized surface sensible and latent heat fluxes. Convective precipitation over land is also sensitive to the soil - moisture / surface - evaporation parameterization.

While the parameterization of the sensible heat flux in most models is of the simple bulk - aerodynamic form based on the temperature difference between the soil surface and the lowest layer air, the parameterization of latent heat transport is complicated by the degree of wetness of the soil and the presence of plants. Most of the NWP models developed in the past utilized a simple "bucket" method pioneered by Manabe (1969). In this method there is a bucket at each grid point and evaporation is reduced from a "potential" value by the ratio of soil water in the bucket (w) and a specified field capacity value (ws). In addition, the potential evaporation is evaluated assuming that the soil is saturated at the model calculated "skin" temperature. It has been pointed out by Dickinson (1983) that this method cannot realistically model the complex biosphere processes.

Toward improvement of the biosphere parameterization, there are some fairly complicated models (e.g. Dickinson, 1984; Sellers et al., 1986) that have been developed recently. Less complicated models such as Pan and Mahrt (19 87) still require several levels of soil moisture and a plant canopy. There is little doubt that we need better understanding of the atmosphere - biosphere interactions. Interactive feedback from the biosphere is probably as important as feedbacks from deep cumulus and from the ocean. The level of complication needed to model this mechanism is, however, a decision that atmospheric modelers must make. This type of question has been asked concern ing other feedback mechanisms as well. A notable example is the choice of the cumulus parameterization schemes. A point against a complex biosphere parameterization scheme is the lack of coordinated atmosphere - biosphere data to validate the many parameters in these schemes. Such experiments as the Hydrological Atmospheric Pilot Experiments (HAPEX) and the International Satellite Land - Surface Climatology Project (ISLSCP) experiments are only beginning to make such data available and should help in the development of such models in the future.

For short- and medium- range weather forecasts, the benefit of a detailed biosphere model is not immediately obvious. Surface evaporation interacts with radiation, clouds, and convection and, yet, each of these processes is currently only crudely parameterized. At the National Meteorological Center (NMC) in the United States of America, we took the approach of a simple surface scheme at the level of complication comparable to the other parameterization schemes in the present generation of models. The Penman - Monteith method (hereafter referred to as PM) recently implemented in the medium - range forecast (MRF) model follows the above philosophy and will be described below.

  1. Method
In the simple bucket method (Manabe, 1969) formerly used in the MRF, latent heat flux from the land surface is parameterized using bulk - aerodynamic formula :

where is the soil moisture availability parameter,Ep is the potential evaporation, is the density of air in the first model layer, L is the latent heat constant of vaporization, Ch is the turbulent exchange coefficient, V is the wind speed of air in the first model layer, qs(Ts) is the saturation mixing ratio at the surface temperature Ts (sometimes referred to as the skin temperature), and qa is the mixing ratio of air in the first model layer. The factor, , is defined as the ratio between the soil moisture content (w) and the field capacity (ws) and has value between 0 and 1. The key assumptions of this formulation in addition to the bulk - aerodynamic ex change concept are:

1) the use of a single parameter, , to simulate the reduction of evaporation when the soil and the vegetation are under stress, and 2) the use of saturation mixing ratio at the 'skin' temperature as the surface mixing ratio.

Even when the soil is wet, there still exist some stomatal resistance in the plants that will reduce the evaporation from the potential value. When the soil is dry, the use of the surface temperature to estimate the surface mixing ratio is too high and will result in an overestimation of the potential value. In both situations, the simple bucket model will overestimate the actual evaporation. As pointed out by Mahrt and Ek (1984), the potential evaporation should really be defined using a temperature different from the surface temperature; a temperature that is compatible with the Penman formulation of the potential evaporation. The Penman potential evaporation formulation can be further modified to include the effect of stomatal resistance for vegetation and one can define a potential evapotranspiration following Monteith (1965). For the modeling of soil water stress situation for the atmosphere - biosphere interface, we have decided to retain the framework of the bucket model but modify it to overcome these weakness. We will describe the methodology in some detail in the remainder of this section.

We will follow the work of Mahrt and Ek (1984) and define potential evaporation as the evaporation that can be realized if the soil is completely wet given the same environmental conditions ,i.e., net radiative flux and ground heat flux. The key point implicit in this definition is the existence of a different skin temperature under the saturated soil condition ( see also Sud and Fennessy, 1981). By allowing the upward longwave radiative flux t o depend upon the skin temperature, we can write the surface energy balance under the saturated soil condition ( =1) as:


where the terms on the left-hand-side are, respectively, the net shortwave radiative flux ( is the albedo), the downward longwave radiative flux, and the upward longwave radiative flux ( is the Boltzmann constant). The terms on the right-hand-side of the equation are the ground heat flux, the sensible heat flux, and the latent heat flux. In contrast, the actual surface energy balance is defined as follows:


where the last term is obtained from Eq. (4.131). The difference in the skin temperature, Ts - Ts', can be larger than 10 K when the soil is dry. In the original simple bucket method, we ignore this difference by using Ts to calculate Ep and thus tend to overestimate the latent heat flux. As a result, there is strong moistening in the lowest model layers during the first 12-24 hours of forecast. At the same time, temperature in the lower layers of the model also tend to be lower during the daytime. This leads to an unrealistic Bowen ratio (the ratio of the sensible and the latent heat fluxes) even over regions of fairly wet forest land (e.g. the Amazon region). We interpret the latter as an additional problem of not parameterizing the presence of vegetation.

In Eq. (4.131), the single unknown variable is the skin temperature Ts'. Mahrt and Ek (1984) rederived the Penman potential evaporation formula starting from (4.131) (neglecting the effect of skin temperature on the upward longwave radiative flux and the ground heat flux as done by Penman, 1948) to obtain


where Rnet is the net radiative fluxes,


Troen and Mahrt (1986) included the effect of the lower skin temperature in Rnet to obtain



To derive (4.133), the linearized form for qs(Ts'):

is used. It should be noted that there is a dependence of Ts in the ground heat flux. We are presently ignoring this effect as it would significantly complicate the potential evaporation formulation. Monteith (1965) and others suggested that, in the presence of plants, the stomatal resistance (rs) must be included. While we follow the idea of a stomatal resistance, we will use a minimum stomatal resistance (defined as the resistance under no water stress) to define a potential evapotranspiration. The unusually small predicted Bowen ratio over tropical rain forest regions using the original bucket method suggests that we should include this effect. The resulting potential evapotranspiration rate when we apply the additional stomatal resistance to the latent heat flux in (2) is shown below:


Finally, the latent heat flux is now defined using the bucket concept as:


In this manner, the water stress in the stomate is parameterized in the same form as the bare soil simple bucket method. The value of stomatal resistance under no water stress is currently set at 60 s m-1.


Baumgartner, A., H. Mayer, and W. Metz, 1977: Weltweite verteilung des rauhigkeitsparameters Zo mit anwendung auf die energie dissipation an der erdoberflache. Meteor. Rdsch., 30, 43-48.

Bhumralker, C. M., 1975: Numerical experiments on the computation of ground surface temperature in an atmospheric general circulation model. J. Appl. Meteor., 14, 1246-1258.

Brutsaert, W. H., 1982: Evaporation into the atmosphere. D. Reidel Publishing Company, Dordrecht, 299 pp.

Businger, J. A., Wyngaard, J. C., Izumi, U., and Bradley, E. F., 1971: Flux-profile relationships in the atmospheric surface layer. J. Atmos. Sci., 28, 181-189.

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

Carson, D. J., and J. R. Richards, 1978: Modelling surface turbulent fluxes in stable conditions. Boundary Layer Meteor., 14, 67-81.

Deardorff, J. W., 1978: Efficient prediction of ground surface temperature and moisture with inclusion of a layer of vegetation. J. Geophys. Res., 83 (4), 1889-1903.

Dickinson, R. E., 1983: Land surface processes and climate-surface albedos and energy balance. Advances in Geophysics, 25, 305-353.

Dyer, A. J., 1974: A review of flux-profile relationships. Boundary Layer Meteor., 7, 165-184.

Fuchs, M., and A. Hadas, 1972: The heat flux density in a non-homogeneous bare loessial soil. Boundary Layer Meteor., 3, 191-200.

Garratt, J. R., 1978: Transfer characteristics for a heterogeneous surface of large aerodynamic roughness. Quart. J. R. Meteor. Soc., 104, 491-502.

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

Garratt, J. R., and R. J. Francey, 1978: Bulk characteristics of heat transfer in unstable baroclinic atmospheric boundary layer. Boundary Layer Meteor., 15, 399-421.

Hanson, J. E. and T. Takahashi, 1984: Modeling evapotranspiration for three-dimensional global climate models. Climate Processes and Climate Sensitivity. American Geophysical Union, Geophysical Monograph, 29, 58-72.

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

Kasahara, A., and W. M. Washington, 1971: General circulation experiments with a six-layer NCAR model, including orography, cloudiness, and surface temperature calculation. J. Atmos. Sci., 28, 657-701.

Mahrt, L., and M. EK, 1984: The influence of atmospheric stability on potential evaporation. J. Clim. Appl. Meteor., 23, 222-234.

Manabe, S., 1969: Climate and the ocean circulation: I, The atmospheric circulation and the hydrology of the earth's surface. Mon. Wea. Rev., 97, 739-774.

Manabe, S., D. G. Hahn, and J. L. Holloway, Jr., 1974: The seasonal variation of the tropical circulation as simulated by a global model of the atmosphere. J. Atmos. Sci., 31, 43-83.

Mellor, G. L., and T. Yamada, 1982: Rev. of Geophys. and Space Phys., 20, 851-875.

Miyakoda, K., and J. Sirutis, 1986: Manual of E-physics. Manuscript, GFDL, Princeton, New Jersey.

Monin, A. S. and A. M. Obukhov, 1953: Basic turbulent mixing laws in the atmospheric surface layer. Trudy Inst. Teoret. Geofiz, Akad. Nauk. S.S.S.R., 24, 163-187.

Monteith, J. L., 1965: Evaporation and environment. Symp. Soc. Exp. Biol., 19, 205-235.

Nickerson, E. C., and V. E. Smiley, 1975: Surface energy budget parameterizations for urban scale models. J. Appl. Meteor., 14, 297-300.

Novak, M. D., and T. A. Black, 1985: Theoretical determination of the surface energy balance and thermal regimes of bare soils. Boundary Layer Meteor., 33, 313-333.

Obukhov, A. M., 1946: Turbulence in an atmosphere with non-uniform temperature. Trudy Inst. Teoret. Geofiz. Akad. Nauk. S.S.S.R., 1, 95-115 (English translation: 1971, Boundary Layer Meteor., 2, 7-29.

Pan, H. -L., and L. Mahrt, 1987: Interaction between soil hydrology6 and boundary-layer development. Boundary-Layer Meteorolo., 38, 185-202.

Penman, H. L., 1948: Natural evaporation from open water, bare soil, and grass. Proc. Roy. Soc. A., 193, 120-195.

Sela, J. G., 1980: Spectral modeling at the National Meteorological Center. Mon. Wea. Rev., 108, 1279-1292.

Sellers, P. J., Y Mintz, Y. C. Sud, and A. Dalcher, 1986: A simple biosphere model (SiB) for use within general circulation model. J. Atmos. Sci., 43, 505-531.

Shukla, J., and Y. Mintz, 1982: Influence of land-surface evapotranspiration on the earth's climate. Science, 215, 1498-1501

Sud, Y. C., and M. Fennessy, 1981: A study of the influence of surface albedo on July circulation in semi-arid regions using the GLAS GCM. J. Climatology, 2, 105-125.

Treon, I., and L. Mahrt, 1986: A simple model of the Atmospheric boundary layer; sensitivity to surface evaporation. Boundary-Layer Meteoeorl., 37, 129-148.

Walker, J., and P. R. Rowntree, 1977: The effect of soil moisture on circulation and rainfall in a tropical model Quart. J. Roy. Meteor. Soc., 103, 29-46.

Warrilow, D. A., 1986: The sensitivity of the UK Meteorological Office atmospheric general circulation model to recent changes in the parameterization of hydrology. ISLSCP, Proceedings of an international conference. ESA, Paris, France, 143-150.

Wierginga, J., 1980: Representativeness of wind observations at airports. Bull. Amer. Meteor. Soc., 61, 962-971.

Yaglom, A. M., 1977: Comments on wind and temperature flux-profile relationships. Boundary Layer Meteor., 11, 89-102.

Yeh, T. -C., R. T. Wetherald, and S. Manabe, and S. Manabe, 1984: The effect of soil moisture on the short-term climate and hydrology change - A numerical experiment. Mon. Wea. Rev., 112, 474-490.


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.123) - (4.125) 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.140 - 4.141 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.195 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.