UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A study of soil temperature and seedling survival in a forest clearcut Stathers, Robert John 1983

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
831-UBC_1983_A6_7 S83.pdf [ 5.43MB ]
Metadata
JSON: 831-1.0095699.json
JSON-LD: 831-1.0095699-ld.json
RDF/XML (Pretty): 831-1.0095699-rdf.xml
RDF/JSON: 831-1.0095699-rdf.json
Turtle: 831-1.0095699-turtle.txt
N-Triples: 831-1.0095699-rdf-ntriples.txt
Original Record: 831-1.0095699-source.json
Full Text
831-1.0095699-fulltext.txt
Citation
831-1.0095699.ris

Full Text

A STUDY OF SOIL TEMPERATURE AND SEEDLING SURVIVAL IN A FOREST CLEARCUT By ROBERT JOHN STATHERS B.Sc.(Agr.)» The University of British Columbia, 1981 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE ln THE FACULTY OF GRADUATE STUDIES Department of Soi l Science We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA / August 1983 (g) Robert John Stathers, 1 9 8 3 I n p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f o r an advanced degree a t the U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by t h e head o f my department o r by h i s o r h e r r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l n o t be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . Department o f So'A Sg-^ ice. The U n i v e r s i t y o f B r i t i s h C o l u m b i a 2075 Wesbrook P l a c e Vancouver, Canada V6T 1W5 Date OCT . l l , R 63 E-6 (2/79) ABSTRACT This thesis reports the results of a study of s o i l temperature and seedling survival on a steep south facing, high elevation forest clearcut on Vancouver Island, B.C., having a history of regeneration failures. During the f i r s t year of the study (1981) there were harsh weather condi-tions during the summer, while during the second year (1982) summer weather conditions were much milder. Lethally high surface s o i l temperatures were found to be the major cause of seedling mortality on this site. Of 2160 Douglas-fir, western hemlock and Pacific silver f i r seedlings planted in April 1981, 87, 33 and 17% respectively had survived by the end of the second growing season. Shade cards significantly increased survival for each of these species. Of an additional 960 Douglas-fir, western hemlock, and Pacific silver f i r seedlings planted in April 1982, 92, 74, and 88% respectively had survived by October 1982. Shade cards, trickle i r r i g a -tion, and shade cards and trickle irrigation combined only significantly increased survival in western hemlock in the seedlings planted in 1982 because of the high overall survival rate resulting from the mild summer weather conditions. A physically based s o i l temperature model was developed to (i) quanti-fy the effects of the factors which influence the s o i l thermal regime and ( i i ) to model the time course of s o i l profile temperatures from easily measured or estimated meteorological data and site specific characteris-t i c s . The model, based upon a f i n i t e difference solution of the heat flow equation, energy and radiation balance, and aerodynamic theory, calculates surface and subsurface temperatures i n soils having thermal properties that I f f vary with depth. A slightly modified version of the atmospheric sta b i l i t y correction method of Paulson (1970) was required to accurately predict surface temperatures under unstable atmospheric conditions. The model was tested with data presented in two other energy balance studies reported in the literature and found to estimate s o i l temperatures very accurately. In addition, good agreement was found between modelled and measured s o i l profile temperatures for 6 and 16 day periods in the forest clearcut mentioned above. Energy balance theory was used to show that maximum surface s o i l temp-eratures in a forest clearcut are not likely to exceed 85°C on a clear hot day. Wind speed, aerodynamic roughness of the surface, and s o i l thermal properties were also shown to strongly influence s o i l temperatures. TABLE OF CONTENTS Page ABSTRACT JL TABLE OF CONTENTS J& LIST OF TABLES J Z D LIST OF FIGURES IZfflt LIST OF SYMBOLS SJL ACKNOWLEDGEMENTS SSS INTRODUCTION ± CHAPTER 1 . EFFECT OF MICROCLIMATE MODIFICATION ON THE SURVIVAL OF 3 SPECIES OF CONIFER SEEDLINGS 1.1 Introduction 4 1.2 Experimental Site $ 1.3 Experimental Design and Procedure £ 1.4 Results and Discussion A. Survival of Seedlings Planted in 1981 JO B. Survival of Seedlings Planted in 1982 13 C. The Effect of Weather on Seedling Survival 73 D. The Effects of Shade and Irrigation on Soil Temperature 2.1 1.5 Conclusions 25 CHAPTER 2 . FACTORS AFFECTING SURFACE SOIL TEMPERATURES IN A FOREST CLEARCUT 2.1 Introduction 29 2.2 Theory A. The Energy and Radiation Balances of a Forest Clearcut Z9 is: Page B. The Soil Heat Flux 31 C. Aerodynamic Equations Describing the Sensible Heat Flux, Latent Heat Flux and Momentum Flux in the Turbulent Boundary Layer 32-2.3 Results and Discussion A. Surface Energy Balance and Clearcut Microclimate ... 4/ a) The Radiation Balance 4/ b) The Soil Thermal Regime ^ £ c) The Soil Water Regime 4.3 d) Boundary Layer Resistance to Heat, Water Vapour and Momentum Transfer 45 B. Relationship Between Surface Energy Balance Components and Surface Temperature 45 C. Comparison of Atmospheric Stability Correction Methods 4g D. Examination of Micrometeorological Factors Affecting Surface Temperature 53 2.4 Conclusions §"J CHAPTER 3. MODELLING SOIL TEMPERATURE IN A FOREST CLEARCUT 3.1 Introduction 50 3.2 Theory 67 3.3 Procedure A. Assumptions Used in the Model 64 B. The Solution Method 65 3.4 Results and Discussion A. Validation of the Model 57 B. Sensitivity Analysis of the Model "77 3.5 Conclusions a-7 Page CONCLUDING REMARKS <\<\ REFERENCES J 01 APPENDIX A. Derivation of the 1-Dimensional Finite Difference Solution of the Soil Heat Conduction Equation jQfa APPENDIX B. Modelling Clear Sky Solar Irradiance on Variable Slopes and Aspects //£. APPENDIX C. Soil Thermal Properties JJQ APPENDIX D. A Listing of the Computer Program Used to Calculate Soil Temperatures and Surface Energy Fluxes 723 VTT LIST OF TABLES CHAPTER 1. TABLE 1.1 TABLE 1.2 TABLE 1.3 TABLE 1.4 Monthly percentage survival of seedlings which were planted in April 1981 Page 77 Percent survival in October of 1981 and 1982 of seedlings planted in April 1981. Numbers followed by the same letter are not significantly different (t = 0.05) according to Duncan's New Multiple Range Test Monthly percentage survival of seedlings which were planted in April 1982 , Percent survival in October 1982 of seedlings planted in April 1981. Numbers followed by the same letter are not significantly different (t = 0.05) according to Duncan's New Multiple Range Test 11 H 15 CHAPTER 2. TABLE 2.1 TABLE 2.2 Thermal Properties of Natural Materials .. Aerodynamic Properties of Natural Surfaces 44 46 3zm CHAPTER 1. LIST OF FIGURES Page FIGURE 1.1 FIGURE 1.2 FIGURE 1.3 FIGURE 1.4 FIGURE 1.5 FIGURE 1.6 Daily values of precipitation, solar irradiance, average air temperature, and average relative humidity during the 1981 growing season J7 Profiles of gravimetric s o i l water content at selected times during the summer of 1981 Daily values of precipitation, solar irradiance, average air temperature, and average relative humidity during the 1982 growing season 20 The diurnal course of s o i l temperature at several depths, surface s o i l temperature (measured with a fine (0.075 mm diameter) wire copper con-stantan thermo-couple) and air temperature at the 1 m height on August 7, 1982 23 The diurnal course of surface s o i l temperature adjacent to shaded, irrigated, shaded and irrigated, and control treatments on August 7, 1982. The course of air temperature at the 1 m height i s also shown The diurnal course of s o i l temperature in the seedling root zone adjacent to shaded, irrigated, shaded and irrigated, and control treatments on August 7, 1982. Numbers above the profile refer to the time of day (P.S.T.) Z4 26 CHAPTER 2. FIGURE 2.1 Integrated diabatic profile correction factors for heat and momentum versus 6 or R^  (indices of atmospheric stability) for unstable, neutral, and stable atmospheric conditions 40 FIGURE 2.2 Surface s o i l temperature versus solar irradiance calculated for a range of values of H + LE + G D from energy balance theory (equation 2.7) assuming a = 0.18, e s = 9.7, c a c = 0.85, and T a = 25°C 41 I X FIGURE 2.3 Surface s o i l temperature versus wind speed calculated using energy balance and aerodynamic theory and the sta b i l i t y correction methods of Paulson (1970), Thorn et al.(1975), Holbo (1971), Fuchs et al,(1968), and Hammel et al.(1982). Values were calculated assuming Rs = 1000 W m-2, G + LE = 150 W m-2, T a = 25°C, ct = 0.18, e s = 0.97, e a = 0.85 and z Q = 0.01 m. The methods of Fuchs et a l and Hammel et al„give identical results (see text).. FIGURE 2.4 Calculated versus measured surface s o i l temperatures on August 9, 1981. Calculated values were obtained from measured Rs, T a, u, GQ and LE/Rg = 0.05 using the stability correction methods of Paulson (1970), Thorn et al,(1975), Fuchs et al.(1968) and Hammel et al e(1982). Surface s o i l temperatures were measured with an infra-red thermometer and a fine wire thermocouple FIGURE 2.5 Surface s o i l temperature versus wind speed at at the 1 m height for various solar irradiances. Values were calculated assuming G Q + LE = 150 W m-2, T a = 25°C, z Q = 0.01 m, a = 0.18, e s = 0.97, and e a = 0.85. FIGURE 2.6 Surface s o i l temperature versus wind speed at air temperatures of 5, 25, and 45°C. Values were calculated assuming Rs = 1000 W m" 2 , G + LE = 150 W m-2, z 0 = 0.01 m, a = 0.18, e s = 0.97, and e a = 0.85. FIGURE 2.7 Surface s o i l temperatures versus wind speed at aerodynamic roughness lengths of 0.001, 0.01, and 0.1 m. Values were calculated assuming Rs = 1000 W m-2, G Q + LE = 150 W m-2, T a = 25°C, a = 0.18, e s = 0.97, and e a = 0.85 CHAPTER 3. FIGURE 3.1 Flow chart showing the solution method used to calculate surface and subsurface s o i l temperatures. FIGURE 3.2 Comparison of modelled and measured (Stearns 1967) temperatures in a dry Peruvian desert s o i l having uniform thermal properties with depth on July 11, 1964 FIGURE 3.3 Comparison o f m o d e l l e d and measured (Novak 1981) t e m p e r a t u r e s and s u r f a c e energy f l u x e s on May 18, 1979 i n a wet s i l t loam s o i l a t A g a s s i z , B.C FIGURE 3.4 Comparison o f m o d e l l e d and measured (Novak 1981) t e m p e r a t u r e s and s u r f a c e energy f l u x e s on J u l y 2 1 , 1979 i n a s i l t loam s o i l a t A g a s s i z , B.C. w h i c h had d r i e d near the s u r f a c e FIGURE 3.5 Comparison o f m o d e l l e d (0) and measured (+) s o i l t e m p e r a t u r e s a t 4 dep t h s i n t h e Cameron V a l l e y f o r e s t c l e a r - c u t e x p e r i m e n t a l s i t e d u r i n g the 6 h o t t e s t days i n 1981 FIGURE 3.6 Comparison o f m o d e l l e d (0) and measured (+) s o i l t e m p e r a t u r e s a t 3 depths i n t h e Cameron V a l l e y f o r e s t c l e a r - c u t e x p e r i m e n t a l s i t e d u r i n g 16 days of v a r i a b l e weather c o n d i t i o n s i n J u l y 1981 FIGURE 3.7a The e f f e c t o f t h e aerodynamic roughness ( z Q ) on m o d e l l e d s u r f a c e s o i l t e m p e r a t u r e . V a l u e s o f z Q used were 5 mm ( x ) , 10 mm (+), and 50 mm (0) w h i l e a l l o t h e r v a r i a b l e s were h e l d c o n s t a n t . FIGURE 3.7b As f o r F i g u r e 3.7a but f o r t e m p e r a t u r e s a t t h e 0.1 m d e p t h FIGURE 3.7c As f o r F i g u r e 3.7a but f o r t e m p e r a t u r e s a t t h e 0.3 m d e p t h FIGURE 3.7d As f o r F i g u r e 3.7a but f o r t e m p e r a t u r e s a t t h e 0.5 m d e p t h FIGURE 3.8a The e f f e c t of the l a t e n t h eat f l u x d e n s i t y ( L E ) on m o d e l l e d s u r f a c e s o i l t e m p e r a t u r e . V a l u e s o f LE used were 0.0 ( + ) , 0.05 R s ( x ) , and 0.2 R s (0) w h i l e a l l o t h e r v a r i a b l e s were h e l d c o n s t a n t . FIGURE 3.8b As f o r F i g u r e 3.8a but f o r t e m p e r a t u r e s a t t h e 0.1 m d e p t h FIGURE 3.8c As f o r F i g u r e 3.8a but f o r t e m p e r a t u r e s a t t h e 0.3 m d e p t h FIGURE 3.8d As f o r F i g u r e 3.8a but f o r t e m p e r a t u r e s a t t h e 0.5 m d e p t h Page FIGURE 3.9a The effects of s o i l thermal properties on model-ed surface s o i l temperature. Values of k s and C s in the 0-0.02, 0.02-0.05 and 0.05-0.8 m s o i l layers were 1.0, 1.2 and 1.4 MJ m-3 °C _ 1 and 0.19 0.29 and 0.67 W m_1 °C _ 1 respectively (X), values of k s are the same as those given above but values of C s are twice those given above (+) and values of k s and C s are twice those given above (0). A l l other variables were held constant. Qfl FIGURE 3.9b As for Figure 3.9a but for temperatures at the 0.1 m depth tyQ FIGURE 3.9c As for Figure 3.9a but for temperatures at the 0.3 m depth 9/ FIGURE 3.9d As for Figure 3.9a but for temperatures at the 0.5 m depth 92 FIGURE 3.10a The effect of the temperature at the lower sur-face of the modelled s o i l slab (0.8 m depth) [lower boundary condition] on modelled surface s o i l temperature. Values of 10°C (+), 12°C (x), and 14°C (0) were used while a l l other variables were held constant 9 3 FIGURE 3.10b As for Figure 3.10a but for temperatures at the 0.1 m depth 9^  FIGURE 3.10c As for Figure 3.10a but for temperatures at the 0.3 m depth 9 5 FIGURE 3.10d As for Figure 3.10a but for temperatures at the 0.5 m depth 9 6 ixn L I S T OF SYMBOLS C volumetric heat capacity (J m °C - ) C cloud cover dimensionless C s volumetric heat capacity of s o i l (J m °C ) F atmospheric s t a b i l i t y c o rrection factor dimensionless G Q s o i l heat f l u x density (W m ) H sensible heat f l u x density (W m ) H' sensible heat f l u x density under neutral (W m ) conditions Kg eddy d i f f u s i v i t y for sensible heat (m s ) KM eddy d i f f u s i v i t y f o r momentum (eddy v i s c o s i t y ) (m s ) Ky eddy d i f f u s i v i t y for water vapour (m s ) L Monin-Obukhov s t a b i l i t y length ( p c p T au* 3/(kgH)) (m) LE latent heat f l u x density (W m ) LE' latent heat f l u x density under neutral (W m ) conditions Ri gradient Richardson number dimensionless — 2 R^ longwave irradiance (W m ) _ o R n net r a d i a t i o n f l u x density (W m ) R s solar irradiance (W m ) T temperature (°C) T a a i r temperature (°C, K) T a mean a i r temperature between z and z Q (K) T s s o i l surface temperature (°C, K) -XZ7J c r a t i o of LE to R s dimensionless C p a s p e c i f i c heat of a i r (kj kg -* ° C - 1 ) C p S s p e c i f i c heat of s o i l (kJ k g - 1 ° C - 1 ) e vapour pressure of water i n a i r (kPa) e Q vapour pressure of water i n a i r at s o i l surfce (kPa) e*(T s) saturation vapour pressure of water (kPa) g g r a v i t a t i o n a l a c c e l e r a t i o n (m s ) k von Karman's constant dimensionless k thermal conductivity (W m - 1 "C""1) k s thermal conductivity of s o i l (W m - 1 ° C - 1 ) u wind v e l o c i t y (m s -*) u* f r i c t i o n v e l o c i t y (m s - 1 ) , z depth below or height above s o i l surface (both p o s i t i v e ) (m) z Q aerodynamic roughness length (m) a shortwave r e f l e c t i v i t y (albedo) dimensionless Y psychrometric constant (Pa ° C _ 1 ) e viscous d i s s i p a t i o n of k i n e t i c energy dimensionless e a atmospheric emissivity dimensionless e a c c l e a r sky emissivity dimensionless G S s o i l surface emissivity dimensionless £ index of atmospheric s t a b i l i t y (z/L) dimensionless q p a density of a i r (kg m ) Ps a T • H •I'M <t>v density of s o i l Stephan-Boltzmann constant momentum flux density diabatic influence function for sensible heat transfer diabatic influence function for momentum transfer diabatic influence function for water vapour transfer (kg m"3) (W m - 2 lCk) (N m-2) dimensionless dimensionless dimensionless ACKNOWLEDGEMENTS I wish to express my sincere gratitude and appreciation to Dr. T.A. Black for his invaluable guidance and supervision during my M.Sc. program. I also wish to thank MacMillan Bloedel Ltd., Woodlands Services D i v i s i o n for funding t h i s research. My appreciation goes to the other members of my committee: Dr. M.D. Novak (Assistant Professor, Department of S o i l Science), Dr. T.M. B a l l a r d (Professor, Department of S o i l Science), and Mr. B.G. Dunsworth (Research S i l v i c u l t u r i s t , MacMillan Bloedel L t d . ) . I extend my thanks to the biometeorology group at U.B.C.: Dr. D.L. Spittlehouse, Mr. F. K e l l i h e r , Mr. D.T. Pr i c e , Mr. D.G. G i l e s , and i n p a r t i c u l a r , Mr. N.J. Livingston, for t h e i r valuable discussions and a s s i s t -ance . I also wish to thank the Department of S o i l Science for t h e i r accept-ance of me into the M.Sc. program, and Ms. D. Shunamon for her excellent word processing of th i s t h e s i s . F i n a l l y , my wholehearted and sincere thanks goes to my family, wife and daughter for the i r encouragement and precious f r i e n d s h i p . INTRODUCTION Soil temperature is one of the major factors affecting conifer seed-ling regeneration in British Columbia. In the north central interior region, low s o i l temperatures inhibit growth (Dobbs and McMinn 1977). In the southwest coastal region, lethally high surface temperatures contribute to mortality (Arnott 1975). The problems associated with seedling regeneration on sites prone to adverse microclimatic conditions have been recognized for a long time (Baker 1929; Maguire 1955) and have prompted considerable research. Various techniques such as shading (Hobbs et a l . 1980; Minore 1970), mulch-ing (Hermann 1964), irrigating (Livingston and Stathers 1982), and modify-ing the s o i l thermal properties (Hermann and Chilcote 1965) have been used to ameliorate harsh environmental conditions. The objectives of this study were to investigate practical methods of improving seedling survival on a high-elevation south facing clearcut on Vancouver Island, B.C. The experimental site, located in the Cameron Valley near Port Alberni, typically experiences harsh environmental condi-tions during the summer and was replanted with different species at least 3 times before an adequate stocking level was obtained. Over the 2 year course of this study, 3 species of conifer seedlings were tested under a number of microclimate modification treatments. These tests prompted a more detailed analysis of the factors influencing seedling microclimate and as a result a physically based computer simulation model was developed. The results of these endeavours form the basis of the three chapters in this thesis. In Chapter 1, the results of seedling survival studies initiated in 1981 and 1982 are reported. The effects of weather and microclimate modi-fication on survival are also discussed. In Chapter 2, the environmental factors influencing surface s o i l temperature are outlined. Relationships between the energy and radiation balances of the surface and surface s o i l temperature are developed. Aero-dynamic methods of calculating the sensible and latent heat fluxes from the surface are also discussed. In Chapter 3, a s o i l temperature model capable of predicting surface temperatures is developed and tested. The model is used to (i) simulate s o i l temperatures over extended time periods using easily measured climate data and ( i i ) simulate the effects of changing certain environmental vari-ables on s o i l temperatures. -3-CHAPTER 1 EFFECT OF MICROCLIMATE MODIFICATION ON THE SURVIVAL OF THREE SPECIES OF CONIFER SEEDLINGS--4-EFFECT OF MICROCLIMATE MODIFICATION ON THE SURVIVAL OF THREE SPECIES OF CONIFER SEEDLINGS 1.1 INTRODUCTION Plantation establishment can be very d i f f i c u l t on the high elevation south facing clearcuts on Vancouver Island. The accumulated winter snow-pack melts late in the spring, leaving very l i t t l e time for newly planted seedlings to acclimatize to intense solar irradiance during the early summer. The shallow, dense, coarse textured soils commonly found at these sites have a low water holding capacity and the seedling root zone dries rapidly after the snow melts. When drought conditions persist for more than a few weeks during the summer, c r i t i c a l l y high surface s o i l tempera-tures and low s o i l water potentials can occur. These conditions are thought to be the major cause of seedling mortality at high elevations i n this region (Arnott 1975). Many operational techniques have been applied to improve seedling survival where moisture and temperature stress occur. The correct tree species, provenance, and stock type selection have reduced seedling mortal-ity; however, under the extreme environmental conditions of summer drought, additional protective measures are necessary. Modifying the seedling microclimate may be the appropriate, cost-effective solution to the regen-eration problems on these sites. The objectives of this study were to assess the effects of moisture and temperature stress reduction on the survival of three seedling species on a high elevation south facing clearcut on Vancouver Island, B.C. The methods used to reduce stress included natural shade, a r t i f i c i a l shade, -5-irrigation, and a r t i f i c i a l shade and irrigation combined. A r t i f i c i a l shade from shade cards has been shown to increase survival on stressful sites in Oregon (Hobbs et al,1980; Minore 1970) and thus was of interest as a pre-scribed operational treatment. Natural self shading by inclining the stem to the south has also been effective in seedling nursery beds (Korstian and Fetherolf 1921) and was considered an inexpensive, easily applied opera-tional procedure. Irrigation has long been recognized as the most amelior-ative procedure for stress reduction (Baker 1929; Maguire 1955); however, i t is usually not a cost effective alternative in remote areas. Irrigation treatments were applied in this study to determine seedling performance under conditions of minimal water stress, and also to determine whether water or thermal stress caused seedling mortality. 1.2 EXPERIMENTAL SITE The experimental site was located at an elevation of 1150 m on a 30° south facing slope in the Cameron Valley near Mt. Arrowsmith, Vancouver Island, British Columbia (49° 13'N, 124° 36'W). The area is in the Mari-time Coastal Western Hemlock biogeoclimatic subzone (Klinka et al,1970). The site was very representative of the surrounding area. The s o i l at this site is an orthic Humo-ferric Podzol. It has devel-oped over compact glacial t i l l and contains a large proportion of coarse fragments («60% > 2mm in diameter by dry weight). It has a sandy loam _ o texture. The bulk density is 1.1-1.3 Mg m near the surface and 1.6-1.8 — 3 Mg m below the 0.5 m depth. These characteristics give i t a very low available water storage capacity. A thin (<0.05 m) discontinuous layer of organic matter is present on the surface. ~6~ The site was harvested in 1974. The stand was comprised predominantly of Douglas-fir (Pseudotsuga menziesii (Mirb.) Franko), western hemlock (Tsuga heterophylla (Raf.) Sarg.) and yellow cedar (Chamaecyparis nootka-tensis (D. Don) Spach). A sparse canopy of fireweed (Epilobium angusti-folium) and salal (Gaultheria shallon) now covers the site. Douglas-fir seedlings were planted at least 3 times prior to this study before an adequate stocking level was achieved. These seedlings were less than 1 m t a l l when this study was initiated. 1.3 EXPERIMENTAL DESIGN AND PROCEDURE In late April 1981, 2160 seedlings were planted on the experimental site. High-elevation provenances of Douglas-fir, western hemlock, and Pacific silver f i r were selected for the t r i a l . The same number of seed-lings for each species were planted. 1-0 container grown styroplugs were used exclusively because of the operational advantages of producing and planting this stock type. Seedlings were maintained in a dormant state i n cold storage for 4 months prior to planting. Nine combinations of tree species and microclimate modification treat-ments were tested. An a r t i f i c i a l shade treatment, a natural self-shade treatment, and an unshaded control treatment were tested for each of the 3 species. The a r t i f i c i a l shade treatment was imposed by placing a 0.2 x 0.3 m waxed cardboard shade card (manufactured by International Forest Supplies, Eugene, Oregon) on the southwest side, 0.05 m from the stem, and sloping approximately 10° back over the seedling. This orientation was necessary to effectively shade the root collar region because of the high solar zenith angle during the middle of the day. The natural self-shade treatment was imposed by inclining the seedling stem to the southwest so -7-that the foliage shaded the root collar region at midday. Unshaded seed-lings were planted vertically and the root collar region was exposed to f u l l sunlight. Microsite selection and planting techniques were similar for a l l species and treatments. Seedlings were planted away from logging debris and in scarified mineral s o i l . Good hydraulic contact was established between the root system and the s o i l by carefully hand-compacting the disturbed s o i l near the roots. The seedlings were laid out on the site in a randomized complete block design. Each of the 9 species-treatment combinations was replicated in 24 blocks. In every block, a treatment consisted of 10 seedlings, planted 1.5 m apart in the downslope direction. Consequently, there were 240 seed-lings per treatment on the site. In 1982, 4 additional blocks of seedlings were planted on the same site. The same 3 species, provenance, stock type, and planting technique were used; however 4 treatments were tested. Treatments included shade cards to reduce heat stress as in 1981, trickle irrigation to reduce plant water stress, a combination of shade and irrigation to minimize water and heat stress, and an unshaded unwatered control. A trickle irrigation system (Drip-Eze Systems, El Cajon, California, U.S.A.) was used to deliver water through drip emitters to individual seed-lings (Livingston 1984). Water was applied when tensiometers on the site indicated that the s o i l water potential was less than -0.05 MPa. A r t i -f i c i a l shade was imposed using the same method described in the 1981 t r i a l . Analysis of variance was used to determine significant differences in survival between blocks and treatments. Duncan's New Multiple Range -8-Test (t = 0.05) was used to rank each treatment. Seedling survival was assessed monthly from May through October i n 1981 and 1982 in order to relate mortality to climatic conditions. A weather station was maintained on the site from May until mid-September in both years since summer conditions were suspected to be the most important factor influencing seedling survival. Solar irradiance, net radiation, air temperature (1 m height), relative humidity (1 m height), wind speed (1 m height), precipitation, s o i l temperature (0.5 m, 0.3 m, 0.05 m, 0.015 m depths, and the surface) and s o i l heat flux (0.03 m) were recorded on battery-powered Campbell Scientific CR21 data loggers (Campbell Scientific, Logan, Utah, U.S.A.). The data loggers were programmed to record hourly averages and store the data on audio cassette tapes. Data loggers were housed in plastic (Tupperware) containers and kept dry with s i l i c a gel desiccant. The plastic containers were kept in watertight metal military ammunition boxes which were locked and chained to an underground concrete pad. The boxes were located in a covered 0.6 m deep hole to protect the data loggers from the harsh climatic conditions at the surface and make them less vulnerable to vandalism. The pyranometer (Model Li-200SZ-05, Li-Cor Inc., Lincoln, Nebraska, U.S.A.) was calibrated with a laboratory standard Kipp and Zonen pyrano-meter at the beginning and end of each season. The net radiometer (Type S-l, Swissteco Pty. Ltd., Melbourne, Australia) calibration was determined by comparison with similar model laboratory standard net radiometer. Net radiation was measured during July and August of both years. Relative humidity was measured with a sulphonated polystyrene sensor (Model 201, Phys-Chemical Research Corp., New York, New York, U.S.A.). - 9 -It was found to have a ±3% accuracy over the range from about 10-97%. Rainfall was measured with a tipping bucket raingauge (Model RG2501, Sierra Misco, Inc., Berkeley, California, U.S.A.) which had a resolution of 1 mm. Wind speed was measured at the 1.0 m. height with a cup anemometer (Model 014A.Met-One, Sunneyvale, California, U.S.A.) which had a stalling speed of 0.5 m/s. Air and deep s o i l temperatures were measured with thermistor (Fenwall UUT-51J1) thermometers which were accurate to ±0.5°C between 0 and 45°C. Temperatures at and near the s o i l surface were measured with fine (0.075 mm diameter) wire copper constantan thermocouples which were refer-enced to a deep s o i l thermistor thermometer. Surface temperatures were also measured with an infrared thermometer (Instatherm Model 14-2200-4 close focus, Barnes Engineering Co., Stanford, Conn., U.S.A.). Soil heat flux at the 0.03 m. depth was measured with 3 heat flux plates connected in series (made by Dr. M.D. Novak, Dept. of Soil Science, U.B.C. following Fuchs and Tanner (1968)). Measurements of the s o i l water content and potential of the seedling root zone were collected on a weekly basis throughout the summer. A neu-tron probe (Model 503 hydroprobe, Campbell Pacific Nuclear, Pacheco, C a l i -fornia, U.S.A.) was used to monitor the s o i l water content over a range of depths from 0.15 to 0.75 m. Gravimetric sampling was used to determine seedling root zone water contents and to calibrate the neutron probe. Soil water potentials were measured on an HR33T dew point hygrometer with 40 so i l psychrometers (Model PT51-10, Wescor Inc., Logan, Utah, U.S.A.) plac-ed in the seedling root zone at the 0.15 m depth. The psychrometers were calibrated in the laboratory using a range of KC1 solutions. Care was tak-en to ensure good hydraulic contact between the s o i l and the horizontally -10-oriented psychrometers. Readings were taken in the morning before large s o i l temperature gradients occurred, using the dew point mode and correct-ing the cooling coefficient (IIV) for temperature variation. Seedling shoot water potentials were measured throughout the summer with a Scholander pressure bomb. Seedling leaf conductances were also measured in situ using a transient state chamber porometer developed for this study and described in Livingston et al,(1983). 1.4 RESULTS AND DISCUSSION A. Survival of Seedlings Planted In 1981 Survival of seedlings planted in April 1981 is summarized on a monthly basis in Table 1.1 for each of the 9 species and treatment combinations. After the f i r s t growing season a significant difference in survival was found between species (Table 1.2). Untreated Douglas-fir (91%) survived significantly better than untreated western hemlock (38%) which in turn survived significantly better than untreated Pacific silver f i r (23%). Shade cards significantly increased survival in western hemlock by 35% and Pacific silver f i r by 23%, however, survival was not significantly increased in Douglas-fir (8%) because of i t s overall high survival rate. The inclined self-shade treatment had no significant effect in the f i r s t year. After the second growing season, the survival of Douglas-fir (76%) was s t i l l significantly better than western hemlock (15%) and Pacific silver f i r (8%). Treatment effects also became more noticeable. Shade cards significantly increased survival by 21% in Douglas-fir, 44% in western hemlock, and 27% in Pacific silver f i r . The inclined self-shade treatment TABLE 1.1 Monthly percentage survival of seedlings which were planted In April 1981. Date DOUGLAS-FIR WESTERN HEMLOCK PACIFIC SILVER FIR Shade Cards Inclined Shade Unshaded Control Total Shade Cards Inclined Shade Unshaded Control Total Shade Cards Inclined Shade Unshaded Control Total Apr 81 100 100 100 100 100 100 100 100 100 100 100 100 May 81 99 98 99 99 99 96 97 98 95 95 94 95 June 81 99 97 98 98 96 94 95 95 91 94 92 93 July 81 99 95 94 96 93 92 90 92 88 87 88 88 Aug 81 99 95 94 96 88 84 84 85 81 77 79 79 Sept 81 98 95 94 96 78 52 43 58 53 34 29 38 Oct 81 98 94 91 94 73 45 38 52 46 25 23 31 May 82 98 91 88 92 63 25 25 38 38 13 11 21 June 82 97 89 79 88 63 25 21 35 38 13 9 20 July 82 97 89 78 87 62 25 19 34 38 10 9 19 Aug 82 97 89 78 87 62 25 19 34 35 10 8 17 Oct 82 97 89 76 87 59 25 15 33 35 10 8 17 -/2_ TABLE 1.2 Percent survival in October of 1981 and 1982 of seedlings planted in April 1981. Numbers followed by the same letter are not significantly different (t = 0.05) according to Duncan's New Multiple Range Test. October 1981 TREATMENT Douglas-Fir Western Hemlock Pacific Silver F i r CONTROL UNTREATED 91 (a) 38 (c) 23 (d) INCLINED STEM SELF SHADED 94 (a) 45 (c) 25 (d) SHADE CARDS ARTIFICIAL SHADE 98 (a) 73 (b) 46 (c) October 1982 > TREATMENT Douglas-Fir Western Hemlock Pacific Silver F i r CONTROL UNTREATED 76 (b) 15 (f) 8 (f) INCLINED STEM SELF SHADED 89 (a) 25 (e) 10 (f) SHADE CARDS ARTIFICIAL SHADE 97 (a) 59 (c) 35 (d) -73" significantly increased survival by 13% in Douglas-fir and 10% in western hemlock; however, survival was increased by only 2% in Pacific silver f i r . By the second year seedlings in the inclined treatment had responded geo-tropically, so there was no practical, difference between this treatment and the control. This suggests that the significant effect of the inclined treatment was fortuitous rather than due to shading. B. Survival of Seedlings Planted i n 1982 Survival of seedlings planted in April 1982 is summarized on a monthly basis in Table 1.3 for each of the 12 species and treatment combinations. At the end of the growing season, trends similar to those observed in 1981 were seen in untreated and shade card treated Douglas-fir and western hem-lock. In contrast, Pacific silver f i r survival was dramatically improved. Survival of untreated Douglas-fir, western hemlock, and Pacific silver f i r seedlings was 90%, 42%, and 69% respectively. Shade cards increased survival by 5%, 53%, and 21% respectively. Irrigation increased survival by 4%, 26%, and 24% respectively. The combination of irrigation and shade cards Increased survival by 0%, 48%, and 29% respectively. The only significant difference between treatments was found in western hemlock (Table 1.4). Other treatment differences were i n s i g n i f i -cant because of the higher overall survival rate and the high standard deviation between blocks where survival was lower. C. The Effect of Weather on Seedling Survival Weather was found to be a major factor influencing seedling survival at the experimental site. The highest rates of seedling mortality occurred TABLE 1.3 Monthly percentage s u r v i v a l of seedlings which were planted l n A p r i l 1982. D O U G L A S - F I R W E S T E R N H E M L O C K P A C 1 : F I C S I L V E R I ' I R Date Unshaded Control Shade Cards I r r i g a t i o n Shade & I r r i g a t i o n Unshaded Control Shade Cards I r r i g a t i o n Shade & I r r i g a t i o n Unshaded Control Shade Cards I r r i g a t i o n Shade & I r r i g a t i o n Hay 82 100 100 100 100 100 100 100 100 100 100 100 100 June 82 96 99 96 93 96 95 93 90 98 91 93 98 July 82 87 98 96 90 75 95 80 90 94 91 93 98 Aug 82 83 95 95 90 42 95 68 90 79 91 93 98 Sept 82 83 95 94 90 42 95 68 90 74 90 93 98 Oct 82 90 95 94 90 42 95 68 90 69 90 93 98 -J5-TABLE 1.4 Percent survival in October of 1982 of seedlings which were planted in April 1982. Numbers followed by the same letter are not significantly different (t = 0.05) according to Duncan's New Multiple Range Test. October 1982 TREATMENT Douglas-Fir Western Hemlock Pacific Silver F i r CONTROL UNTREATED 90 (a) 42.(b) 69 (a) SHADE CARDS 95 (a) 95 (a) 90 (a) IRRIGATION 94 (a) 68 (a) 93 (a) SHADE CARDS AND IRRIGATION 90 (a) 90 (a) 98 (a) '16-during extended periods of hot dry weather when the s o i l profile water content was low. The 1981 growing season was characterized by cool cloudy wet weather from the time the seedlings were planted until mid-July. From mid-July unt i l September less than 10 mm of rain f e l l and exceptionally hot clear weather occurred. Mid-day maximum air temperatures ranged from 25-37°C during this period. A summary of precipitation, solar radiation, air temp-erature, and relative humidity during the 1981 growing season is shown in Figure 1.1. The course of s o i l profile water contents in 1981 showed very l i t t l e change in root zone water storage during the early part of the summer because of the high r a i n f a l l . From mid-July to July 27 a rapid drawdown occurred in the top 0.4 m (Figure 1.2). This depletion continued at a slower rate during the hot weather in August unt i l minimum s o i l water contents of <5% by dry s o i l weight were measured in the top 0.1 m. Changes in storage were observed to a depth of 0.75 m. Soil water potentials in the seedling root zone (0-0.15 m) were found to be as low as -1.8 MPa during this period. Very l i t t l e measurable difference in s o i l water poten-t i a l was found between the untreated, shaded, and inclined treatments. Daytime seedling water potentials as low as -3.4 MPa were also measured during this period (Livingston and Stathers 1981). The majority of seedling mortality occurred during the hot weather and dry s o i l conditions in August. The seedling assessment after this period showed a 41% decline in untreated western hemlock survival and a 50% decline in untreated Pacific silver f i r survival (Table 1.1). Douglas-fir withstood these harsh conditions very well and i t s survival was reduced by only a few percent. -11' PRECIPITATION (MM) ll, , I.I J 60' eo RELATIVE HUMIDITY (%) SO 20-10 TOTAL SOLAR RADIATION (MT/ M* ) 40 30 20 AVERAGE AIR TEMPERATURE C O 10 « ; T* 19 24 29 4 9 1 4 19 24 29 4 9 14 19 JUNE JULY AUGUST FIGURE 1.1 Da l ly values of p r e c i p i t a t i o n , solar i r rad iance , average a l r temperature, and average r e l a t i v e humidity during the 1981 growing season. -18-"I [ 1 1 T L....._""l AU6..6. r - J .1 jJ~ JUL 27 rl r 0.2 JE 0.4 Q. UJ O LJ _J fe 0.6 DC CL AUG 14 r [ j APR 24 (field capacity ) O CO 0.8 _ L _ L 5 10 15 20 GRAVIMETRIC SOIL WATER CONTENT (%) 25 FIGURE 1.2 Profiles of gravimetric s o i l water content at selected times during the summer of 1981. '79' During the winter of 1981-82, seedling mortality did not change dramatically in the Douglas-fir treatments. Western hemlock and Pacific silver f i r survival were reduced by about 10%. The summer of 1982 was not as intensely hot and dry as in 1981. A summary of the weather during the 1982 growing season is given in Figure 1.3. Cool overcast conditions occurred through most of the summer with the exception of a week in June and two weeks in mid-August. In June, the s o i l water content was s t i l l high enough to buffer the effects of the warm clear weather. Less than 30 mm of rain f e l l during June, July and August, so the clear hot weather in August caused more seedling stress because of the dry s o i l profile. Even though considerably more rain f e l l in 1981 than 1982 so i l water contents were similar in August of both years because of the r a i n f a l l distribution, the low water holding capacity of the s o i l , and the different distribution of potential evapotranspiration rates in each summer. Measured root zone water potentials were also similar in August of both years. The duration of stressful conditions was much shorter in 1982 than in 1981. Even though s o i l water contents were similar, the overcast skies in the summer of 1982 resulted in very few days when surface s o i l temperatures reached lethally high levels. The higher seedling survival rate in 1982 reflects these milder conditions. The fact that s o i l water contents and potentials in the seedling root zone were very similar in August of each year, while at the same time the intensity and duration of high surface temperatures was higher in 1981 when significantly more mortality was observed, suggests that thermal stress was the major cause of seedling mortality for shade tolerant western hemlock and Pacific silver - 2 0 -FIGURE 1.3 Daily values of precipitation, solar irradiance, average air temperature, and average relative humidity during the 1982 growing season. f i r seedlings. The results of the 1982 experiment (Table 1.4) provide further evi-dence of this. The null hypothesis of the 1982 experiment was that there would be no difference in survival between treatments. There was evidence to disprove the null hypothesis for western hemlock, i.e. a significant improvement in survival was obtained by shading, irrigating, and combined shading and irrigating. If water stress had been the sole cause of mortal-ity, then the unwatered treatments should have shown lower survival; how-ever, the shade treated seedlings survived as well as the seedlings in the 2 treatments which received water. If thermal stress had been the sole cause of mortality, then unshaded treatments should have shown lower sur-vival. The unshaded-unwatered treatment did show significantly lower survival, but the unshaded-watered treatment did not; however, surface s o i l temperatures did not reach lethal levels in the irrigated treatment because of the wetter s o i l . This suggests that thermal stress was the major cause of seedling mortality. D. The Effects of Shade and Irrigation on S o i l Temperatures Irrigation, shade, and combined irrigation and shade treatments were shown to significantly increase seedling survival at the experimental site. Shade cards were also shown to be more effective as the severity of drought increased. In this section the effects of each of these treatments on s o i l temp-eratures is illustrated using measurements collected on August 7, 1982. This particular day was representative of summer drought conditions at the experimental site. The sky was clear, the air was warm, and the s o i l pro--22-f i l e was very dry. Similar conditions had been observed for the previous week. The diurnal course of air temperature at the 1.0 m height, s o i l sur-face temperature (measured with a fine wire thermocouple), and s o i l temp-eratures at several depths below the untreated control treatment are shown in Figure 1.4. Although air temperature only reached 26.5°C, the surface temperature remained above 50°C for more than 7 hours. A maximum surface temperature of 64°C was recorded at 1500 P.S.T. At the 0.05 m depth, the maximum temperature was 32.5°C. At the 0.15 m depth, the minimum and maxi-mum s o i l temperatures were 15 and 20°C. At the 0.3 and 0.5 m depths, the temperatures remained constant at approximately 13 and 15°C respectively. These results show a strong contrast between temperatures at the sur-face and those below 0.15 m. Deep s o i l temperatures were low for the opti-mal growth of Douglas-fir seedlings (Brix 1967; Heninger and White 1974); however, at the same, time surface temperatures were lethally high (Maquire 1955; Silen 1960; Baker 1929; Cleary 1978). Operational procedures which would increase the heat flow from the surface into the profile would be beneficial on this site. This would reduce the daytime surface maximum temperature, increase the nighttime surface minimum temperature, and increase the mean s o i l temperature at depth. Factors affecting the s o i l thermal regime are discussed in more detail in Chapter 2. The diurnal course of air temperature and surface temperature for the untreated control, shade, irrigation, and shade and irrigation combined treatments are shown in Figure 1.5. The combined shade and irrigation treatment gave the greatest surface temperature reduction; about 10°C less -23-- — i 1 1 " I AIR T E M P E R A T U R E ( 1.0 m height) SURFACE TIME ( hr* P.S.T. ) FIGURE 1.4 The diurnal course of s o i l temperature at several depths, surface s o i l temperature (measured with a fine (0.075 mm diameter) wire copper constantan thermo-couple) and air temperature at the 1 m height on August 7, 1982. r 1 1 1 1 1 > r TIME P.S.T. (hrs) FIGURE 1.5 The diurnal course of surface s o i l temperature adjacent to shaded, irrigated, shaded and irrigated, and control treatments on August 7, 1982. The course of air tempera-ture at the 1 m height is also shown. -25-than air temperature and up to 50°C less than the untreated control. The irrigation treatment reduced midday surface temperatures by up to 35°C and the shade card treatment by up to 40°C. The dramatic rise in temperature after 1500 P.S.T. in the shaded treatments was a consequence of the shade card's shadow moving to the east of the seedling. This rapid temperature rise even at this late hour in the day indicates the importance of the correct shade card placement for effective shading. The corresponding seedling root zone temperatures for each of these treatments are shown in Figure 1.6. Very l i t t l e temperature difference was measured between treatments at depths below 0.05 m, even though very large differences were found at the surface. This i s explained by the fact that the treatments affected a relatively small proportion of the so i l surface, and therefore the treatment effect was overwhelmed at depth. A l l of the higher temperature amelioration treatments were effective in keeping surface temperatures below the c r i t i c a l limit of 54°C (Cleary 1978; Silen 1960) during periods of hot weather when the so i l was dry. This is reflected in the seedling survival statistics; particularly the 1981 t r i a l where significant differences in survival were seen between untreated and shaded seedlings. 1.5 CONCLUSIONS Very low survival of planted 1-0 containerized seedling stock was observed on a steep south facing high elevation clearcut on Vancouver Island. At the end of the summer with harsh weather conditions, survival of untreated Douglas-fir seedlings was found to be significantly higher than that of untreated western hemlock and Pacific silver f i r seedlings. -Z6-E o CL UJ Q O CC CL CO 1 — 1 1 1 1 1 1 1 1 \ T~7]~~ N 0 TREATMENT " \ / / / SURFACE T E M P E R A T U R E V / HOURS ABOVE 5 0 °C - 7 ' HOURS A B O V E 30°C - 9 " . i i 1 1 1 1 — ; *S 1 1 1 i 20 1 7 1 3 *~~yr7r^ NO SHADE - WATER i*r if < i »# if * / •i i / \ / / SURFACE T E M P E R A T U R E / / HOURS ABOVE 50°C - 0 ' * * HOURS A B O V E 30°C - 3 i | l i 1 1 1 1 — - V . 1 \ / 1 1 1* 1 I I I 1 1 ' ' ^ S H A D E " WATER SX M /I SURFACE T E M P E R A T U R E / / HOURS A B O V E 50°C - 0 * ** HOURS A B O V E 30°C - 1 • • _ i 1 — 1 1 1 1 1 W / £ ' SHADE - NO WATER \ I II SURFACE T E M P E R A T U R E \ 1/ HOURS ABOVE 5 0 ° C - 0 5 - ** ** HOURS A B O V E 30°C 6 • t 1 1 J-SOIL TEMPERATURE ( ° C ) FIGURE 1.6 The d iurnal course of s o i l temperature i n the seedling root zone adjacent to shaded, i r r i g a t e d , shaded and i r r i -gated, and control treatments on August 7, 1982. Numbers above the p r o f i l e refer to the time of day ( P . S . T . ) . -27-This was also true for these seedlings at the end of the second milder growing season. For seedlings planted at the beginning of the milder grow-ing season, untreated Douglas-fir and P a c i f i c s i l v e r f i r seedlings survived s i g n i f i c a n t l y better than the untreated western hemlock seedlings. During the harsh summer, shade cards s i g n i f i c a n t l y increased s u r v i v a l of western hemlock and P a c i f i c s i l v e r f i r seedlings; however, a f t e r the second growing season following planting, shade cards s i g n i f i c a n t l y increased the s u r v i v a l of a l l three species. During the harsh summer, an i n c l i n e d self-shade treatment only increased s u r v i v a l marginally. For seedlings planted at the beginning of the milder growing season, shade cards, t r i c k l e i r r i g a t i o n and the combination of these two treatments only s i g n i f i c a n t l y increased west-ern hemlock s u r v i v a l . The majority of seedling mortality occurred during hot droughty periods i n the summer when root zone s o i l was very dry. It was t e n t a t i v e l y concluded that thermal stress due to very high surface s o i l and root seed-l i n g c o l l a r temperatures was the major cause of mo r t a l i t y . -£8~ CHAPTER 2 FACTORS AFFECTING SURFACE SOIL TEMPERATURE IN A FOREST CLEARC UT FACTORS AFFECTING SURFACE SOIL TEMPERATURE IN A FOREST CLEARCUT 2.1 INTRODUCTION In Chapter 1 coniferous seedling mortality on a high elevation, south facing forest clearcut on Vancouver Island, B.C. was found to be associated with high surface s o i l temperatures during periods of summer drought. The objective of this chapter i s to discuss the factors which cause the high surface s o i l temperatures that ultimately lead to seedling mortal-i t y and to c a l c u l a t e the e f f e c t of these factors on surface temperature i n a forest c l e a r c u t . A better understanding of the physical processes which determine the surface s o i l temperature w i l l aid i n the development of ameliorative treatments that can improve seedling s u r v i v a l under harsh weather and s o i l conditions. To meet t h i s objective, the physics describing the combined r a d i a t i o n and energy balances of a surface are outlined. The aerodynamic approach to estimating sensible and latent heat fluxes i s also discussed. This theory i s applied i n the following chapter where i t i s used i n conjunction with a s o i l heat flow model to estimate surface and subsurface s o i l temperatures from e a s i l y measured meteorological data. 2.2 THEORY A. The Energy and Radiation Balance of a Forest Clearcut The energy balance of a s o i l surface can be written as follows: R n - H - LE - G Q = 0 (1) where Rn, H, LE and G Q are the net r a d i a t i o n , sensible, l a t e n t , and -30-s o i l heat f l u x densities respectively at the surface. The net r a d i a t i o n f l u x density can be expressed as the sum of the net solar and long wave irradiances as follows: R N = Rg(l - a) + e S ( R L - aT!sk) (2) where R S i s the solar irradiance, a i s the solar r e f l e c t i o n c o e f f i c i e n t (albedo), e s i s the long wave emissivity of the surface, a i s the Stephan-Boltzmann constant, T s i s the s o i l surface temperature i n K, and R L i s the long wave irradiance from the sky. RL can be estimated as follows: RL = £a a T a k (3) where T a i s the a i r temperature i n K at screen height and e a i s the e f f e c t i v e emissivity of the sky. The emissivity i s a function of wavelength; however, i t can be treated as a constant equal to some average value for f a i r l y large wavebands (Campbell 1977). Under clear skies, e a i s strongly dependent upon the vapour density of the atmosphere which i n turn i s well correlated with a i r temperature. The clear sky emissivity ( e a c ) i s c l o s e l y approximated as follows (Idso and Jackson 1969): e a c = 1 - 0.261 exp (-7.77 x 10~^ T a 2 ) (4) where T a i s i n °C. The atmospheric emissivity normally ranges from 0.75 to 0.85 under clear skies. Under cloudy skies the atmospheric emissivity i s higher be-cause the emissivity of clouds i s nearly 1.0. The atmospheric emissivity on cloudy days can be estimated by adding the energy emitted by the clear portions of the sky to the energy emitted by the clouds as follows - 3 1 -(Campbell 1977): e a " eac + c C1 ~ eac " 4AT/Ta) (5) where C is the fraction of sky covered by cloud, AT i s the temperature difference between screen height and the cloud base air temperature and T a is in K. From inspection of (5) i t can be seen that estimates of e a under cloudy conditions are quite insensitive (±10%) to errors in the e s t i -mation of C and AT. Note when C = 0 (ie. clear skies), (5) reduces to e a = e . -ac By combining equations (1) and (2) as follows: R s(l-ct) + e s (R L - oTg^) - H - LE - G D = 0 (6) i t becomes evident that the absorbed solar and long wave radiation are dis-sipated at the s o i l surface as either emitted long wave radiation, sensible heat, latent heat, or s o i l heat. Upon rearrangement of (6) the s o i l surface temperature can be expres-sed as follows: T s = [R s ( l - a) + es RL - H - LE - G Q]/e s a ) 0 , 2 5 (7) B. The Soi l Heat Flux The surface s o i l heat flux density can be found from Fourier's Law of heat conduction as follows: G D = -k s (3T/3z)| z = o (8) where k s is the apparent thermal conductivity of the s o i l (which includes the combined effects of pure conduction and vapour d i s t i l l a t i o n due to temperature gradients), z i s s o i l depth, and 3T/8z|z = o is the s o i l - 3 2 -temperature gradient at the surface. The negative sign indicates that the flux is directed into the s o i l (G 0 > 0) when 3T/3z < 0 (i.e. temperature i s decreasing with depth). Soil heat flow is quite accurately described by (8) because radiative and convective heat transfer are negligible in the s o i l . Heat transfer due to vapour movement is similarly insignificant because the thermal and isothermal vapour fluxes tend to cancel during the day and are very small at night (Kimball et a l . 1976). A quantitative method of calculating GQ from a numerical solution of the heat flow equation w i l l be presented in the following chapter; however, to a f i r s t approximation, daytime average G Q is reasonably approximately by 15% of daytime average Rs, while during the nighttime GQ is well approximated by Rn (Novak 1981). C. Aerodynamic Equations Describing the Sensible Heat Flux, Latent Heat  Flux, and Momentum Flux i n the Turbulent Boundary Layer Sensible and latent heat transfer occur predominantly by convection in the turbulent boundary layer. Therefore, they are closely related to the momentum flux to the surface. Equations describing the sensible heat, latent heat, and momentum flux densities are respectively: where K^ , Ky and are the eddy di f f u s i v i t i e s for heat, vapour and momentum respectively, pa is the density of air, Cpa is the specific heat of air, y is the psychrometric constant, e is the vapour pressure of water in air, u is the windspeed, and z is the height above the ground. The sign H = - p a c p a K H (3T/3z) LE = - ( p a c p a/y) K v (3e/3z) T = p a K M (3u/3z) (10) ( I D (9) -33" convention used above is determined from the direction of transfer. During the day the gradients are negative (lapse) and H and LE are directed away from the surface (positive fluxes). At night the gradients are positive (inversion) and H and LE are directed toward the surface (negative fluxes). The momentum flux, also known as the shearing stress, is always directed toward the surface. Using aerodynamic theory, the eddy di f f u s i v i t i e s in (9), (10) and (11) are defined as follows: K H = k u* z/>H (12) K v = k u^ z/(|)V (13) % = k u^ z/cJ)M (14) where k i s von Karman's constant, u^ (= ( T / p a ) ^ " 5 ) is the f r i c t i o n velocity (the tangential velocity of an eddy in the turbulent boundary layer (Thorn 1975)), and <j>H, <[>y and cj)M are the diabatic influence functions for heat, vapour and momentum respectively. The diabatic i n f l u -ence functions account for the effects of atmospheric stability on the turbulent transfer of H, LE, and T . The atmosphere is in a state of neutral equilibrium when the vertical temperature gradient is exactly equal to the rate of cooling (due to expan-sion) of a parcel of air rising adiabatically (-0.01°C m - 1). During the day, when the s o i l surface is warmed by solar radiation, the air tempera-ture near the surface decreases rapidly with height. The atmosphere is then in a state of unstable equilibrium because parcels of air which are displaced upward continue to ascend since they are less dense than the surrounding a i r . At night when the s o i l surface temperature is lower than the air temperature because of longwave radiative cooling, buoyancy forces - 3 4 -oppose any original air displacement and the atmosphere is considered to be in a state of stable equilibrium. In unstable conditions turbulence is enhanced by buoyancy and in stable conditions i t is suppressed. Under neutral conditions the diabatic influence functions equal 1, under unstable conditions they are less than 1, and under stable conditions they are greater than 1 (Webb 1965). The state of the atmosphere can be described quantitatively by the relation between the production of turbulent kinetic energy by buoyancy and mechanical forces. This i s described by the turbulent kinetic energy balance equation as follows: — — g H e " p a 3z + p a c p a T a u ^ where g is the gravitational acceleration, e is the viscous dissipation of kinetic energy and T a is in degrees K. The f i r s t term in this equation describes the mechanical production of turbulent kinetic energy and the second term describes the buoyant production of turbulent kinetic energy. The ratio of the buoyant to mechanical turbulent kinetic energy pro-duction terms is used as the basis of two different indices of atmospheric stability; the Richardson number, Ri, (Richardson 1920) and the Monin-Obukhov scaling or stability length, L, (Monin and Obukhov 1954). The gradient Richardson number, derived from (15), can be written in differential form as follows: Ri = (g/T a) (3T/3z) / (3u/3z) 2 (16) and is calculated in f i n i t e difference form as follows: Ri = (g/T a)[(T s-T a)/u 2] [ln (z/z 0)]z (17) where T a is the mean air temperature between the measurement heights z -35-and z Q. The surface roughness length, z Q, is the distance above the ground where the wind speed extrapolates to zero when plotted against the logarithm of height. It can also be thought of as the size of the smallest eddy in the turbulent boundary layer (Thorn 1975). The Monin-Obukhov scaling length, also derived from (15), is calculat-ed as follows: L = (p a c p a T a u*3) / (k g H). (18) L is defined in terms of the sensible heat flux and is therefore independ-ent of height whereas Ri is defined in terms of the air temperature gradi-ent and varies with height. A more commonly used form of L is the s t a b i l i -ty index, £ (= z/L), which has been found to be related to the Richardson number as follows (Dyer and Hicks 1970): Ri = c (W^M2)- (19) Measurements by Dyer and Hicks suggest that <$>y[ ) is nearly unity (<|>jl/(j)M = 1.00 ± 0.06) over a very wide range of atmospheric stab i l -i t y , which implies that Ri « £. This has also been suggested by Thorn et al,(1975). Paulson (1970) and Businger (1972, 1975) suggest that Ri = c, for unstable conditions (£ or Ri < 0). For stable conditions (£ or Ri > 0) they suggest that the relationship between Ri and £ is approximated as follows: C - Ri / (1 - 4.7 Ri). (20) Barker and Baxter (1975) used a similar method of calculating £ from Ri. The indices of stability are used to calculate the diabatic influence functions as outlined below. It should be noted that these are semi-empir-i c a l relationships which have been derived from the analysis of many large - 3 6 -data sets as well as dimensional analysis. Physically quantitative relationships describing turbulence do not yet exist. For unstable atmospheric conditions the relationship between the dia-batic influence functions and £ is described as follows (Businger (1966); Pandolfo (1966); Dyer and Hicks (1970)): *H * *V = <f>M2 = (1-16C)" 0 - 5 : Z < 0. (21) For stable atmospheric conditions, Businger et a l . (1971) found that: *H ™ *V = *M = ( ° - 7 4 + 4 ' 7 5) : C > 0. (22) The applicability of (21) i s questionable under extremely unstable condi-tions where free convection is the dominant mode of heat transfer (Dyer and Hicks 1970). Many different methods of incorporating diabatic influence functions in the aerodynamic equations have appeared in the literature in the past few decades. There are basically three different approaches. Two of these are integrated profile methods and the third is based upon the ratios of the eddy d i f f u s i v i t i e s . The earliest of the integrated profile methods was introduced by Lettau (1962) and was used by Fuchs et a l . (1968) to determine the sensible heat flux density. Assuming similarity between momentum and sensible heat transfer (KJJ = K^) Fuchs et a l . calculated the sensible heat flux density as follows: H = P a C p a k 2 u (T s - T a) / ln(z/z Q) - lfa) 2 (23) where ty^, the integrated diabatic profile correction factor for momentum transfer is defined as: -37-z *M = / (4»M " 1) (z + z 0 ) _ 1 d z - (24) o Lettau numerically integrated (24) using: • M = (1 - 18 R i ) " 0 - 2 5 (25) as suggested by Panofsky, Blackadar, and McVehil (1960. The second integrated diabatic profile approach, developed by Paulson (1970), differs from Lettau's method in that separate diabatic influence functions are used for heat and momentum. By substituting (12) into (9) and integrating with respect to z from z Q to z, the sensible heat flux density is expressed as follows: H = Pa cpa k u* ( Ts" Ta) / ln(z/z Q) - ^ H ) (26) where T |» h = / (1 - $H) r 1 dr.. (27) o The f r i c t i o n velocity i s obtained similarly by combining (11) and (14) to give: = (6u/6z) k(z - z Q) / <|>M (28) and integrating from z Q to z to give: u 4 = k u / In (z/z 0) - I } J m (29) where C Tpn = J (1 - <|)M) V, d? , (3 0) o Integrated solutions of (27) and (30) i n terms of x, were obtained for unstable atmospheric conditions (using 21) by Paulson (1970) as follows: f H = 2 In [(1+x2) 12] : c < 0 (31) - 3 8 -ifa = 2 £n [(1+x) / 2] + In [(1+x2) / 2] - 2tan _ 1x + tr/2 : c. < 0 (32) where x = (1 - 16?)°' 2 5. (33) For stable atmospheric conditions Businger (1975) suggested that the diabatic profile correction factors for heat and momentum were identical and could be calculated as follows: = = 4 - 7 5 : 5 > 0. (34) There are several variations on each of the two integrated diabatic profile methods quoted in the literature. Some of the more notable recent work in this f i e l d has been done by Berkowicz and Prahm (1982), Hammel et a l . (1982), Businger (1972, 1975), and Barker and Baxter (1982). Berkowic and Prahm used (31), (32), and (33) to estimate H from routine meteorologi cal data. Hammel et a l . assumed ipjj = ipjj and used (31) and (33) to model s o i l temperatures in a wheat f i e l d . It should be noted that both Businger and Hammel et a l . have lef t out the multiplication factor of 2 in (31). This was probably a typographical error since integration of (27) shows that Paulson's work is correct. Berkowicz and Prahm also confirm this. The third basic approach used to calculate sensible and latent heat fluxes was developed independently by Holbo (1971) and Thorn et a l . (1975). It was used by Leckie et a l . (1981) to estimate sensible heat flux from remotely sensed surface temperature. This non-integrated profile method i based upon the stability dependent multiplicative correction factor, "F", which is incorporated in the basic sensible heat flux equation as follows: H = p a c p a k 2 u (T s - T a) F / ln (z/z Q) 2 . (35) -39-The correction factor F is defined as follows: F = (chH C^M)-1 = (R n - GQ) / (H'+LE') (36) where H' and LE' are the sensible and latent heat flux densities calculated assuming neutral atmospheric conditions. Using the same diabatic influence functions from equations (21) (with Ri * ?) and (22), Thorn (1975) has shown that: F = ( 1 - 1 6 R i ) 0 - 7 5 : Ri < 0 (37) and F = (1 - 5 R i ) 2 : Ri > 0 (38) whereas Holbo (1971) found that for a bare pumice s o i l in Oregon, U.S.A.: F = ( 1 - 3 4 R i ) 0 - 5 5 : Ri < 0 (39) which Thorn et a l . (1975) claim differs from their F values by less than 10% through the moderately unstable range (-0.6 < Ri < 0). The integrated diabatic correction factors used in a l l of the inte-grated methods are shown in Figure 2.1. The diabatic influence functions used by Thorn et al.(1975) in their non-integrated method are identical to those used by Paulson (1970). It is of interest to note that Lettau's function = Fuchs et al.'s ^ function « Businger's and Hammel et al.'s i^ H function = Paulson's ipg function * 2. The latent heat flux is defined similarly in a l l 3 methods since ij'H = <()Y and Kg = Ky. Using the integrated diabatic profile method of Paulson (1970) the latent heat flux density is calculated as follows: LE = p a c p a k u^ (e Q - e) / [ln (z/z 0) - i|>v] y (40) where e Q - e is the vapour pressure difference between the s o i l surface and the atmosphere. This equation is easily applied when the s o i l i s moist - 4 0 --2.5 -2.0 -1.5 INDEX OF STABILITY 1.0 -0.5 ( Rl or £ ) FIGURE 2 . 1 Integrated diabatic profile correction factors for heat and momentum versus 6 or Ri (indices of atmospheric stability) for unstable, neutral, and stable atmospheric conditions. since e Q = e*(T s), i.e. the saturation vapour pressure at surface temperature, but i t becomes very d i f f i c u l t to estimate e Q as the s o i l surface dries. 2.3 RESULTS AND DISCUSSION A. Surface Energy Balance and Clearcut Microclimate The energy balance of a forest clearcut is a function of the net solar and longwave irradiances at the surface, the thermal and moisture characteristics of the s o i l profile, the atmospheric gradients of temperature, vapour pressure and wind, and the surface resistance to heat and vapour exchange. a) The Radiation Balance The net solar and longwave irradiances have the greatest single influence upon the microclimate of a forest clearcut. Solar radia-tion is highly variable in intensity throughout the year, with l a t i -tude, and with slope and aspect. South-facing slopes receive s i g n i f i -cantly more radiation than other aspects in the northern hemisphere, particularly when the solar elevation angle is low in the winter. Increasing slope angle increases the irradiance on south-facing slopes and decreases i t on north-facing slopes. The total daily irradiance on east and west facing slopes is similar, but the irradiance on the west facing slope lags the east-facing slope by a few hours. Solar radiation is attentuated and backscattered by clouds and particulate matter in the atmosphere. In British Columbia there are typically -4Z-fewer hours of sunshine in June than July because June usually has more cloudy days (Hay 1979). The solar reflection coefficient (albedo) of the surface i n f l u -ences the net solar irradiance and can quite easily be used to manipu-late the energy balance. Typical albedos of soils range from 0.06 to about 0.6 depending on the angle of the sun, the colour of the surface and the s o i l water content. The net longwave irradiance depends upon the temperatures and emissivities of the sky and the surface. Surface emissivities of s o i l usually range from 0.95 to 0.98, increasing with s o i l water content. Atmospheric emissivities range from about 0.75 to 1.0 depending upon sky cloud cover and the temperature and humidity of the lower atmos-phere . b) The So i l Thermal Regime The s o i l thermal regime is dependent upon (1) the surface heat flux density; G Q (the quantity of heat passing through a unit area of surface s o i l per unit time), (2) the s o i l thermal conductivity; k s (the s o i l heat flux density per unit temperature gradient), and (3) the volumetric heat capacity of the s o i l ; C s (the change in heat content of a unit s o i l volume required to cause a 1°C temperature change). The s o i l thermal conductivity and volumetric heat capacity have a very significant influence upon s o i l profile temperatures. Both of these thermal properties are functions of the s o i l bulk density and the relative proportions of the air, organic, mineral, and water - 4 3 -fractions of the s o i l . Table 2.1 shows the values of the thermal properties of each of these materials. Since the thermal conductivity and volumetric heat capacity of water are so much higher than the other s o i l constituents, particularly air, the s o i l water content has a large influence upon the s o i l thermal regime. The values of both thermal properties for organic s o i l are also much lower than for mineral s o i l and as a result the thermal response of a layered forest s o i l can be significantly modified by removing the upper organic horizon. c) The So i l Water Regime Soil water indirectly influences the surface energy balance in 3 separate ways. As mentioned above, the values of the thermal conduct-i v i t y and volumetric heat capacity of water are much higher than those of a i r . As the s o i l drains, pore water is replaced by air and conse-quently the thermal conductivity and volumetric heat capacity of the s o i l decrease. The latent heat flux density is also dependent upon the surface s o i l water content. During hot dry periods when there is a large evaporative demand, the surface s o i l dries, forming a mulch which restricts vapour movement. The water content of the s o i l profile and the s o i l hydraulic properties determine the supply of water to the surface. The third effect occurs when precipitation having a temperature different from that of the s o i l percolates through the s o i l profile. This is an example of heat transport resulting from mass flow. TABLE 2.1 Thermal properties of natural materials (after DeVries 1963) Material Remarks Density Mg m-3 C MJ m-3 C _ 1 k W m-1. C - 1 Water 4°C s t i l l 1.00 4.18 0.57 Air 10°C s t i l l 0.0012 0.0012 0.025 turbulent 0.0012 0.0012 «125 Sandy Soil dry 1.6 1.28 0.3 40% pore space saturated 2.0 2.96 2.2 Clay Soil dry 1.6 1.42 0.25 40% pore space saturated 2.0 3.10 1.58 Peat Soil dry 0.3 0.58 0.06 80% pore space saturated 1.1 4.02 0.50 -45-d) Boundary Layer Resistance to Heat, Water Vapour and Momentum  Transfer Heat, water vapour, and momentum transfer are strongly influenced by the resistance of the laminar and turbulent boundary layers of the clearcut surface. Transfer in the laminar boundary layer occurs by molecular diffusion. Beyond the laminar boundary layer turbulence causes convective transfer to become the dominant transfer process. The resistance of the turbulent boundary layer depends upon the aerodynamic roughness of the surface and the wind speed. Extended flat surfaces which are aerodynamically smooth develop relatively thick boundary layers. As the height of the roughness elements increases eddies induced by turbulence near the surface break down the boundary layer. The aerodynamic roughness lengths of various natural surfaces are given in Table 2.2. As the wind velocity increases, turbulent mixing increases heat and vapour transfer. The air tempera-ture, water vapour pressure and wind speed directly above the surface depend to a certain extent on the nature of the local surface; how-ever, they are also a function of macroclimate and upstream weather processes. B. Relationship Between Surface Energy Balance Components and Surface  Temperature Equation (7) can be used to study the relationship between the surface energy balance and surface temperature. This is illustrated in Figure 2.2 where s o i l surface temperatures are shown for a range of -46" TABLE 2.2 Aerodynamic properties of natural surfaces (after Oke 1978). Surface Remarks Roughness length (m) Water S t i l l - open sea 0.1 - 10.0 x 10_ -5 Ice smooth 0.1 x 10~h Snow 0.5 - 10.0 x 10" Sand-desert 0.0003 Soils 0.001 - 0.01 Grass 0.02 - 0.1 m 0.25 - 1.0 m t a l l t a l l 0.003 - 0.01 0.04 - 0.10 Agricultural Crops 0.04 - 0.2 Orchards 0.5 - 1.0 Forests 1.0 - 6.0 Sources: Sutton (1953), Szeicz et a l (1969), Kraus (1972). -47-H + L E + G (Wm-*) SOLAR IRRADIANCE (Wm* 1) FIGURE 2.2 Surface s o i l temperature versus solar irradiance c a l c u l a t -ed for a range of values of H + LE + G Q from energy balance theory (equation 2.7) assuming o = 0.18, e 8 -0.97, e a c - 0.85, and T a = 25°C. solar irradiances from 0-1100 W m when the combined sensible, latent and so i l heat flux densities increase from 0 to 600 W m . A constant air temperature of 25°C, albedo of 0.18, atmospheric emissivity of 0.85 and surface emissivity of 0.97 were used for these calculations. The maximum surface temperature attainable under these conditions is represented by the upper line, where only radiative heat loss occurs. These calculations show that the maximum surface temperatures attainable for solar irradiances of 200, 400, 600, 800, and 1100 W m-2 are 40, 60, 80, 95, and 110°C respect-ively. Even under conditions of an extremely dry, thick LFH surface and calm atmospheric conditions surface temperatures would f a l l short of these values. As the combined sensible, latent, and s o i l heat flux densities increase by 100 W m the surface temperature drops by approximately 10°C under these conditions. The combined total of H, LE, and GQ would rarely be less than 300 W m during a mid-summer day, hence surface temperatures should not exceed about 85°C. C. Comparison of Atmospheric Stability Correction Methods As indicated earlier, there are three methods of stability correcting aerodynamic estimates of the sensible heat flux density. The following comparison of these methods was undertaken because the most recent, scien-t i f i c a l l y rigorous approach suggested by Paulson (1970) requires a complex and time-consuming iterative procedure to estimate separate diabatic cor-rection factors for heat and momentum. Because the other integrated pro-f i l e method, f i r s t derived by Lettau (1962) to describe momentum transfer and later used by Fuchs et a l . (1968) to describe heat transfer, assumes that the diabatic correction factors for heat and momentum are identical i t is much easier to implement. The method developed by Thorn et a l . (1975) - 4 9 -and Holbo (1971) i s appealing because of i t s simplicity and ease of appli-cation; however, because i t is a non-integrated profile method that was developed for use at a single point some distance from the surface where gradients are smaller, its application near a surface subject to very unstable conditions was questioned. The results of a comparison of these methods applied under identical environmental conditions are shown in Figure 2.3. In each case, the following variables remained constant: T a = 25°C, Rs = 1000 W m-2, G + LE = 150 W m-2, a = 0.18, e s = 0.97, e a = 0.85, z = 1.0 m, and z Q = 0.01 m while the wind speed was allowed to drop from 10 to 0.5 m s - 1 so that performance under very unstable atmospheric conditions could be deter-mined. The upper solid line in Figure 2.3 shows the calculated surface temp-erature when sta b i l i t y correction was not used (i.e. (7) solved iteratively using (26) with \|>JJ = 0) . It i s evident from the large difference between the corrected and uncorrected methods that stability correction is import-ant, i.e. that free convection is a very important heat transfer process under unstable conditions. The integrated profile methods of Paulson, Fuchs et a l , and Hammel et al, a l l predict surface temperatures within at least a few degrees of each other over the whole range of wind speeds. The method outlined by Paulson was solved using (7), (26), (29), (31), (32), and (33). Hammel et al.'s method similarly used (7), (26), (29) and (33); however i t assumes that ipj^  = i|>g and that the function equals (31) without the multiplication factor of 2. Fuchs et al.'s method ut i l i z e s (7), (23), (24), and (25) which gives a \JJ m function very similar to that used by Hammel et a l . (as seen in Figure 2.1). The non-integrated methods of -50-I • I - I i i 1 1 1 1 1— 0 2 4 6 8 10 WIND SPEED ( m s-' ) FIGURE 2.3 Surface s o i l temperature versus wind speed calculated using energy balance and aerodynamic theory and the s t a b i l i t y cor rec t ion methods of Paulson (1970), Thom et a l (1975), Holbo (1971), Fuchs et a l (1968), and Hammel et a l (1982). Values were calculated assuming Rg = 1000 W m - 2 , G + LE - 150 W m" 2 , T a - 25°C, o - 0.18, e s = 0.97, e a «• 0.85 and z Q = 0.01 m. The methods of Fuchs et a l and Hammel et a l give i d e n t i c a l resu l t s (see t e x t ) . -5/-Thom et a l . and Holbo were calculated using (7), (35), (37), (38), and (39). The non-integrated profile method used by Thom and Holbo breaks down under very unstable conditions near the surface. At low wind speeds sur-face temperatures were underestimated by up to 20°C. This occurred because the "F" factor increased too quickly for negative Richardson numbers, hence the sensible heat flux was overestimated and as a result the surface temp-erature was seriously underestimated. This method, though well suited to applictions away from the surface, does not accurately predict fluxes near the surfce where gradients are highly non-linear under unstable conditions. A further test of each method involved the comparison of measured and modelled surface temperatures at the experimental site described in Chapter 1. Measurements of surface temperature were taken throughout August 9, 1981 with an infrared thermometer and a fine wire thermocouple. Surface temperatures were modelled using measured air temperature, solar i r r a d i -ance, wind speed, and s o i l heat flux density. The latent heat flux density was estimated to be similar to that measured by Ballard et a l . (1977) on another dry south facing clearcut in southwestern B.C. A 1:1 plot showing the correspondence between modelled and measured surface temperatures using each of the methods outlined above is shown in Figure 2.4. The integrated profile methods of Paulson, Fuchs et a l . , and Hammel et a l . satisfactorily predicted surface temperatures throughout the range of atmospheric stability conditions encountered on this day. The non-integrated diabatic profile methods of Thom et a l . and Holbo success-ful l y predicted temperatures under stable atmospheric conditions, but as in the previous test they underestimated surface temperatures by up to 15°C under unstable atmospheric conditions. -SZ-- 7 0 o LU tr I-< rr _ _ LU 50 LU LU o 2 rr 0 0 o LU _l -J LU Q O 30 10 1 1 1 PAULSON METHOD A TH0M/H0LB0 METHODS o FUCHS ET AL. METHOD — a HAMEL ET AL. METHOD • - • ALL METHODS IDENTICAL- • . A 6 oo hi LINE 10 30 50 70 MEASURED SURFACE TEMPERATURE (°C ) FIGURE 2.4 Calculated versus measured surface s o i l temperatures on August 9, 1981. Calculated values were obtained from measured Rg, T a, u, G Q and LE/R 8 = 0.05 using the s t a b i l i t y c o r r e c t i o n methods of Paulson (1970), Thorn et a l (1975), Fuchs et a l (1968) and Hammel et a l (1982). Surface s o i l temperatures were measured with an i n f r a - r e d thermometer and a fi n e wire thermocouple. -S3' Further micrometeorological studies are required to determine whether Paulson's method can be improved upon for sloping clearcuts. The non-integrated profile method of Thorn et a l . and Holbo is not recommended for this application. D. Examination of Micrometeorological Factors Affecting Surface  Temperature In this section, the effects of changes in solar irradiance (energy input), windspeed, aerodynamic roughness and air temperature on surface temperature are compared. The method of calculation makes use of the Paulson integrated diabatic profile correction procedure tested and recom-mended in the previous section. As indicated earlier, this involves the interative solution of (7), (26), (29), (31), (32) and (33). The results of the calculations are shown in Figures 2.5, 2.6 and 2.7 where the effects of solar irradiance, air temperature, and surface roughness on surface temperature are shown for wind speeds ranging from 0.5 to 10 m s -*. The upper, middle and middle curves of Figures 2.5, 2.6 and 2.7 respectively are the same as the Paulson curve of Figure 2.3. Figure 2.5 shows that higher wind speeds are required to reduce surface temperatures as the solar irradiance increases. Figure 2.6 shows that the surface temperature approaches air temperature as the wind speed increases. This is the result of decreasing turbulent boundary resistance with increasing wind speed. To a good approximation, i f a l l other variables are constant, a rise in air temperature results in the same rise in surface temperature. Figure 2.4 shows that higher wind speeds are required to reduce surface temperatures as the aerodynamic roughness of the surface decreases. The roughness of the -54-80 h o I. 6 0 -LU WIND SPEED (m/s) FIGURE 2.5 Surface s o i l temperature versus wind speed at at the 1 m height for various solar i r rad iances . Values were calcu la ted assuming G 0 + LE = 150 W uT , T a - 25°C, z Q -0.01 m, a - 0.18, e 8 - 0.97, and e a - 0.85. -55-80h a i , 60 UJ cc < CC £ 40 ui Ul o cc CO 20 AIR TEMPERATURE = 45 °C AIR TEMPERATURE = 25 °C AIR TEMPERATURE =5°C WIND SPEED 6 (m/s) 8 10 FIGURE 2.6 Surface s o i l temperature versus wind speed at a i r temp-eratures of 5, 25, and 45°C. Values were ca l c u l a t e d assuming R s - 1000 W m"2, G + LE » 150 W m~2, z Q -0.01 m, a » 0.18, e 8 = 0.97, and e a - 0.85. -S6' WIND SPEED (m/s) FIGURE 2.7 Surface s o i l temperatures versus wind speed at aero-dynamic roughness lengths of 0.001, 0.01, and_0.1 m. Values were calculated assuming Rg » 1000 W m- , G 0 LE - 150 W nT 2, T a - 25°C, a = 0.18, e 8 = 0.97, and e a » 0.85. -57-clearcut in this study was estimated to be between 0.01 and 0.07 m. At low wind speeds, surface temperature rises sharply in response to decreased roughness. 2.4 CONCLUSIONS By combining energy and radiation balance theory, an equation was derived that enabled surface temperatures to be related to the major energy balance components in a forest clearcut. This expression showed that maxi-mum surface temperatures occur when a l l absorbed radiation is dissipated radiatively (i.e. H + LE + G 0 =0). An increase of 10 Wm in the sum H, LE, or G Q was found to cause a surface temperature reduction of of approximately 1°C. A method was developed for calculating surface temperatures which combined the above expression with aerodynamic theory of sensible heat transfer through the turbulent boundary layer. The aerodynamic theory used to calculate H was shown to require integrated diabatic profile correction functions to correct for the strongly unstable atmospheric conditions in clearcuts. While 3 different approaches tested gave satisfactory results, the integrated diabatic profile method of Paulson (1970) is recommended for estimating fluxes from near surface gradients. The non-integrated diabatic profile methods of Thorn et a l . (1975) and Holbo (1971) were shown to over-estimate the sensible heat flux density under very unstable atmospheric conditions. It should be noted that the analysis in this section is somewhat simplistic in section D since a change in one variable (such as Rs, T a, H, LE, or GQ) is generally accompanied by changes in one or more other -56-variables. These interactive effects w i l l be accounted for in the model described in Chapter 3. -59-CHAPTER 3 MODELLING SOIL TEMPERATURE IN A FOREST CLEARCUT -60-MODELLING SOIL TEMPERATURE IN A FOREST CLEARCUT 3.1 INTRODUCTION The importance of s o i l temperature on the survival and growth of coni-ferous seedlings was discussed in Chapter 1. In Chapter 2 energy balance, radiation balance and aerodynamic theory were used to describe the factors affecting surface s o i l temperatures. In this chapter the above theory is used, in conjunction with a numerical solution of the heat flow equation, to model the time course of s o i l profile temperatures. A physically based model of the s o i l thermal regime has an important application in determin-ing whether high surface s o i l temperatures would occur under certain envi-ronmental conditions and to assess the effects of microclimate modification on seedling root zone temperatures. The availability of a validated s o i l temperature model would also greatly reduce the need for expensive long term programs of s o i l temperature measurement. In the past two decades s o i l temperature modelling has received con-siderable attention in the literature. Fourier analysis was used by Stearns (1967), Ballard (1972), and MacKinnon (1975) to describe s o i l temp-eratures in desert, alpine, and agricultural soils respectively. A f i n i t e difference solution of the heat flow equation was used by Hanks et a l . (1971) and Wierenga and de Witt (1970) to describe s o i l temperatures in agricultural soils; however, as in the previous method i t was necessary to specify the time course of surface s o i l temperature. Gupta et a l . (1981) used the method described by Hanks et a l . and an empirical relationship between air and surface temperature to model temperatures of a s o i l covered with various thicknesses of corn mulch with reasonable success. More -61-recently, Hammel et a l . (1982) developed a physically based method of e s t i -mating the surface s o i l temperature which required the f i n i t e difference solution of the heat and water flow equations for s o i l , and aerodynamic theory (see Chapter 2); however, their midday modelled and measured surface temperatures differed by as much as 20°C. The objectives of this chapter are (1) to develop a physically based model for calculating the time course of surface and subsurface tempera-tures from easily measured or estimated meteorological and s o i l data; (2) to validate the model with measured data from two other energy balance studies reported in the literature as well as s o i l temperatures measured in the forest clearcut described in Chapter 1; and (3) to use a sensitivity analysis approach to show how various environmental factors influence the s o i l thermal regime. 3.2 THEORY The partial differential equation describing the flow of heat in s o i l can be written as follows: - 3G/3z = p s c p s 3T/3t ' (1) where G is the s o i l heat flux density, T i s the s o i l temperature, z is the s o i l depth, t is the time, p s is the s o i l density (solid plus water), and C p S is the specific heat of s o i l . The equation is derived from the heat balance of a unit volume of s o i l where the change in heat storage (p s C p S 3T/3t) is equated with the divergence of the heat flux density (-3G/3z). The s o i l heat flux density is described by Fourier's Law as follows: G = -k s 3T/3z (2) where k s is the apparent thermal conductivity of the s o i l . The negative -61' sign implies that the heat f l u x i s p o s i t i v e (downward) when the temperature gradient i s negative i f z increases with depth. Combining equations (1) and (2) gives: C s 3T/3t = 3(k s 3T/3z)/3z (3) where C s (= p s C p S ) i s the volumetric heat capacity of the s o i l . P h i l i p (1957) suggested a more general heat flow equation incorporat-ing the e f f e c t s of vapour movement on temperature; however, his approach was not used because r e l a t i v e l y l i t t l e heat flows by t h i s process, p a r t i c u -l a r l y when the s o i l i s dry. Kimball et a l . (1972) suggest that heat t r a n s f e r due to vapour movement i s i n s i g n i f i c a n t because the thermal and isothermal vapour fluxes roughly cancel during the day and are very small at night. Non-conductive heat transfer processes can also be accounted for by using the apparent thermal conductivity of the s o i l rather than the true thermal conductivity. The i m p l i c i t f i n i t e difference s o l u t i o n of the heat d i f f u s i o n equation (Richtmeyer 1957) involves the d i s c r e t i z a t i o n of the (z,t) s o i l p r o f i l e plane into a rectangular gri d of points j = 1,2,3... L along the z axis and n = 1,2,3... along the t a x i s . The distance between the v e r t i c a l nodes i s Az and between the time steps i s At. The temperature at any given grid point i s T n j . In f i n i t e d i f f e r e n c e form, (3) becomes: n n-1 T - T At 1_ Az V/1 n-1 n n-1 n 1 / T + T - T - T 2Az \ J+l i+1 J 1 - k, V/2 (n-1 n T + T J J n-1 n - T - T (4) -63-for a l l interior nodes j = 2 to L - l . When the upper and lower boundary conditions T n and Tn, are specified a system of L equations in L unknowns is established. These equations form a tridiagonal matrix which can be solved quite efficiently using Gaussian elimination. The derivation and method of solving (4) are outlined in more detail in Appendix A. The advantages of using a f i n i t e difference approximation rather than an analytical solution of the heat flow equation are that (i) the s o i l thermal properties can vary with depth and in time and ( i i ) the upper and lower boundary conditions can be arbitrary rather than periodic functions of time. The major d i f f i c u l t y with this type of model involves the specification of the upper and lower boundary conditions. The lower boundary condition is quite easily simulated by modelling the s o i l profile to a depth where the temperature does not change appreciably with time. Over a period of a few weeks this depth is about 1 m. On an annual basis a profile depth of a few meters would be required. The surface boundary condition is much more d i f f i c u l t to estimate because i t depends upon many factors. Using the radiation and energy balance theory outlined in Chapter 2 the s o i l surface temperature can be calculated as follows: T s = ( R s ( l - a) + e s RL - H - LE - GQ) / e s a 0. 25 (5) By (i) estimating H by the method described in Chapter 2, ( i i ) approximating LE using a simple empirical method (described later), and ( i i i ) calculating GQ using (2) in f i n i t e difference form after solving (4), and (iv) knowing the solar irradiance, air temperature, windspeed, albedo, and atmospheric and surface emissivities, i t is possible to -64-iteratively calculate the surface temperature. 3.3 PROCEDURE The s o i l temperature simulations were run on an 8 bit microcomputer (Campbell Scientific Inc., Logan, Utah, U.S.A., Model C2000 with a 4 MHz Z80 CPU) having only 64 K bytes of random access memory. The model was written in the basic programming language (Microsoft Mbasic) and run in a compiled form to increase the operating efficiency. The microcomputer con-tained just enough memory to simulate a s o i l profile depth of 0.8 m. It took between 10 to 15 minutes of computer time to simulate one actual day using 24 one hour time steps. The model could be easily transferred to a mainframe computer where the execution time would increase by several orders of magnitude and memory limitations would not determine the depth of the s o i l profile. The program source code is listed in Appendix D. A. Assumptions Used i n the Model Three assumptions were made to simplify the solution of the surface so i l temperature. The f i r s t was that the latent heat flux density could be modelled empirically using: LE = c Rs (6) where c i s an experimentally determined coefficient. This relationship i s supported by the measurements of LE in a forest clearcut by Ballard et a l . (1975). During a period of summer drought the latent heat flux density was approximately in phase with the solar irradiance during the day and was negligible at night. Their data suggested that c « 0.05 under dry surface so i l conditions. Under moist surface s o i l conditions, the potential -65-evaporation formula of Jensen and Haise (1963) (LE = 0.025°C (T a + 3°C)RS) could be used to give a more accurate estimate of LE. The empirical modelling of LE is a shortcoming in the model at the present time; however, the model is set up in such a way that this problem can be overcome. A numerical s o i l water flow model running concurrently could be used to calculate the vapour pressure of the air at the so i l sur-face (Hammel et a l . 1982). This would allow an aerodynamic estimate of the latent heat flux density to be made using the diabatic profile correction factors calculated for sensible heat transfer and an additional iteration routine. The second assumption was that the semi-empirical diabatic influence functions described in Chapter 2 could be used to correct the aerodynamic estimate of the sensible heat flux density under non-neutral conditions on a sloping forest clearcut. Although the diabatic influence functions are semi-empirical, relationships, which have been derived from many large data sets obtained over generally horizontal land surfaces and dimensional analysis, they quite accurately correct for the effects of buoyancy on a sloping clearcut as shown in Chapter 2. A third assumption was that equations (4) and (5) in Chapter 2 could adequately estimate the atmospheric emissivity under cloudy skies. Very reasonable estimates of e a were obtained using this method and even a relatively large margin of error (eg. ± 0.05) would not seriously affect the surface energy balance. B. The Solution Method The solution involves the sequential estimation of s o i l temperatures -66-at time (t = 2) based upon the temperature distribution in the s o i l profile at time (t = 1) and estimates of the upper and lower boundary temperatures at time (t =2). The model is initiated by representing the s o i l profile as a series of nodes in a two-dimensional array; one dimension for depth and the other for time. Each of the nodes at time (t = 1) is given an i n i t i a l temperature (uniform through the profile), thermal conductivity and volumetric heat capacity. Values of the thermal conductivity and volumetric heat capacity are also specified at time (t = 2). Solar irradiance, air temperature, wind speed, and % sky cloud cover are then read from a data f i l e for time (t = 2). The latent heat flux density and atmospheric emissivity are semi-empirically determined from these variables. Rough estimates of the surface temperature and s o i l heat flux density are also made. Using an iterative procedure, convergence between the aerodynamic and energy balance estimates of the sensible heat flux density i s found at a unique surface temperature. This temperature is used as the upper boundary condition in the s o i l temperature model at time (t = 2). Soil temperatures between the upper and lower boundaries and the near surface s o i l heat flux density (using the f i n i t e difference form of (2)) are then calculated. If the estimated, (i.e. the value used in (5)), and calculated s o i l heat flux densities are not equal, a new estimate of the s o i l heat flux density i s determined and the process is repeated. When the estimated and calculated s o i l heat flux densities converge, the correct s o i l temperature distribu-tion is attained. Pertinent data is then stored or printed, the s o i l temperature d i s t r i -bution at the next time step is set equal to that at the present time step, new meteorological data is read from a f i l e , and the process is repeated. The flow chart in Figure 3.1 shows this solution method. 3.4 RESULTS AND DISCUSSION A. Validation of the Model Extensive testing was undertaken to determine the validity of the model. The finite difference s o i l heat flow component was tested f i r s t by comparison with the analytical solution for heat flow through a i m thick slab having uniform thermal properties. Using a sinusoidal upper boundary condition, very good agreement (T = ±0.1°C) was found between the two solution methods over the range of thermal properties and surface tempera-ture amplitudes observed in so i l s . The fi n i t e difference s o i l heat flow component of the model was evaluated using Fourier analyzed s o i l temperatures reported by Stearns (1967). These data were collected on July 11, 1964 in a dry sandy Peruvian desert s o i l having approximately constant thermal properties through the s o i l profile. In this test, the temperature measured at the 1 mm depth by Stearns was used as the upper boundary in the f i n i t e difference solution. Modelled and measured temperatures at the 0.02, 0.1, and 0.2 m s o i l depths are compared in Figure 3.2. Excellent agreement was found through the whole diurnal period at each of the depths considering (1) measurement errors of ±0.5°C, (2) errors in the estimated thermal properties and in the so i l homogeneity assumption, (3) errors in the temperature probe location (±1 cm), and (4) errors associated with the Fourier transcription of Steam's data (±0.5°C). ~6B-INITIALIZE MODEL T ( r . l ) , ks(z , t ) , C B ( i , t ) 1 INPUT Ra, T a and u at t-2 I CALCULATE p a , c p a , LE, and c a ESTIMATE G 0 ( e e t ) , T B CALCULATE R L, eB<TT8<>, Rn, and H (energy balance) Estimate New G(est) CALC T(z,2) and G 0 Set T ( r , l ) - T(t ,2 ) END FIGURE 3.1 Flow chart showing the so lu t ion method used to ca lcu la te surface and subsurface s o i l temperatures. -69-i 8 12 TIME (hrs) FIGURE 3.2 Comparison of modelled and measured (Stearns 1967) temp-eratures i n a dry Peruvian desert s o i l having uniform thermal properties with depth on Ju ly 11, 1964. -70' These two tests suggest that the s o i l heat flow component of the model is very accurate. It was assumed that the fi n i t e difference model worked equally well for a s o i l profile having thermal properties that varied with depth. The complete model, using the energy balance and aerodynamic theory to estimate the surface s o i l temperature, was tested with measured energy fluxes and s o i l temperatures reported by Novak (1981) for a large agricul-tural f i e l d at Agassiz, B.C. Tests were conducted on two diurnal periods; one on May 18, 1979, a clear day when the s o i l water content was near f i e l d capacity, and the other on July 21, 1979, another clear day after the s o i l had dried considerably near the surface. Novak (personal communication) supplied hourly values of the solar irradiance,°air temperature, wind speed, and latent heat flux density in addition to the s o i l thermal proper-ties, albedo, surface emissivity, and temperature at the bottom of the modelled s o i l profile (0.8 m). The surface roughness length was estimated to be 0.005 m using values quoted in the literature for a smooth agricul-tural f i e l d . The i n i t i a l profile temperature distribution was obtained by running the model for one day, starting from a temperature profile set equal to the lower boundary condition. This procedure gave an i n i t i a l temperature distribution very similar to that measured by Novak. The comparison of modelled and measured energy fluxes and s o i l temp-eratures is shown in Figure 3.3 for the May 18 test. In general, excellent agreement was found between modelled and measured values except for a few hours near midday. During this time the modelled sensible heat flux was overestimated by up to 30 W m and as a result the modelled surface temperature was underestimated by about 2.5°C. Modelled and measured -71' 1000 T r - i — i — i — - T — i — r m i l «. • IS M * — «. 0 -200 55 • • i I U J L J I L. 45 r-T — i — r SOIL DEPTH MODELLED MEASURED • 0.0 m O 0.025 m -0.1 m • — — • • 0.2 m • — — 0.5 m A — 1.0 m A? ^J.I--U-i.'.'-'.^L 12 18 2 4 T I M E ( HRS. RS.T. ) FIGURE 3.3 Comparison of modelled and measured (Novak 1981) tempera-tures and surface energy fluxes on May 18, 1979 i n a wet s i l t loam s o i l at Agassiz , B .C. -7Z' temperatures at the 0.1, 0.2 and 0.5 m depths were in near perfect agree-ment. There was a discrepancy of about 3°C between the midday modelled and measured temperatures at the 0.025 m depth. This was partially due to the error in the modelled surface temperature but i t also could have been due to an overestimate of the volumetric heat capacity in this layer. The results of the July 21 test are shown in Figure 3.4. Excellent agreement between modelled and measured energy fluxes was obtained through-out the entire day. Modelled and measured surface temperature were also in close agreement during the daytime; however, modelled values were up to 3°C less than measured values in the early morning. Subsoil temperatures at the 0.1, 0.2, and 0.5 m depths were a l l modelled to within less than 1°C through the whole day. The temperature at the 0.025 m depth was underesti-mated during the daytime and overestimated during the nighttime, suggesting as in the May 18 test, that the volumetric heat capacity just below the surface was incorrectly estimated. The results shown in Figures 3.3 and 3.4 suggest that the model is capable of accurately simulating the time course of energy fluxes, surface and subsurface temperatures from very few data inputs. The relatively small differences between modelled and measured s o i l temperatures are probably less than the spatial variability that could be expected within a few meters of the temperature sensors. A l l of the variables used as model inputs except the latent heat flux density are routine climate station measurements. The objective of the next test was to determine i f the simple empirical LE model (6) could be used to model LE with sufficient accuracy to model the surface tempera-ture. Theory outlined in Chapter 2 suggests that a change in the latent -73-- 2 0 0 55 45 35 ec < BC z i i ! 25 o CO 15 o' 1 / A / * // // 9% 0 / / A L _ f as. / -A A \o< V • ' o v N ° pTTr _ _ _ _ • o*.-SOIL DEPTH MODELLED MEASURED ao m Ojoes m 0.1 m 0.2 m 0.9 m 1.0 m O A o • • t I 1 1 1 L 6 12 18 24 TIME (HR& R&T.) FIGURE 3.4 Comparison of modelled and measured (Novak 1981) tempera-tures and surface energy fluxes on July 21, 1979 In a s l i t loam s o i l at Agassiz , B . C . which had dr ied near the surface. -74-heat flux density of 10 W m- would cause roughly a 1°C surface temperature change. It should be possible then, to use an empirical LE model to calcu-late the surface temperature to within at least a few degrees. The test of the simple LE relationship for calculating s o i l tempera-tures from standard climate data was conducted on the clearcut site de-scribed in Chapter 1 for a period of the six hottest clear sunny days in August 1981. The so i l profile had dried considerably from the surface to a depth of well over 0.3 m. To account for this effect on the s o i l thermal properties in the model, the s o i l was divided into 4 layers. The thermal conductivity and volumetric heat capacity in the 0-0.02, 0.02-0.05, 0.05-0.25 and 0.25-0.8 m layers were 0.3, 0.3, 0.5, 0.7 W m _ l oC _ 1 and 1.2, 1.6, 2.0, 2.2 MJ m - 3 oC - 1 respectively. The temperature measured at the 0.8 m depth was 15.0°C and did not change significantly during the six days. This was used as the lower boundary condition. A roughness length of 0.05 m was used to characterize the surface. Measured solar irradiance, air temperature, and wind speed were used as inputs to the model. Since the s o i l surface was dry, c was assumed to be 0.05 in (6). The surface emissivity and albedo used were 0.95 and 0.16 respectively. The results of this simulation are shown in Figure 3.5. Modelled and measured temperatures at the surface and the 0.15, 0.3, and 0.5 m s o i l depths generally showed agreement to within less than 2°C. The increase in mean temperature with time was also followed well by the model. The model-led amplitude at the 0.15 m depth was underestimated by approximately 2°C; however, the modelled mean daily temperature at this depth showed good agreement with the measured values. The modelled amplitude at the 0.3 m -75-u e LU or z> »-< or LU LU O CO 80 60 40 20 0 40 30 20 10 0 40 30 20 10 0 40 30 20 10 SOIL SURFACE r * T r 0.15 m DEPTH "St. «" -i 1 1 1 r 0.3 m DEPTH T 1 1 1 0.5 m DEPTH 24 48 72 96 TIME ( hrs RS.T.) 6 7 8 9 DATE (AUGUST, 1981) 120 144 10 FIGURE 3.5 Comparison of modelled (0) and measured (+) s o i l tempera-tures at 4 depths in the Cameron Valley forest clear-cut experimental site during the 6 hottest days in 1981. 76-depth were slightly greater than the measured values, but again the mean daily temperature was well modelled. The discrepancy between modelled and measured amplitudes at these depths suggests that the estimated s o i l thermal properties were slightly in error. These results show that a simple function describing the latent heat flux density can be incorporated in the model so that easily measured climate station data can be used to simulate s o i l temperatures under dry s o i l conditions. Another objective was to test the model for a longer simulation period and under variable weather conditions. A sixteen day interval beginning on July 2, 1981 was chosen for this test. The weather for the week preceding July 2 was warm and dry. From July 2 to the July 5 i t was hot and dry. In the early morning of July 6 a frontal storm approached and cloudy cool weather occurred until the July 11. Less than 10 mm of rain f e l l during this period. From July 12 un t i l the July 17 a gradual warming trend occurred and no rain f e l l . As in the last test, four layers were used to represent the change in thermal properties with s o i l depth. The thermal conductivity and volumetric heat capacity in the modelled 0-0.02, 0.02-0.05, 0.05-0.25, and 0.25-0.8 m layers were 0.19, 0.23, 0.76, 0.96 W m _ 1C _ 1 — 3 — 1 and 0.78, 0.95, 1.4, 1.6 MJ m C respectively. The temperature measured at the 0.8 m depth was 8.0°C and the surface roughness, albedo, and surface emissivity remained the same as in the previous test. Measured solar irradiance, air temperature, and wind speed were used as model inputs. The latent heat flux density was estimated to be 5 percent of the solar i r r a d i -ance, as in the previous test. The cloudy sky atmospheric emissivity (equations (4) and (5), Chapter 2) function was used throughout the -77-interval and i t was assumed that the sky was completely cloudy from July 5 until July 11. The results of the test are shown in Figure 3.6. Modelled and measur-ed s o i l temperatures at the surface, 0.1, and 0.2 m s o i l depths showed very good agreement. Modelled surface temperature showed an accurate response to the change in the weather. Modelled subsoil temperatures also reflected the pattern of heat storage, heat loss, then heat storage during the con-secutive clear hot, cloudy cool, then clear hot periods. Modelled and measured temperatures differed at the 0.1 and 0.2 m depths for the f i r s t two days. This reflected inaccuracies in the i n i t i a l temperature distribu-tion in the profile rather than a deficiency in the model. The results shown in Figure 3.6 further suggest that standard climate data can be used to model s o i l surface and subsurface temperatures with very good accuracy. On days when the latent heat flux density is higher than in the two cases modelled here either an increase in the c values used in (6) or the method of Jensen and Haise (1963) could be used to more accurately simulate the course of LE. More research is necessary to determine the effect of s o i l moisture content on the relationship between LE and Rs in clearcuts. B. Sensitivity Analysis of the Model The effects of the aerodynamic roughness of the surface, the latent heat flux density, the s o i l thermal properties, and the temperature at the bottom of the s o i l profile on modelled s o i l profile temperatures were determined using a sensitivity analysis. Changes in the s o i l thermal regime due to microclimate modification and the necessary accuracy of e s t i -mated site dependent variables can be determined from these tests. -76-60 40 — 20 O -i r SURFACE < * o W 30 Q_ _1 I u _ l I 1 L . -i 1 1 r 0.2 m depth _ i i i 1 1-6 8 10 12 14 TIME (days ) JULY 1981 16 FIGURE 3.6 Comparison of modelled (0) and measured (+) soil tempera-tures at 3 depths in the Cameron Valley forest clear-cut experimental site during 16 days of variable weather conditions in July 1981. -19-A six day test period was used so that both cumulative and immediate effects could be observed. The test was conducted by selecting a set of standard conditions and allowing the variable of interest to change while holding a l l others constant. The standard conditions used in a l l of the tests (and represented by the (x) symbol in Figures 3.7-3.10) were as f o l -lows: measured solar irradiance, air temperature, and wind speed from August 5-10, 1981, c = 0.05, z Q = 0.005 m, T(0.8 m) = 12.0°C, a = 0.16, e s = 0.97 and the volumetric heat capacity and thermal conductivity of the 0.0-0.02, 0.02-0.05, 0.05-0.8 m s o i l layers were estimated to be 1.0, 1.2, 1.4 MJ m-3^"1, and 0.19, 0.29, 0.67 W m - 1 0^ 1 respectively. The effect of the aerodynamic roughness length on s o i l profile temp-eratures is shown in Figures 3.7a-d. Increasing z Q from 5 to 10 mm reduced midday surface s o i l temperatures by about 5°C. Increasing z Q from 5 to 50 mm reduced midday surface temperatures by about 15°C. Very l i t t l e effect was noticed at the surface at night because z Q influences sensible and latent heat transfer and these fluxes are small at night. The aerodynamic roughness had a large influence on so i l warming at depth. Decreasing the roughness length caused less of the available energy at the surface to be partitioned as sensible heat, and as a result more so i l warm-ing occurred. This effect was particularly noticeable at the 0.3 and 0.5 m depths as shown in Figures 3.7c and 3.7d. Hammel et a l . (1982) also noted the significance of surface roughness and suggested that i t was the most important factor involved in the effect of tillage on s o i l temperature. Management practices which increase the aerodynamic roughness of the surface would be beneficial on forest clear-cuts prone to high surface s o i l temperatures. -80-The effect of the latent heat flux density upon s o i l profile tempera-tures is shown in Figures 3.8a-d. Soil temperatures were calculated using c = 0.0, 0.05 and 0.2 in (6) to give a wide range of s o i l evaporation rates. Surface s o i l temperatures decreased as the latent heat flux increased. Midday surface s o i l temperatures were 10°C lower when LE was _ o increased from 0 to about 200 W m under these conditions. This appears to disagree with the observation in Chapter 2 that a 100 W m~ decrease in H + LE + G causes a 10°C increase in T s. Both observations are correct. The simulation in this chapter shows that while LE is decreased by 200 W m-2, H + G Q has increased by about 100 W m~2. Novak (1981) observ-ed this effect in his energy balance-soil temperature study of a bare s o i l . Increasing the latent heat flux density also had the effect of reducing the total amount of available energy at the surface for s o i l heat transfer so that subsurface s o i l temperatures did not increase as rapidly. The effect on s o i l profile temperatures of doubling the volumetric heat capacity or doubling both the volumetric heat capacity and the thermal conductivity are shown in Figures 3.9a-d. Roughly, doubling C s and k s would be the result of increasing the volumetric water content from about 10 to 20% in a coarse textured mineral s o i l . These changes in thermal properties had significant effects upon s o i l temperatures (up to 10°C). Figure 3.9 shows that increased moisture content at the surface decreases T s mainly because of i t s effect of increasing k s rather than C s assuming LE is constant. Increasing the thermal conductivity increased both the rate of s o i l warming and the amplitude of the temperature wave at depth. While increasing the volumetric heat capacity slightly reduced surface s o i l temperatures i t , markedly decreased temperature at depth. At T o 60 UJ CC H 401 cc UJ OL M W M o o o o oo • • o o o o U J o CO 20 "8 0 o - o j * Ml SOIL SURFACE 7 8 9 DATE (days) AUGUST 1981 10 FIGURE 3.7a The effect of the aerodynamic roughness (z Q) on modelled surface s o i l temperature. Values of z Q used were 5 mm (x) , 10 mm (+), and 50 mm (0) while a l l other variables were held constant. 0 6 0 O L U or H 4 0 or L U o. O to 0.1 m DEPTH «:::« 5 0 °°Sssgss° • " * * • . « o°°°o0*?« 7 8 DATE (days) AUGUST 1981 10 FIGURE 3.7b As f o r Figure 3.7a but for temperatures at the 0.1 m depth. T E M P E R A T U R E > 01 o 0Q c 1-1 a OJ » a* c rt rti O H § (D » rt C M (0 CD » rt n> o S a. n T3 -£0-TEMPERATURE 5 6 7 8 DATE (days) AUGUST 1981 FIGURE 3.8a The effect of the latent heat flux density (LE) modelled surface s o i l temperature. Values of LE were 0.0 (+), 0.05 Rg (x), and 0.2 Rj, (0) while other variables were held constant. 7 8 9 DATE (days) AUGUST 1981 FIGURE 3.8b As f o r Figure 3.8a but f o r temperatures at the 0.1 m depth. 0 60 © U J cc 3 £ 40 cc U J Q . UJ o CO 20 gggSoSoSSSooSsgsaigScS^ gjSSoOoooooogggggJgggoooO )88«88i 6 7 8 9 DATE (days) AUGUST 1981 FIGURE 3.8c As f o r Figure 3.8b but f o r temperatures at the 0.3 m depth. DATE (days) AUGUST 1981 FIGURE 3.8d As f o r Figure 3.8b but for temperatures at the 0.5 m depth. 0 6 0 o 111 CC H 4 0 cr UJ CL £20 o • • **«»** » «o **»* I oa *M 0 o o2 ! :8s SOIL SURFACE 7 8 DATE (days) AUGUST 1981 FIGURE 3.9a The effects of s o i l thermal properties on modeled sur-face s o i l temperature. Values of k g and C s i n the 0-0.02, 0.02-0.05 and 0.05-0.8 m s o i l layers were 1.0, 1.2 and 1.4 raj m - 3 ° C - 1 and 0.19 0.29 and 0.67 W m - 1 °C~ l respect ively (X) , values of k 8 are the same as those given above but values of C s are twice those given above (+) and values of k s and C g are twice those given above (0 ) . A l l other variables were held constant. o 60 O UJ cc 3 < CC U l CL 2 U J o C O 40 ***** ***** gooogg 20 8 8 7 8 9 DATE (days) AUGUST 1981 10 FIGURE 3.9b As f o r Figure 3.9a but for temperatures at the 0.1 m depth. 0.3 m DEPTH O60 O L U cr cr O CO I N O « 6 7 8 9 DATE (days) AUGUST 1981 FIGURE 3.9c As f o r Figure 3.9a but for temperatures at the 0.3 m depth. 0 60 O UJ cc ? 40 cc UJ o. I 20 o CO 0.5 m DEPTH 7 8 9 DATE (days) AUGUST 1981 10 FIGURE 3.9d As f o r Figure 3.9a but for temperatures at the 0.5 m depth. o 60 o • UJ CC 3 I -< cc U l Q_ 2 UJ o CO 40 20 7 8 9 DATE (days) AUGUST 1981 SOIL SURFACE i 10 VTCURE 3 10a The effect of the temperature at the lower surface of FIGURE 3.10a The ^ ^ ^ ^ d e p t h ) [ l o w e r b o u n d condition] on modelled surface s o i l temperature. Values of 10°C (+), 12°C (x), and 14°C (0) were used while a l l other variables were held constant. o 60 O 40 tr I-< tr LU Q. 2 LU _ _ {_ 20 o CO r *0 Q *o «99 9 ' t o ' 8 8 W 6 7 8 DATE (days) AUGUST 1981 FIGURE 3.10b As for Figure 3.10a but for temperatures at the 0.1 m depth, NO DATE (days) AUGUST 1981 FIGURE 3.10c As for Figure 3.10a but for temperatures at the 0.3 m depth. DATE (days) AUGUST 1981 FIGURE 3.10d As for Figure 3.10a but for temperatures at the -97-depth there is a compensating effect of increasing both k s and C s, since increasing k s w i l l tend to increase the mean temperature and increasing C s w i l l tend to decrease i t . Thermal properties obviously have a large influence upon s o i l thermal regimes. On sites which are prone to high surface temperatures, management practices, such as scalping off the surface organic layer or mixing the surface organic layer with the mineral s o i l below which tend to increase C s or k s w i l l reduce lethally high daytime surface temperatures. The f i n a l sensitivity test involved the effect of the lower boundary s o i l temperature upon temperatures at other depths in the s o i l profile. The test was conducted because of the d i f f i c u l t y in obtaining measured deep so i l temperatures, particularly on remote sites where estimates may be in error by a few degrees. The test was conducted by using 10, 12, and 14°C as the lower boundary temperatures and noting the effects at other depths. The results of the test are shown in Figures 3.10a-d. They show that the modelled lower boundary temperature has no effect upon the surface tempera-ture and very l i t t l e effect upon root zone s o i l temperatures. Provided that the lower boundary is deep enough in the s o i l , good estimates of s o i l temperature can be made even i f the lower boundary temperature is slightly in error. 3.5 CONCLUSIONS A f i n i t e difference numerical s o i l temperature model was developed to calculate surface and subsurface temperatures in soils having thermal pro-perties which vary with depth. The model accurately predicted s o i l temp-eratures in a forest clearcut over extended periods of time using easily '98' measured climate station data. Excellent agreement was also found using data from two other energy balance studies reported in the literature. A sensitivity analysis showed that the s o i l thermal properties, the latent heat flux density, and the aerodynamic roughness of the surface and the lower boundary s o i l temperature could be estimated with sufficient accuracy from relationships in the literature when measurements of these variables were not available from the site. -99-CONCLUDING REMARKS Shade tolerant western hemlock and Pacific silver f i r were not well adapted to the severe heat and water stress which occurred during the summer on a high elevation south facing slope. Douglas-fir, a moderately shade intolerant species, survived much better under these conditions. Modifying the seedling microclimate with shade cards, irrigation, or a combination of these two treatments significantly increased survival. Presumably, these treatments w i l l be more effective during hot dry summers when there is inadequate water in the s o i l to supply plant requirements and suppress surface s o i l and seedling root collar temperatures. The cost effectiveness of microclimate modification should be carefully examined for different types of sites. The use of good stock, careful planting, and attention to micro-site conditions may also substantially increase survival on stress prone sites. Using the theory developed in Chapter 2 i t is possible to speculate about alternative microclimate modification techniques. High surface s o i l and root collar temperatures can be decreased by either reducing the amount of solar radiation absorbed by the surface or dissipating the absorbed energy more efficiently as sensible, latent, or s o i l heat. Blocking the solar irradiance is probably the most effective means of reducing the energy input. This can be accomplished to a lesser degree by modifying the solar re f l e c t i v i t y of the surface. Evaporative cooling of the surface can be substantially increased by irrigating dry s o i l ; however, this is not a very practical operational procedure. The sensible heat flux density is strongly influenced by the aerodynamic roughness of the surface wo-and the wind speed. Designing cut blocks to maintain air flow and keeping the surface as rough as possible w i l l help reduce high surface tempera-tures. The seedling thermal regime can also be improved by planting in mineral s o i l rather than s u r f i c i a l organic material. Scarifying increases both the thermal conductivity and volumetric heat capacity of the s o i l and as a result surface temperature fluctuatons are reduced and s o i l tempera-tures at depth respond more rapidly to conditions at the surface. The procedures for increasing the sensible and s o i l heat flux densities are probably the most inexpensive and easily applied methods of improving the s o i l thermal regime; however, these measures may not be effective enough on some sites. The computer model described in Chapter 3 can be used to predict the effects of changing various site characteristics and weather conditions on s o i l temperatures. Standard climate data was used to successfully model s o i l temperatures for extended periods of time in the clearcut. With minimal additional work the model could be used to model s o i l temperatures at remote sites based upon existing climate and soils data. The model could also be used to provide a s o i l temperature growing degree day index and could provide the means for a probablistic analysis of planting success. Other practical applications w i l l become apparent as the model receives more use. -707-R E F E R E N C E S Arnott, J.T. 1975. Field performance of container-grown and bare root trees in coastal British Columbia. Can. J. For. Res. 5: 186-194. Barker, E.H. and T.L. Baxter. 1975. A note on the computation of surface layer fluxes for use in numerical modelling. J. Appl. Meteorol. 14: 620-622. Baker, F.S. 1929. Effect of excessively high temperatures on coniferous reproduction. J. For. 27: 949-975. Ballard, T.M. 1972. Subalpine s o i l temperature regimes in southwestern British Columbia. Arctic and Alpine Research 4: 139-146. Ballard, T.M., T.A. Black and K.G. McNaughton. 1977. Summer energy balance and temperatures in a forest clearcut in south-western British Columbia. In: Energy, Water, and the Physical Environment of the S o i l . Report of the 6th B.C. Soil Science Workshop, Richmond, B.C. April 20-21, 1977. Bercowicz, R. and L.P. Prahm. 1982. Sensible heat flux estimated from routine meteorological data by the resistance method. J. Appl. Meteorol. 21: 1845-1864. Brix, H. 1967. An analysis of dry matter production of Douglas-fir seed-lings in relation to temperature and light intensity. Can. Jour. Bot. 45: 2063-2072. Businger, J.A. 1975. Aerodynamics of Vegetated Surfaces. In: Heat and Mass Transfer in the Biosphere, Part 1 - Transfer Processes in the Plant Environment, (Devries, D.A. and N.H. Afgan, eds.), Scripta Book Co., John Wiley and Sons, New York. pp. 5-28. Businger, J.A., J.C. Wyngard, Y. Izume and E.F. Bradley. 1971. Fulx-pro-f i l e relationships in the atmospheric surface layer. J. Atmos. Sci. 28: 181-189. Businger, J.A. 1966. Transfer of momentum and heat in the planetary boundary layer. Proceedings of the Symposium on Arctic Heat Budget and Atmospheric Circulation (The Rand Corporation), pp. 305-332. Campbell, G.S. 1977. An Introduction to Environmental Biophysics, Springer-Verlag, New York. Carslaw, H.S. and J.C. Jaeger. 1959. Conduction of Heat in Solids, 2nd Ed., Oxford University Press, London. Cleary, B.D., R.D. Graves and R.K. Hermann. 1978. Regenerating Oregon's Forests. Oregon State University Extension Service, Corvallis, Oregon. - 7 0 2 -Cochran, P.H. 1969. Thermal properties and surface temperatures of seed-beds. Pacific Northwest Forest and Range Expt. St., USDA Forest Service, Portland, Oregon. De Vries, D.A. 1963. Thermal Properties of Soils. In: Physics of Plant Environment, (Van Wijk, W.R., ed.), North Holland Publishing Co., Amsterdam, The Netherlands. Dobbs, R.C. and R.G. McMinn. 1977. Effects of scalping on s o i l tempera-ture and growth of white spruce seedlings. In: Energy, Water and the Physical Environment of the Soil. Report of the 6th B.C. Soil Science Workshop, Richmond, B.C., April 20-21, 1977. Dyer, A.J. and B.B. Hicks. 1970. Flux-gradient relationships in the constant flux layer. Quart. J.R. Met. Soc. 96: 715-721. Dyer, A.J. 1967. The turbulent transfer of heat and water vapour in an unstable atmosphere. Quart. J. Roy. Meteor. Soc. 93: 501-508. Freeze, R.A. and J.A. Cherry. 1979. Groundwater, Prentice-Hall, Inc., Englewood C l i f f s , New Jersey, p. 544-545. Fuchs, M. and CB. Tanner. 1968. Calibration and f i e l d test of s o i l heat flux plates. Soil Sci. Soc. Amer. Proc. 32: 326-328. Fuchs, M., CB. Tanner, G.W. Thurtell and T.A. Black. 1968. Evaporation from Drying Surfaces Using the Combination Method. Agronomy Journal 61: 22-26. Geiger, R. 1965. The Climate Near the Ground. Harvard Univ. Press, Cambridge, Mass. Gupta, S.C, J.K. Radke and W.E. Larson. 1981. Predicting temperatures of bare and residue covered soils with and without a corn crop. Soil Sci. Soc. Am. J. 45: 405-412. Hammel, J.E., R.I. Papendick and G.S. Campbell. 1981. Fallow Tillage Effects on Evaporation and Seedzone Water Content in a Dry Summer Climate. Soil Sci. Soc. Am. J. 45: 1016-1022. Hanks, R.J., D.D. Austin and W.T. Ondrechen. 1971. Soil temperature e s t i -mation by a numerical method. Soil Sci. Soc. Amer. Proc. 35: 665-667. Heninger, R.L. and D.P. White. 1974. Tree seedling growth at different s o i l temperatures. For. Sci. 20: 363-367. Hermann, R.K. 1964. Paper mulch for reforestation in Southwest Oregon. J. For. 62: 98-101. Hermann, R.K. and W.W. Chilcote. 1965. Effect of seedbeds on germination and survival of Douglas-fir. For. Res. Lab. Oregon State University, Corvallis, Oregon, Research Paper 4, 28 p. -703-Hobbs, S.D., R.H. Byars, D.C. Henneman and CR. Frost. 1980. First year performance of 1-0 containerized Douglas-fir seedlings on droughty sites in southwest Oregon. Forest Research Laboratory. Oregon State University, Corvallis, Oregon. Holbo, H.R. 1971. Development of a stability correction for estimating Corrective Transfer. Technical Report No. 71-7, Department of Atmospheric Services, Oregon State University, Corvallis, Oregon. Idso,' S.B. and R.D. Jackson. 1969. Thermal radiation from the atmos-phere. J. Geophys. Res. 74: 5397-5403. Jackson, R.D. and S.A. Taylor. 1965. Heat Transfer. In: Methods of Soil Analysis, Part 1 - Physical and Mineralogical Properties Including Statis-tics of Measurement and Sampling. C.A. Black, ed. American Society of Agronomy, Inc., Madison, Wisconsin, pp. 349-360. Jensen, M.E. and H.R. Haise. 1963. Estimating evapotranspiration from solar radiation. J. I r r i g . Drain. Div., Am. Soc. Engr. 89: 15-41. Kimball, B.A., R.D. Jackson, R.J. Reginato, F.S. Nakayama and S.B. Idso. 1976. Comparison of f i e l d measured and calculated s o i l heat fluxes. Soil Sci. Soc. Amer. J. 40: 18-25. Klinka, K., F.C Nuszdorfer and L. Skoda. 1979. Biogeoclimatic units of central and southern Vancouver Island. Prov. of British Columbia, Min. of Forests. 120 pp (1 map). Korstian, C F . and N.J. Fetherolf. 1921. Control of stem girdle of spruce transplants caused by excessive heat. Phytopathology. 2: 485-490. Lavender, D.P. and W.S. Overton. 1972. Thermoperiods and s o i l tempera-tures as they affect growth and dormancy of Douglas-fir seedlings of dif-ferent geographic origin. For. Res. Lab. Sch. For., Oregon State Univer-sity, Corvallis, Oregon, Research Paper 13, 26 p. Leckie, D.C, T.A. Black and P.A. Murtha. 1981. Development and testing of a method of estimating sensible heat flux from natural surfaces using remotely sensed surface temperatures. Can. Jour, of Chem. Eng. 59: 189-194. Lettau, H.H. 1962. Notes on theoretical models of profile structure in the diabatic influence layer. Final Rep. Contract DA-36-039-SC-80282. USEPG, Ft. Huachuca, Arizona, and Dept. of Meteor., University of Wisconsin, Madison, Wisconson. pp. 195-226. L i s t , R.J. 1971. Smithsonian Meteorological Tables, 6th Ed., Washington, D.C, Smithsonian Institution Press. Livingston, N.J. and R.J. Stathers. 1982. Annual Report for Seedling Establishment Project #316-1. MacMillan Bloedel Ltd., Inhouse Report, 130 pp. -/Of-Livingston, N.J. and R.J. Stathers. 1983. Annual Report for Seedling Establishment Project #316-1. MacMillan Bloedel Ltd., Inhouse Report, 35 pp. Livingston, N.J., T.A. Black, D. Beames and B.G. Dunsworth. An instrument for measuring the average stomatal conductance of conifer seedlings. Can. J. For. Res. (submitted). Lumley, J.L. and H.A. Panofsky. 1964. The Structure of Atmospheric Turbu-lence. Interscience Monographs and Texts in Physics and Astronomy, Vol. XII. Wiley, New York. MacKinnon, J.C. 1975. Modelling the thermal regime of a Nova Scotian s o i l under oats. Can. Agric. Eng. 18: 41-45. Maguire, W.P. 1955. Radiation, surface temperature, and seedling survi-val. For. Sci. 1(4): 277-285. Minore, D. 1970. Shade benefits Douglas-fir in Southwestern Oregon cutover area. U.S.D.A. Forest Service Research Note. 2 p. Pacific Northwest Forest and Range Exp. St., Corvallis, Oregon. Monin, A.S. and A.M. Obukhov. 1954. Dimensionless characteristics of turbulence in the surface layer. Akad. Nauk. SSSR. Geofiz. Inst., Tr. 24: 163-187. Monteith, J.L. 1973. Principles of Environmental Physics. Edward Arnold, London. Munn, R.E. 1966. Descriptive Micrometeorology. Academic Press, New York. Novak, M.D. 1981. The moisture and thermal regimes of a bare s o i l in the Lower Fraser Valley during spring. Ph.D. Thesis, University of British Columbia, Vancouver, B.C. Oke, T.R. 1978. Boundary Layer Climates, Methuen and Co. Ltd., New York. Pandolfo, J.P. 1966. Wind and temperature profiles for constant-flux boundry layers in lapse conditions with a variable eddy conductivity to eddy viscosity ratio. J. Atmos. Sci. 23: 495-502. Panofsky, H.A. 1963. Determination of stress from wind and temperature measurements. Quart. J.R. Meteor. Soc. 89: 85-94. Panofsky, H., A.K. Blackadar and G.E. McVehil. 1960. The Diabatic Wind Profile. Quart J. Roy. Meteor. Soc. 86: 390-395. Paulson, CA. 1970. The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer. J. of Appl. Meteorol. 9: 857-861. '105-Philip, J.R. 1957. Evaporation, and moisture and heat fields in the Priestley, C.H.B. 1959. Turbulent Transfer in the Lower Atmosphere, Univ. Chicago Press, Chicago. Richardson, L.F. 1920. The supply of energy from and to atmospheric eddies. Proc. Roy. Soc, London. A97: 354-373. Richtmyer, R.D. 1957. Difference Methods for I n i t i a l Value Problems. Wiley-Interscience, New York. Silen, R.R. 1960. Lethal Surface temperatures and their interpretation for Douglas-fir. Ph.D. Thesis, Oregon State University, Corvallis, Oregon. Smith, F.H. and R.R. Silen. 1963. Anatomy of Heat Damaged Douglas-fir Seedlings. For. Sci. 9: 15-32. Stearns, CR. 1966. Micrometeorological Studies in the Coastal Desert of Southern Peru, Ph.D. Thesis, Univ. of Wisconsin, Madison, Wisconsin. Steel, R.G. and J.H. Torrie. 1960. Principles and Procedures of Statistics With Special Reference to the Biological Sciences. McGraw H i l l , New York. pp. 99-109. Thom, A.S. 1975. Momentum, Mass, and Heat Exchange of Plant Communities. In: "Vegetation and the Atmosphere", Vol. I, (Monteith, J.L., ed.), Aca-demic Press, London, pp. 57-109. Thom, A.S. 1972. Momentum, mass and heat exchange of vegetation. Quart. J.R. Meteor. Soc. 98: 124 Thom, A.S., J.B. Stewart, H.R. Olver and J.H.C. Gash. 1975. Comparison of aerodynamic and energy budget estimates of fluxes over a pine forest. Quart. J.R. Meteor. Soc. 101: 93-10.5. Van Wijk, W.R. and W.J. Derekson. Sinisoidal Temperature Variation in a Layered So i l . In: Physics of Plant Environment, (Van Wijk, ed.), North Wierenga, P.J. and C.T. deWit. 1970. Simulation of heat transfer i n soi l s . Soil Sci. Soc. Amer. Proc. 34: 845-848. -106-A P P E N D I X A D E R I V A T I O N O F T H E ONE D I M E N S I O N A L F I N I T E D I F F E R E N C E S O L U T I O N O F T H E S O I L H E A T FLOW E Q U A T I O N -107-APPENDIX A Derivation of the One Dimensional Finite Difference Solution  Of the So i l Heat Flow Equation The partial differential equation describing the diffusion of heat i n a s o i l profile can be written as follows: -3G = p s c p s 3T_ (1) 8z 9t where T i s the temperature, t is the time, z i s the depth, p s is the s o i l bulk density, C p S is the specific heat of s o i l , and G is the s o i l heat flux density. The equation is derived from the heat balance of a unit volume of s o i l where the change in heat storage ( p s c p s 6T/6t) is equat-ed with the divergence of the s o i l heat flux density (-6G/Sz). The s o i l heat flux density is given by Fourier's Law: G = -k s <5T/6z (2) where k s is the s o i l thermal conductivity and 6T/6z is the s o i l tempera-ture gradient. The negative sign implies a positive heat flux into the s o i l when z is measured positively downward from zero at the surface. Combining equations (1) and (2) gives: C s 9T = 9_ k s 9T 9t 3z 3z (3) where C s (= p s c p s ) is the volumetric heat capacity of the s o i l . The implicit f i n i t e difference solution of (3) (Richtmeyer 1957) involves the discretization of the (z,t) plane in the s o i l profile into a rectangular grid of points j = 1,2,3...L along the z axis and n = 1,2,3... along the t axis. The distance between the vertical nodes is At. The solution at any given grid point is T n j . The partial differentials in equation (3) can be written in difference form as: 3T/3t n-1/ /2 n- 1/ 2 n n-1 = C s (T - T )/At S j j 3 and 3 [ks3T/3z.] Tz j n- 1/ 2 ^ [ ^ j + i / 2 ( " / a . ) j + 1 / 2 l where n- 1/ 2 , n- 1/ 2 [ k s . _ 1 / 2 (3T/3z)_ 1 / 2 ] 3T 3z 3+1/2 1_ Az n-1 n T + T j+1 j+1 n-1 n T + T j 5 and 8T 3z 1/2 n- / 1_ Az f n-1 n T + T j i n-1 n T + T " ( 3-1 3-1 The completed f i n i t e difference form of equation (3) becomes: n n-1 T - T 3 3 At = ! _ < k Az 'j+V2 n-1 n n-1 / T + T - T 2Az V j+1 j+1 j n- 1/ /2 - k. n-1 n n-1 n 1 / T + T - T - T Collecting similar terms and rearranging gives: -T n- 1/ 2 At 2 A z ' - T + T n- 1/ 2 At S j + 1 / 2 n- 1/ 2 At ;s . 1/2 J-V' + Cc n- 1/ /2 2 A z ' n-1 +T n- 1/ 2 At s,a/2 n- 1/ 2 At :s. 1,2 'j+V 2 A z ' 2 A z ' n- 1/ 2 At V / 2 2 A z Z n-1 = T + Ct j+l 2 A z ' n- 1/ 2 At ..1,2 Az' n-1 + T 2 A z 2 (9) At Equation (9) can also be expressed as: n n n - A T + B T - C T = D (10) j j+l j J J j " 1 3 where n- 1/  /2 A = k " A t / 2 A z 2 ( I D 1 ,2 j D j - V 1 /2 n-V Cj = k s j + l / 2 A t / 2 A Z .2 (12) B = A + C + C s (13) j 3 J 3 n-1 / n-D = A T + C s j j j+l V s j .1/2 (U) -770-When equations for the upper and lower boundary conditions ( j = 1 and j = L) are incorporated, the f i n i t e difference approximations in equation (10) constitute a system of L linear equations in L unknowns.The coeffi-cients A,B,C and D vary with the node (j = 1, j = 2 to L -1, j = L), with the time step (n = 1, n = 2, n > 2), and with changes in the heat capacity and thermal conductivity of the profile. The variables on which A,B,C and D depend are the s o i l thermal properties, the node size in the grid, and the temperatures at the upper and lower boundary conditions. The values of T n j are calculated from the recurrence relation: n n T = E T + F j j J+1 J where and E = A / j J B - C E j j J-1 D + C F j j J-1 / B - C E j J J-1 (15) (16) (17) The E and F coefficients are calculated from nodes j = 1 to j = L, and the T's are back-calculated from j = L to j = 1. The coefficients A-F at node j = 1 and A-F at node j = L are defined at each new time step as follows: A (1) = 0 B (1) = 1 C (1) = 0 D (1) = T (t, 1) E (1) = A (1) / B (1) F (1) = D (1) / B (1) ( l g ) A (L) = 0 B (L) = 1 C (L) = 0 D (L) = T ( t , L) E (L) = A (L) / B (L) F (L) = D (L) / B (L) -7J2-APPENDIX B MODELLING CLEAR SKY SOLAR IRRADIANCE ON VARIABLE SLOPES AND ASPECTS -773-APPENDIX B. Modelling Clear Sky Solar Irradiance on Variable Slopes and Aspects 1. INTRODUCTION Potential solar irradiance is the maximum short wave irradiance poss ble for any given observation site and time. It only occurs under clear skies and varies with latitude, altitude, slope, aspect, and time of the day and year. It can be calculated from the geometrical relationship between the sun and earth and an empirical atmospheric radiation attenua-tion function as follows: ^t slope = b^ slope + $d slope where St si 0pe> b^ slope> a n c^ $d slope a r e t n e t o t a l solar i r r a d i -ance, direct and diffuse components respectively. For a horizontal surface, the direct component is calculated as follows (List 1971): sb horiz = sp c o s z (2) where z is the zenith angle of the sun. Sp, the direct solar beam irradiance can be empirically determined as follows (List 1967): S p = a* S p o (3) where Sp D is the extraterrestrial direct solar beam irradiance (solar constant = 1360 W m ), a is an atmospheric transmission coefficient, and is the optical airmass number. The atmospheric transmission coefficient varies from 0.9 for a very clear atmosphere to 0.6 for a very hazy or smoggy atmosphere. A typical value of (a) for clear days would be 0.84 (Campbell 1977). The optical airmass number, the ratio of the slant path length of solar radiation to the zenith path length in the atmosphere, for elevation angles greater than 10 degrees is given as follows: m = (P/PQ) / (cos z (4) where (P/PQ) , the atmospheric pressure at the observation site divided by sea level atmospheric pressure, corrects m for altitude effects. The zenith angle of the sun (z) can be found from: cos z = sin c|> sin e + cos § cos e cos 15 (t - t Q) (5) where tj> is the latitude, e is the solar declination, t is the time of day, and t 0 is the time of solar noon. The horizontal direct short wave irradiance can be corrected for slope and aspect as follows (Hay 1979): Sb slope = Sb froriz cos l cos z where i , the angle of incidence of the sun's rays on the slope, is found from: cos i = cos s cos z + sin s sin z cos (az-b) where s is the inclination angle of the slope, b is the azimuth of the slope (aspect) and az, the azimuth of the sun, is found from: cos az = sin § cos z - sin e ^ cos <j> sin z To a good approximation the clear sky diffuse radiation can be found from half the difference between the irradiance on a horizontal surface -7/5-below and above the atmosphere as follows (List 1971): S d = 0.5 S p o (1 - a m) cos z (8) Sjj i s almost a constant fraction of S p o over most of the range of zenith angles (Campbell 1977), hence i t does not have to be corrected for slope and aspect, i.e. S^  slope " ^d* Reflected radiation might also increase the total solar irradiance by a few percent depending upon the sky view factor of the observation site. -us-1200 900 E ui < 600 Q < tr cr < o to 300 SOUTH 50* HORIZONTAL \ EAST 30* * WEST SO* \ NORTH SO* \ / J s L 100 200 JULIAN DAY 300 Calculated solar irradiance at solar noon on h o r i z o n t a l , south 30°, Ea and West 30°, and North 30° slopes. -117' L i s t of Symbols _ o S t siope = total solar irradiance on a sloping surface W m _ 2 Sb slope = direct solar irradiance on a sloping surface W m _ 2 sb horiz = direct solar irradiance on a horizontal surface W m _ 2 S,j = diffuse solar irradiance W m _ 2 ^d slope = diffuse solar irradiance on a sloping surface W m _2 Sp = direct solar beam irradiance W m Spo = solar constant (1360 ± 40 W m~2) W m - 2 i = angle of incidence of sun1s rays on a slope degrees z = zenith angle of sun degrees s = inclination of slope from horizontal degrees b = azimuth of slope (aspect) degrees az = azimuth of sun degrees p = latitude of observation site degrees e = declination of sun degrees P = atmospheric pressure at observation site Pa P Q = atmospheric pressure at sea level Pa m = optical airmass number dimensionless a = atmospheric transmission coefficient dimensionless t = time of day h t n = solar noon h APPENDIX C SOIL THERMAL PROPERTIES APPENDIX C Soil Thermal Properties The s o i l thermal properties used in the simulations presented in this thesis were chosen on the basis of s o i l texture and s o i l water content. The volumetric heat capacity was calculated from the following equa-tion (De Vries 1963): c s = V^m + xo co + V^w + x a c a ^ where the subscripts m, o, w and a represent mineral s o i l , organic matter, water, and air, respectively, x is the volume fraction and C is the volu-metric heat capacity. Volumetric heat capacities of mineral soils having differing bulk densities and water contents are shown in Figure 1. The thermal diffusivity of a mineral s o i l sample taken from the Cameron Valley clearcut site was measured using the method described by Sterling and Taylor (1965). Measurements taken at volumetric water con-tents of 0.0 and 0.2 agreed very well with the d i f f u s i v i t i e s of sand quoted by van Wijk and de Vries (1963) (Figure 2). The thermal conductivity was calculated from the volumetric heat capacity and the thermal diffusivity as follows: k s = K SC S (2) where k s is the thermal conductivity, Ks is the thermal diffusivity, and C s is the volumetric heat capacity. The calculated values of thermal conductivity are plotted against s o i l water content in Figure 3. -\io-T 1 -i r S o i l volumetric heat capacity versus s o i l volumetric water content for a mineral s o i l . 1.0 » 0.8 E (0 I o 0.6 >-CO Z> u. t 0.4 o cc LU X 0.2 0.1 02 1 I 1 1 1 1 1 — A A A A A i O — — O o O 0 o O O — o O O O < — — o — o 1 1 1 1 1 1 . 1 SOIL VOLUMETRIC WATER CONTENT (mVm3) Soil thermal dif f u s i v i t y versus volumetric water content of a sandy s o i l having a bulk density of 1.0 ( Q ) , 1.3 ( O )» and 1.6 ( A ) Mg m~3. 0 0.1 0.2 03 0.4 SOIL VOLUMETRIC WATER CONTENT (mVm5) Soil thermal conductivity versus volumetric water content of a sandy having a bulk density of 1.0 (Q>» 1 , 3 ( O )» a n d 1.6 ( A ) Mg m -3 c> APPENDIX D COMPUTER PROGRAM USED TO CALCULATE SOIL TEMPERATURES AND SURFACE ENERGY FLUXES 200 230 240 250 270 280 290 310 320 330 340 350 360 380 390 400 410 420 500 REM 510 REM 515 518 520 522 524 526 530 535 538 540 542 544 546 548 571 575 576 580 590 600 610 615 620 630 780 790 800 810 820 830 840 REM SET UP MODEL INITIAL CONDITIONS AND CONSTANTS****************** OPEN"I",#1,"SOILTEST.DAT" • TO READ IN CLIMATE DATA OPEN"0",#2,"LONG.DAT" • TO OUTPUT SOILTEMP DATA DIM A(161),B(161),C(161),D(161),E(161),F(161),COND(161),RHOCP(161 TEMP(2,161),DEPTH(161) INPUT"MEAN DAILY SURFACE SOIL TEMPERATURE "jTBAR INPUT"SURFACE ROUGHNESS LENGTH 20 (M) 0.001-lMM "jZO INPUT"LE • X * RS SIGMA ESURF ALBDO GRAV K Zl DT NODES DZ T PI SET UP 5.67E-08 .95 .16 9.8 .41 1 60 161 .005 2 3.14159 •;HC INPUT X "; ALF STEPHAN BOLTZMANN CONSTANT SURFACE EMISSIVITY SURFACE ALBEDO GRAV ACCEL VON KARMAN CONSTANT MEASUREMENT HEIGHT MINUTES BETWEEN DATA SAMPLES NUMBER OF NODES IN THE PROFILE 80 CM PROFILE WITH .5 CM NODES TIME VARIABLE • 4159 ' PI  f INITIAL AND BOUNDARY CONDITIONS OF THE SOIL PROFILE THERMAL PROPERTIES OF THE PROFILE PRINT"FROM THE SURFACE DOWN TO THE 2 CM DEPTH " INPUT"ENTER THE VOLUMETRIC HEAT CAPACITY 1.6 (MJ/M3 K) ";HC INPUT"ENTER THE THERMAL CONDUCTIVITY 1.32 (W/M K) "jTC FOR Z-NODES TO (NODES-4) STEP -1 COND(Z) - TC 'J/S M C RHOCP(Z)- HC *lE+06 'J/M3 C NEXT Z PRINT"FROM THE 2 CM SOIL DEPTH DOWN TO THE 5 CM DEPTH " INPUT"ENTER THE VOLUMETRIC HEAT CAPACITY 1.6 (MJ/M3 K) INPUT"ENTER THE THERMAL CONDUCTIVITY 1.32 (W/M K) ";TC FOR Z-(NODES-5) TO (NODES-10) STEP -1 COND(Z) - TC 'J/S M C RHOCP(Z)- HC *lE+06 «J/M3 C NEXT Z PRINT"FROM THE 5 CM SOIL DEPTH DOWN TO THE 25 CM DEPTH INPUT"ENTER THE VOLUMETRIC HEAT CAPACITY 1.6 (MJ/M3 K) INPUT"ENTER THE THERMAL CONDUCTIVITY 1.32 (W/M K) "jTC FOR Z-(NODES-ll) TO (NODES-50) STEP -1 COND(Z) - TC *J/S M C RHOCP(Z)" HC *lE+06 'J/M3 C NEXT Z PRINT"FROM THE 25 CM SOIL DEPTH DOWN TO 80 CM " INPUT"ENTER THE VOLUMETRIC HEAT CAPACITY 1.6 (MJ/M3 K) INPUT"ENTER THE THERMAL CONDUCTIVITY 1.32 (W/M K) "»TC FOR Z- (NODES-51) TO 0 STEP -1 COND(Z) - TC 'J/S M C RHOCP(Z)- HC*lE+06 'J/M3 C NEXT Z DD • SQR(COND(15)/RHOCP(15)*864001/PI) PRINT'DAMPING DEPTH (CM) « "jDD Z • -DZ jHC ;HC -IZ5-850 POR I • NODES TO 1 STEP -1 860 Z - Z • DZ 870 DEPTH(I) - Z 880 TEMP(1,I)-TBAR • +10*EXP(-DEPTH(I)/DD)*SIN(2*PI - DEPTH(I)/DD) 885 TEMP(2,1)-TBAR 890 REM IF I<200 THEN PRINT USING"##.#### "jI,DEPTH(I)jTEMP(1,1),COND(I),RHC 900 IF I > NODES THEN LPRINT USING"If.f ";TEMP(l,I)j 910 NEXT I 911 LPRINT 912 PRINT 914 LPRINT" " j 915 LPRINT USING"\ \";"TIME">"TA"j"U"j"RI"> 917 LPRINT OSING"\ \";"RS";"RN";"G";"H";"LE"j 920 LPRINT USING"\ \";"T0";"T5";"T10">"T15";"T20"j"T25";"T30";"T35";"T40" 925 LPRINT 1000 REM*********************MAIN PROGRAM STARTS HERE**********"""** 1010 WHILE NOT EOF(l) 1020 INPUT!1,DAY,TIME,T0,T1,T10,T15,T30,T50,TA,RS,U 1030 PRINT" T« "jTIME" SOLAR - "RS» 1050 PRINT"TAIR - "TA" WIND • "U 1060 GOSUB 2500 'CALCULATE VARIABLES THAT DEPEND ON RS,TA,U 1120 WHILE ABS(GO-GA) > 1 1140 IF GA>G0 THEN GA=GA-(ABS(GA-G0))/2 'ESTIMATE A NEW 1145 IF G0>GA THEN GA«=GA+ (ABS (GA-G0)) /2 'VALUE OF GA 1160 GOSUB 3000 'ITERATE FOR H AND TS 1180 GOSUB 4000 'ITERATE FOR GO 1182 PRINT "RN,LE,GA,H1,TS"; 1183 PRINT USING"!!!!.!! ";RN,LE,GA,Hl,TS 1185 PRINT "RN,LE,G0,H1,TS"j 1187 PRINT USING"####.#t RN,LE,GO,HI,TEMP (T-l,NODES) 1192 PRINT "RI ,PHIH,PHIM,USTAR"j 1193 PRINT USING"!!I.lttt ";RI,PHIH,PHIM,USTAR 1194 REM PRINT TWOGASAGO,LASTGA,GA 1195 IF TWOGASAGO - GA THEN GA - (GA+LASTGA)/2 1196 TWOGASAGO - LASTGA 1197 LASTGA - GA 1199 REM PRINT TWOGASAGO,LASTGA,GA 1200 PRINT 1205 WEND 1210 PRINT 1220 REM CHANGE FROM ONE HOUR TO NEXT WHEN RN-H+LE+G**** 1230 FOR Z - 1 TO NODES 1240 TEMP (1,Z) - TEMP(2,2) 1250 NEXT Z 1 3 0 0 REM PRINT DATA TO FILE AND PRINTER 1320 LPRINT USING"!!.t ">TIME,TA,U,RI} 1330 LPRINT USING"!!!!.! "jRS,RN,H1,GA,LE; 1360 REM FOR Z - NODES TO 143 STEP -1 -IZ6-1370 1380 1385 1387 1390 1400 1405 1407 1410 1500 2400 2440 2500 2520 2540 2560 2580 2620 2640 2650 2660 2680 2740 2750 2760 2880 3000 3010 3020 3030 3040 3050 3060 3070 3080 3100 3110 3120 3140 3160 3180 3200 3210 3400 REM LPRINT USING"it.it ";TEMP(T,Z); REM NEXT Z REM LPRINT:LPRINT" PRINTi 2,USING"ft"»TIME; FOR Z « NODES TO 61 STEP -10 LPRINT USING" ##.I";TEMP(T,Z); PRINT42,USING" it"jDEPTH(Z)*100j PRINTi2,USING" f(.#"jTEMP(T,Z); NEXT Z LPRINT:PRINTi2,"" WEND CLOSE : END REM SUBROUTINE TO SET VARIABLES****************************** - 1.00055-(TA*(-.0001516)) 'DENSITY OF AIR « 1300.05 - (TA*4.07) 'SPECIFIC HEAT OF AIR • RHOA*CPA 'HEAT CAPACITY OF AIR - 1-(.261*EXP((TA*2)*-.000777)) 'SKY EMISSIVITY • EAIR*SIGMA*(273.2+TA)"4 'DOWNWARD LONGWAVE FLUX • 0.15*RS-40 'ESTIMATED GO - GA-10 - .015*(TA+3)*RS » ALF*RS 'NOTE ALF RANGES FROM 0 TO 1 - 100 'RESET HI » 20 'RESET H2 • TA -5 "RESET TS RETURN 'TO MAIN PROGRAM REM SENSIBLE HEAT TRANSFER SUBROUTINE a************************** WHILE H2 < HI GOSUB 3500 'ESTIMATE H2 IF HI > H2 THEN TS - TS + 1 WEND RHOA CPA RHOCPA EAIR RLD GA GO REM LE LE HI H2 TS RETURN WHILE H2 > HI GOSUB 3500 'ESTIMATE H2 IF HI < H2 THEN TS = TS - .25 WEND WHILE H2 < HI GOSUB 3500 'ESTIMATE H2 IF HI > H2 THEN TS - TS + .05 WEND WHILE H2 > HI GOSUB 3500 'ESTIMATE H2 IF HKH2 THEN TS • TS - .01 WEND RETURN ' TO MAIN PROGRAM -111-3500 REM SUBROUTINE TO ESTIMATE H2************************************ 3520 TK - (TS+TA)/2 +273.15 3540 RLD - ESURF"EAIR*SIGMA*(273.15+TA)"4 3560 RLU • ESURF*SIGMA*(273.15+TS)~4 3580 HI •> RS*(1-ALBDO)+RLD-RLU-LE-GA 3600 RN - RS*(1-ALBDO)+RLD-RLU 3610 RI - -(GRAV/TK)*(TS-TA)*Z1/U~2 • REM TO MULT LOG(Z/Z0) TO Z 3612 IF (16*RI)>11 THEN RI - .0625 3615 X - <l-16*RIp.25 3620 PHIH - 2* LOG((l+X~2)/2) 3622 IF RI > 0 THEN PHIH - 4.7* ( RI/(1-4.7*RI)) 'PHIH • 4.7 ZETA 3625 PHIM - 2*LOG((l+X)/2) + LOG((1+X~2)/2) - 2* ATN(X) + PI/2 3627 IF RI > 0! THEN PHIM - PHIH 3630 USTAR - (K*U) / ( LOG(Z1/Z0) - PHIM) 3635 HTRANS - K*USTAR/(LOG(Z1/Z0) - PHIH) 3700 H2 - RHOCPA*HTRANS*(TS-TA) 3730 RETURN * TO SUBROUTINE ITERATION 4000 REM SUBROUTINE SOILTEMP CALCULATES GO AND PROFILES OF TEMPERATURE 4020 TEMP(T.NODES)- TS 4060 DZZ - 2* (DZ~2) 4080 DTT - 3600 'SECONDS PER HOUR 4100 A(l) * 0 4120 B(l) - 1 4140 C(l) - 0 4160 D(l) « TEMP(T,1) 4180 E(l) - A(l)/B(l) 4200 F(l) - D(l)/B(l) 4220 A(NODES)- 0 4240 B (NODES)- 1 4260 C(NODES)- 0 4280 D(NODES)- TEMP(T.NODES) 4300 E(NODES)- A(NODES)/B(NODES) 4320 F(NODES)- D(NODES)/B(NODES) 4340 FOR Z - 2 TO NODES-1 4360 C(Z) - COND(Z-l) *DTT / DZZ 4380 A(Z) • COND(Z) *DTT / DZZ 4400 B(Z) - RHOCP(Z) + A(Z) + C(Z) 4420 M - A(Z) * TEMP(T-1,Z+1) 4440 N » (RHOCP(Z) - A(Z) - C(Z)) * TEMP(T-1,Z) 4460 O - C(Z) * TEMP(T-1,Z-1) 4480 D(Z) - M + N + O 4500 NEXT Z 4520 FOR Z - 2 TO NODES-1 4540 BCE - B(Z) - ( C(Z)*E(Z-1) ) 4560 E(Z)- A(Z) / BCE 4580 P(Z)- (D(Z) + (C(Z) » F(Z-l))) / BCE 4600 NEXT Z 4620 FOR Z - NODES-1 TO 2 STEP -1 4640 TEMP(T.Z) •( E(Z) • TEMP(T,Z+1) ) + F(Z) 4660 NEXT Z 4740 GO - COND(NODES)* (TEMP(T.NODES) - TEMP(T,NODES-1)) / DZ 4750 GI - COND(NODES-1)*(TEMP(T,NODES-1) - TEMP (T.NODES-2)) / DZ 4760 REM PRINT"G0 - "jGO" GAIR » "»GA" TS - "»TS" T(T,TOP) • "% TEMP(T,1) 4800 RETURN 'TO MAIN PROGRAM 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0095699/manifest

Comment

Related Items