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.
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/L
2.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 Ri
Ric.
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 Z
100Zo for heat and (presumably) moisture. This disconcerting
result implies that for Zo
2m, 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.
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.
, 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
and
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 or
must 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,
and
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.
, 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 functions
may be considered to be 'correction' terms that account for the
departure from neutral. Both
and
(
) 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.
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.
) 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
and
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.
The turbulent kinetic energy equation (see, for example, Brutsaert, 1982) requires that
the kinematic heat flux, w q , be replaced by w q + 0.61 q w q 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
and
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.
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,
or
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:
(4.122)
(4.123)
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,
(4.124a)
(4.124b)
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) 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
(4.125)
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,
(4.126)
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,
(4.127)
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:
(4.128)
This available energy is used to raise the temperature Ts of the interface slab or thickness ,
(4.129)
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.
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.
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:
(4.131)
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:
(4.132)
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
(4.133)
where Rnet is the net radiative fluxes,
and
Troen and Mahrt (1986) included the effect of the lower skin temperature in Rnet to obtain
(4.134)
where
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:
(4.135)
Finally, the latent heat flux is now defined using the bucket concept as:
(4.136)
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.
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.
.
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:
where
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
or
and
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,
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,
and
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
and
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
and
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.
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.