Agricultural and Forest Meteorology, 47 (1989) 259-271 Elsevier Science Publishers B.V., Amsterdam - - Printed in The Netherlands
259
INVESTIGATIONS WITH A H I G H E R - O R D E R C A N O P Y TURBULENCE MODEL INTO MEAN SOURCE-SINK LEVELS AND BULK CANOPY RESISTANCES
KYAW THA PAW U
Department of Land, Air and Water Resources, University of California, Davis, CA 95616 (U.S.A.) TILDEN P. MEYERS
NOAA, Air Resources Laboratory, Atmospheric Turbulence and Diffusion Division, P.O. Box 2456, Oak Ridge, TN 37831 (U.S.A.) (Received October 17, 1988; revision accepted March 21, 1989)
ABSTRACT Paw U, K.T. and Meyers, T.P., 1989. Investigations with a higher-order canopy turbulence model into mean source-sink levels and bulk canopy resistances. Agric. For. Meteorol., 47: 259-271. The legitimacy of using bulk aerodynamic and canopy resistances in surface energy budget equations is examined in this paper. Specifically, the variation of the effective source heights for momentum, water vapor, and heat is analyzed by estimating zero plane displacements for a modeled soya bean canopy, when model leaf surface resistances are changed, causing varying stability conditions. Also, the changes in bulk canopy resistances are examined as a function of the leaf surface resistances. The model used is a linked soil-plant-atmosphere model based on higherorder closure principles. The model confirms previous reports of major soil contributions to evapotranspiration under certain conditions. The erratic behavior of the zero plane displacements for water vapor and heat as a function of leaf surface resistance demonstrates that the concept of a single effective source-sink height is not easily applied to plant canopies. The zero plane displacement for momentum was found to be consistent with previous results, independent of leaf surface resistance. Canopy resistance changes qualitatively match leaf surface resistance changes, but quantitative differences can be large. The canopy resistance changes less quickly than the leaf surface resistance even when soil evaporation is not a factor.
INTRODUCTION The estimation of evapotranspiration from vegetation and the soil has long been an issue in plant biometeorology. In the simplest models, the latent ene r g y f l u x d e n s i t y a s s o c i a t e d w i t h e v a p o t r a n s p i r a t i o n is s o l v e d w i t h a n e n e r g y budget. Within such models, the water vapor flux must pass through two bulk resistances: the aerodynamic resistance from the atmosphere to the height of
260
the effective source-sink of water vapor and heat, and the canopy resistance representing the physiological and soil-based resistance to water transport to this source-sink height. Sensible heat flux is assumed to pass through only the aerodynamic resistance. The zero plane displacement has been interpreted as the height of the effective sink for momentum (Thom, 1971 ) and has also been used in the same way for heat and mass fluxes (Brutsaert, 1982). Jackson {1981 ) showed that Thorn's (1971) equation for the zero plane displacement could be derived for flows when buoyancy influences were neglected. Brutsaert (1982) presented a theoretical methodology for estimating differing values of the roughness length for heat and mass by assuming that the zero plane displacement was the same for heat, mass, and momentum. Using experimental data and similarity arguments, others have noted that because the source-sinks of momentum, heat, and mass flux are usually at different heights, a different zero plane displacement should be identified for each flux (Garratt, 1978; Hicks et al., 1979; Raupach et al., 1979). Some of these researchers have also explained that measurements within the roughness layer can lead to erroneous zero plane displacements estimates (Raupach et al., 1980). The issue of using the same zero plane displacement and roughness length for momentum, mass and heat fluxes is important because these parameters are frequently employed to derive bulk transfer coefficients such as the aerodynamic resistance. This resistance and the bulk canopy resistance are used in energy budget and Penman-Monteith equations (Brutsaert, 1982; Choudhury et al., 1986). Such equations are included for the surface boundary conditions in many regional and general circulation models (Sellers and Dorman, 1987). If a unique zero plane displacement cannot be identified, the usage of bulk aerodynamic resistances will result in errors. In this paper we report the modification of a layered, plant-canopy energybudget model coupled with a higher-order turbulence closure solution scheme to address these issues. Our first objective is to examine if one can estimate the height of the effective source-sinks of momentum, heat and mass fluxes for a short crop (soya bean) canopy. This topic has been studied by Hicks et al. (1979) and Raupach et al. (1979) for tall (tree) canopies, but has not been considered in detail for shorter crop canopies. This is because the instrumentation used in the above studies is not easily placed in short, closed crop canopies. The second objective is to examine if bulk surface resistances (used in the Penman-Monteith equation, along with the aerodynamic resistance) can be estimated using leaf surface resistances. Similar work has been reported by others such as Cowan (1968), using first-order closure (K-theory) methods. More recently, first-order closure has been shown to be inappropriate for withincanopy turbulence descriptions (Shaw and Pereira, 1982; Meyers and Paw U, 1987).
261 MODEL DEVELOPMENT AND VALIDATION
We modified a model described by Meyers (1985) and Meyers and Paw U (1987) to include more of the atmosphere above the plant canopy. Our model is a higher-order closure method, where equations for first, second, third and fourth moments of the turbulent fields of wind velocity, temperature, and water vapor are solved. At the fourth moment level, closure of the equations is obtained by assuming quasi-Gaussian distributions (Meyers and Paw U, 1987). The turbulence model is linked to an energy budget and leaf physiology model, described briefly below. The original model, which used 41 nodes, was validated using the data of Baldocchi ( 1982 ) for soya bean (Glycine max) canopies. Validation of a version for forest canopies has been carried out by Meyers et al. (in press). For the version used in this paper, we extended the domain of the original model from 41 nodes to 61 nodes. This was to produce a smoother, more numerically stable output of the temperature, humidity and wind speed gradients near the soya bean canopy surface, and to keep the boundary of the domain at an acceptable distance from the nodes used in calculations in this study. The ground and plant canopy occupy the first 21 nodes, with the atmosphere described in the next 40 nodes. Within each plant canopy layer, the full energy budget is solved exactly (Paw U, 1987; Paw U and Gao, 1988). First-order approximate solutions to the energy budget, used previously (Meyers, 1985; Paw U et al., 1985; Meyers and Paw U, 1987) cannot be used here because the extreme physiological conditions we modeled (i.e. very high and very low surface and leaf resistances) could create large leaf-air temperature differences. Leaf surface resistances (stomatal resistances) were modeled by assuming they were inversely proportional to the visible radiation flux density (Meyers and Paw U, 1987). The soil boundary condition was approximated using a crude energy budget with a surface resistance rs1 for the soil, separating a perfectly moist soil interior from the atmosphere. The leaf and soil radiative and biophysical characteristics have been described in Meyers and Paw U (1987). The radiation within the plant canopy is described by a modified version of the Norman (1979) model, with an assumed spherical leaf distribution (other realistic leaf distributions could be used, but we found this distribution provided an adequate approximation). This relatively simple model was shown to work well earlier (Meyers and Paw U, 1987). The turbulent transfer of momentum, heat, and water vapor was solved directly using the higher-order closure scheme described by Meyers ( 1985 ), Paw U et al. (1985) and Meyers and Paw U (1987). The general solution relies upon applying Reynolds decomposition and then averaging combinations of the Eulerian Navier-Stokes equations of motion for each velocity component and mass/flux balance equations for each scalar (heat and mass). The fourth moments of the turbulent quantities resulting from the above analysis are as-
262
sumed quasi-Gaussian; the adequacy of this assumption for plant canopies is given by Meyers and Paw U (1987) and Shaw and Seginer (1987). In the present study, the model differs from that of Meyers and Paw U (1987) to include all effects of buoyancy terms, some of which were neglected originally when they were found to be small during the validation process. We included these terms because of the extremes in buoyancy effects induced by the extreme physiological conditions we modeled. METHODOLOGY FOR ZERO PLANE DISPLACEMENT ESTIMATION
The zero plane displacement was estimated in several ways. The "center-ofpressure" method (Thorn, 1971, modified by Shaw and Pereira, 1982) involves determining the centroid of the Reynolds stress gradient by integration ht-
f z Ouw - ~z- dz dc=h~
f
" Ouw O~z dz+uw(o)
(1)
0
where z is the height above the ground, ~ is the Reynolds stress, h~ is the canopy height, and d~ is the zero plane displacement. The centroid of the heat flux gradient dh¢ and the vapor flux gradient dye were calculated in a similar manner. A problem can exist when the gradients of the scalar fluxes change sign between the bottom and top of the canopy. Another way to estimate the zero plane displacement is to use assumed logarithmic similarity equations, and model outputs of the wind speed, temperature, and humidity profiles. We used the similarity equations of Brutsaert ( 1982 ) to estimate the zero plane displacements ds, dhs, and dw for momentum, heat, and vapor, respectively, for both stable and unstable conditions; the solution schemes for these equations are given in the Appendix. The roughness heights for momentum, Zo, and heat, Zohwere also determined using similarity equations given in the Appendix. BULK CANOPY RESISTANCE ESTIMATION
The feasibility of estimating a meaningful bulk canopy resistance was determined by using a "big-leaff' P e n m a n - M o n t e i t h equation with the higherorder closure model. This was done in several ways. The change in the estimated bulk canopy resistance rc was compared to the change in the input stomatal and soil resistances. The change in latent energy flux density was estimated as a function of changes in the input stomatal and soil resistances, and the estimated bulk canopy resistance.
263 The P e n m a n - M o n t e i t h equation for the latent energy flux density associated with evapotranspiration is
LE zl(Rn- G ) +flCp((~e)/ra
(2)
A+7[(ra+r~)/ra] where Rn is the net radiation flux density, G is the ground heat flux density (flux densities in W m - e ) , p is the air density, Cp is the mass specific heat, 5e is the vapor pressure deficit, 7 is the psychrometric constant (approximately 66 Pa °C-1 ), zl is the slope of the saturation vapor pressure curve (Pa °C-1 ), ra is the aerodynamic resistance (s m -1) and rc is the canopy resistance (s m - 1). This equation is based on the ability to identify a single plane which is the source-sink of heat and water vapor. In one test, the modeled value of latent energy flux density for the canopy and soil was used to solve for the canopy resistance re; the aerodynamic resistance ra was estimated by either assuming a neutral atmosphere (Rosenberg et al., 1983, with d--0.63 h~ and zo=O.13h~) or by using the mean wind speed at 2h~(Uh) divided by the friction velocity u* squared (no stability assumptions made). All of the individual stomatal resistances were then multiplied by a set factor, and the canopy resistance solved as described above was then compared with the original re. Presumably changes in r~ should match stomatal resistance changes. A second, more direct way to determine the feasibility of bulk canopy resistance estimates was to calculate a canopy resistance independently, and then substitute this into the P e n m a n - M o n t e i t h equation along with the aerodynamic resistance r a calculated as described above. The common technique of approximating the bulk canopy resistance by parallel weighting of the surface {stomatal) resistances by the leaf area index was used 1 rc~ = ~ LAIi i
(3)
rs i
where rsi is the surface resistance for the ith leaf class, LAIi is the L A I for that leaf class, and rcav is the average rc so calculated. A different version of rcav was estimated (rcavg), where the ground resistance rs] was included in eq. 3 defining the ground as a "big leaf" with an L A I of 1. The L E from the P e n m a n - M o n teith equation was estimated using eq. 3, and compared to the L E estimated from the higher-order closure model. Problems with weighting leaf area contributions have been discussed by Shuttleworth (1975, 1976) and Monteith (1977) using resistance theory and first-order closure, but no work has been done with a more accurate description of within-canopy turbulence. Usage of first-order closure is associated with
264 TABLE 1 Base values of i n p u t variables Reference height (Zr) Air t e m p e r a t u r e at zr H u m i d i t y at zr W i n d speed at z~ Net radiation flux density Solar radiation flux density Soil heat flux density
1.75 m 26.81°C 14.37 g k g - 1 3 . 2 7 m s -~ 399 W m - 2 584Wm ~ 19.0 W m -~
more serious qualitative and quantitative errors than usage of higher-order closure (Shaw and Pereira, 1982; Meyers and Paw U, 1987). MODEL INPUTS
The base set of inputs to the model are given in Table 1, and were taken from data gathered by Baldocchi (1982) in a soya bean crop. The base input to the stomatal physiology model is given in Meyers and Paw U (1987). The stomatal resistance of every leaf was increased or decreased in incremental steps by a factor of 2 or 4, from the resistances arising from the base input to the physiology model. At each stomatal resistance step, four soil resistances (0.0625, 100, 1000, and 15 000 s m -1 ) were used to simulate a very wet soil, drying soils, and a completely dry soil, respectively. The changes in the physiological nature of the canopy and the soil wetness varied the source-sink distribution of sensible heat and latent energy flux densities. At the same time, the stability of the atmosphere above the canopy changed because of the leaf surface and soil resistance changes. This allowed examination of the purely physiological effects of plants on the estimates of physical characteristics (i.e. the zero plane displacement and bulk surface resistance ). Model solutions included calculations of the sensible heat, latent energy and momentum fluxes and flux divergences, in addition to the profiles of temperature, humidity and wind speed. This allowed estimation of the zero plane displacement, roughness length, and canopy resistances as described earlier. RESULTS AND DISCUSSION
Zero plane displacements and roughness lengths On the x-axis of Fig. 1 (a), the leaf surface resistances are divided by their base values. The base values are determined by running the simulation using the base inputs (Table 1 ) and the physiological model used by Meyers and Paw
265 4.
3 ~" 2
3
oi
*
o
0
-~ o -~ '"i~ -"
o
1 0
4
~
oo
10
10 2
oo °
-0.1
a +3000
65
0.0 Zref//L
0.1
0.5 z~
2-
o
*
z~
m
a
~ o
i
*
0.0
N 0 L~
-1/ -0.1
0 19.-2oo 0.0 Zref/L
0/1
0.5 -0.10
-0.05
0.00 Zref/L
0.05
0.10
Fig. 1. Zero plane displacements and roughness lengths, with a soil surface resistance of 0.0625 s m ~ (wet soil). In all parts, stars depict displacements or lengths for momentum flux, circles for vapor flux, and triangles for heat flux. (a) Zero plane displacements determined by centroid method (Thorn, 1971 ), as a function of the leaf surface resistance factor a. When a = 1, the leaf surface resistances are equal to those calculated with the base inputs (Table 1 ); when a--16, the leaf surface resistances are 16 times those calculated with the base inputs. (b) Zero plane displacements determined by centroid method, plotted versus stability, which is a function of the leaf surface resistances. (c) Zero plane displacement estimated using similarity equations, plotted versus stability, which is changed by changes in the leaf surface resistance. (d) Roughness length estimated using similarity equations.
U (1987). A value of a = 1 on the x-axis indicates model output calculated with the base inputs. Stomatal and soil resistance changes did not greatly affect the zero plane displacement (Fig. 1 (a)). The figure shows that although the leaf surface resistances changed by several orders of magnitude, the zero plane displacement of momentum, estimated from the model output using eq. 1, did not change much from a value between 0.6 and 0.7he. This is consistent with previous studies. To facilitate comparison with previous reports of zero plane displacements shifts associated with stability changes (Raupach, 1979), we have plotted the zero plane displacement against stability in Figs. 1 (b) and 1 (c). It should be
266 remembered, however, that the leaf and soil resistance changes are causing the changes in the zero plane displacements and the stability, simultaneously, and that the stability itself is not causing the zero plane displacement changes. Also, for the sake of brevity, in this section only a wet soil case is discussed in detail; similar results were found for the drier soil cases. The zero plane displacement for water vapor, estimated with eq. 1, was not as independent of leaf surface resistance as that for momentum (Fig. l ( a ) ) . This was not unexpected because the relative contribution of water vapor sources changed as the stomates were opened or closed. The zero plane displacement (dvc) decreased with increasing leaf surface resistance {Fig. 1 (a)), because closing of the stomates allowed the soil to become the major source of water vapor. An apparent dependency on stability is created because both the source-sink distribution and stability are changed when the leaf surface resistance is changed (Fig. 1 (b)). The sensible heat flux zero plane displacement, estimated by eq. 1, was strongly affected by the changing leaf surface resistance (Fig. 1 {a) ). This also results in an apparent dependence on stability (Fig. 1 (b)). The displacement height was sometimes above the canopy height. This was because the sensible heat flux was negative at the soil surface and lower part of the canopy, and sometimes positive at the canopy top. The denominator of eq. 1, when applied to sensible heat flux, could have a value close to zero, resulting in an unrealistically large value for the displacement height. The similarity-derived zero plane displacement for momentum was generally independent of leaf surface resistance changes and therefore stability (Fig. 1 (c) ). The values were consistent with previous results of 0.6-0.7. For the sake of brevity, we have not shown a few exceptional cases (high soil resistance) when the displacement height was unrealistic. The similarity-derived zero plane displacement for humidity was generally equal to that for momentum (Fig. 1 (c)) when the soil was wet. However, the similarity-derived zero plane displacement for heat was far from constant, departing from 0.6 for most stabilities, because of the leaf surface resistance changes. Although the canopy model has been validated for real field conditions, there is the possibility that under the conditions tested here, the profiles of the mean variables were not accurate enough to use for determination of displacement heights. In addition, within the "roughness sublayer" the usual forms of the fluxgradient relationships, based on similarity, may not hold (Raupach, 1979; Raupach et al., 1980). If this were the problem, one would expect similar problems appearing for all fluxes, and not mainly for the heat flux, however. Raupach's conclusions (Raupach et al., 1980) are based on the identification of horizontal inhomogeneity and wake effects, which are not directly modeled above the canopy in our numerical simulation. A related explanation is that similarity analysis includes only the zero plane displacement to account for the complicated heat flux source distribution within a canopy and non-local turbulent
267
processes which transport heat from these sources. As shown here, the zero plane displacement for heat flux cannot be considered a variable independent from the physiologically induced source-sink distribution, as is assumed in similarity analysis. The momentum roughness length ranged from 0.11 to 0.16h¢, which is well within typical values found in the literature (Fig. 1 (d)). However, the roughness length for heat is not reasonable, being less than zero under stable conditions and over 0.3h¢ for unstable conditions. This is caused by the calculated displacement height exceeding the canopy height and the reference height.
Analysis of bulk canopy resistance Figure 2 depicts the dependence of bulk canopy resistance rc on the individual surface resistances of the leaves. On the x-axis, the leaf surface resistances are divided by their base values, in the same way as in Fig. 1 (a). On the y-axis, the resulting bulk canopy resistance re, calculated using eq. 2, is divided by the bulk canopy resistance evaluated at the base value of the leaf surface resistances. A value of 1 on both axes occurs when the simulation is run with the base radiation input (Table 1 ) and the leaf physiology model of Meyers and Paw U (1987). The figure shows that when the soil resistance is high (15000 s m - l ) the
U
0.1
',r
0.1
.
.
.
.
. ,,,i
. . . . . . . .
w
1 10 G Fig. 2. Relative changes in the derived canopy resistance r~ (shown by the ratio of r~/r~.b, where r~,~ is the canopy resistance estimated at base i n p u t conditions) as a function of the leaf surface resistance factor a, for differing values of the soil surface resistance (solid line represents a soil resistance r~--15000 s m - ' , long dashed line, r~l = 100 s m-~, a n d short dashed line, r~l--0.06 s m ~). W h e n a-- 1, the leaf surface resistances are equal to those calculated with the base inputs (Table 1 ); when a = 16, the leaf surface resistances are 16 times those calculated with the base inputs.
268
00oj
o
.......
0
0
/
300q
o°/°
1
o ,~,~" ~ o
:0
Ioo
.o ? ~ . °
oo
. . . . . .
200 300 400 LE ( S i m u l a t e d )
5oo
600
Fig. 3. Relationship between the latent energy flux densities estimated using the Penman-Monteith equation and using higher-order closure, when the canopy resistance is calculated using a weighted stomatal resistance (eq. 3). Circlesrepresent canopy resistance without the ground (r,av), and the diamonds with the ground included (r, av~). curve is close to a 1:1 line (although not equal to it). The shape of the curve is caused by the canopy resistance changing less than the leaf surface resistance. The leaves and soil tend to humidify the local canopy environment, as the transpiration increases, which then decreases the leaf vapor pressure deficit. The above process is coupled with the complex turbulence within the plant canopy. W h e n the soil resistance is lower (i.e. the soil more "wet"), the curve departs greatly from a 1:1 relationship. This demonstrates that soil evaporation can become an important component of evapotranspiration, reducing the control by the plant canopy physiology. These results agree qualitatively with the experimental data of Brown and Covey (1966) and D e n m e a d (1973). W h e n the canopy resistances are estimated with eq. 3, and the aerodynamic resistance estimated from the simulation outputs, large errors in estimated latent energy flux density can occur (Fig. 3 ). The errors appear even when the ground surface contribution is small, and are related to the method of estimating rear. The parallel resistance weighted by leaf area index is not an accurate measure of transfer in the canopy, where canopy source-sinks and turbulent transport interact. CONCLUSIONS
Numerical "experiments" with our higher-order plant canopy closure model reveal that while the zero plane displacement for m o m e n t u m flux is relatively
269
constant with changes in stomatal conductance, the displacement heights for water vapor flux and especially heat flux are not. The changing source-sink relationships of a plant canopy, controlled by plant physiology and soil parameters, create a complex picture of the effective source height within the canopy. The reversal of the sensible heat flux within the canopy is part of the problem. The roughness length for heat also exhibits problems similar to that for the zero plane displacement. These problems result in difficulties in using similarity theory a n d / o r bulk aerodynamic resistance descriptions of canopies. The bulk canopy resistance description of canopies is also shown to be problematic. The problems are least when the soil is dry. Even so, when the canopy resistance is estimated using leaf-area weighted parallel resistances, serious errors can arise. ACKNOWLEDGEMENTS
This research was supported in part by National Aeronautics and Space Administration grant number NAG5-937 and National Science Foundation grant number ATM-8703832. A portion of this work was conducted as a contribution to the National Acid Precipitation Assessment Program. Discussions with Dr. Dennis Baldocchi aided our research. Revision of the manuscript was greatly aided by comments from several reviewers, including Drs. Don Aylor and Roger Shaw. Word processing help was given by Ms. Cherie Felsch. APPENDIX
The similarity equations of Brutsaert (1982) were inverted to solve for the zero plane displacement and the roughness length. Under unstable conditions, the similarity equation for m o m e n t u m is
u*
Oz-L
\ L ]
(A1)
where k is von Karman's constant, u* is the friction velocity, (9U/c)z is the shear of the mean wind profile, L is the Monin-Obukhov length and ds is the zero plane displacement for momentum. Under stable conditions, the equation for ds is
ds = z -
-~-~;
(A2)
with similar equations for heat flux and water vapor (see Brutsaert, 1982). The roughness length may be determined from these equations after the zero plane displacement has been estimated
270
zo= (z-ds)exp[-(~+ ~")l
(A3)
where ~ is the integrated stability function (Brutsaert, 1982).
REFERENCES Baldocchi, D.D., 1982. Mass and energy exchanges of soybeans: Microclimate-plant architectural interactions. Ph.D. Dissertation, CAMaC Progress Report 82-3, University of Nebraska, Lincoln, NE, 265 pp. Brown, K.W. and Covey, W., 1966. The energy-budget evaluation of the micrometeorological transfer processes within a cornfield. Agric. Meteorol., 3: 73-96. Brutsaert, W., 1982. Evaporation into the Atmosphere. Reidel, Boston, 299 pp. Choudhury, B.J., Reginato, R.J. and Idso, S.B., 1986. An analysis of infrared temperature observations over wheat and calculation of latent heat flux. Agric. For. Meteorol., 37: 75-88. Cowan, I.R., 1968. Mass, heat and momentum exchange between stands of plants and their atmospheric environment. Q. J. R. Meteorol. Soc., 94: 523-544. Denmead, O.T., 1973. Relative significance of soil and plant evaporation in estimating evapotranspiration. In: Plant Response to Climatic Factors. UNESCO, Paris, pp. 505-511. Garratt, J.R., 1978. Flux profile relations above tall vegetation. Q. J. R. Meteorol. Soc., 104: 199211. Hicks, B.B., Hess, G.D. and Wesely, M.L., 1979. Analysis of flux-profile relationships above tall vegetation - an alternative view. Q. J.R. Meteorol. Soc., 105: 1074-1077. Jackson, P.S., 1981. On the displacement height in the logarithmic velocity profile. J. Fluid Mech., 111: 15-25. Meyers, T.P., 1985. A simulation of the canopy microenvironment using higher order closure principles. Ph.D. Dissertation, Purdue University, West Lafayette, IN, 153 pp. Meyers, T.P. and Paw U, K.T., 1987. Modelling the plant canopy micrometeorology with higherorder closure principles. Agric. For. Meteorol., 41: 143-163. Meyers, T.P., Huebert, B.J. and Hicks, B.B., in press. HNO3 deposition to a deciduous forest. Wellus. Monteith, J.L., 1977. Resistance of a partially wet canopy: Whose equation fails? Boundary-Layer Meteorol., 379-383. Norman, J.M., 1979. Modeling the complete crop canopy. In: B. Barfield and J. Gerber (Editors), Modification of the Aerial Environment of Crops. Am. Soc. Agric. Eng., St. Joseph, MI, pp. 249-277. Paw U, K.T., 1987. Mathematical analysis of the operative temperature and energy budget. J. Therm. Biol., 12: 227-233. Paw U, K.T. and Gao, W., 1988. Applications of solutions to non-linear energy budget equations. Agric. For. Meteorol., 43: 121-145. Paw U, K.T., Shaw, R.H. and Meyers, T.P., 1985. Evapotranspiration as modeled by a higher order closure scheme. In: Advances in Evapotranspiration. Am. Soc. Agric. Eng., St. Joseph, MI, pp. 43-50. Raupach, M.R., 1979. Anomalies in flux-gradient relationships over forests. Boundary-Layer Meteorol., 16: 467-486. Raupach, M.R., Stewart, J.B. and Thom, A.S., 1979. Comments on the paper 'Analysis of fluxprofile relationships above tall vegetation - an alternative view'. Q. J. R. Meteorol. Soc., 105: 1077-1078.
271 Raupach, M.R., Thorn, A.S. and Edwards, I., 1980. A wind-tunnel study of turbulent flow close to regularly arrayed rough surfaces. Boundary-Layer Meteorol., 18: 373-397. Rosenberg, N.J., Blad, B.L. and Verma, S.B., 1983. Microclimate: The Biological Environment. Wiley, New York, 495 pp. Sellers, P.J. and Dorman, J.L., 1987. Testing the simple biosphere model (SiB) using point micrometeorological data. J. Climate Appl. Meteorol., 26: 622-651. Shaw, R.H. and Pereira, A.R., 1982. Aerodynamic roughness of a plant canopy: a numerical experiment. Agric. Meteorol., 26: 51-65. Shaw, R.H. and Seginer, I., 1987. Calculation of velocity skewness in real and artificial plant canopies. Boundary-Layer Meteorol., 39: 315-332. Shuttleworth, W.J., 1975. The concept of intrinsic surface resistance: Energy budgets at a partially wet surface. Boundary-Layer Meteorol., 8: 81-99. Shuttleworth, W.J., 1976. Experimental evidence for the failure of the Penman-Monteith equation in partially wet conditions. Boundary-Layer Meteorol., 10: 91-94. Thorn, A.S., 1971. Momentum absorption by vegetation. Q. J. R. Meteorol. Soc., 97: 414-428.