A solar radiation algorithm for ecosystem dynamic models

Ecological Modelling, 61 (1992) 149-168

149

Elsevier Science Publishers B.V., Amsterdam

A solar radiation algorithm for ecosystem dynamic models Nedialko T. Nikolov and Karl F. Zeller US Forest SerL,ice, Rocky Mountain Forest and Range Experiment Station, 240 West Prospect, Fort Collins, CO 80526, USA (Accepted 12 November 1991)

ABSTRACT Nikolov, N.T. and Zeller, K.F., 1992. A solar radiation algorithm for ecosystem dynamic models. Ecol. Modelling, 61: 149-168. A general algorithm is described for estimating average monthly solar radiation in cal cm-2 day- l received on mountain slopes that uses basic topographic and climatic information for input, i.e. latitude, elevation, slope aspect and orientation, and average monthly data for ambient temperature, relative humidity and total precipitation. The algorithm is an extension of the methodology developed by Lui and Jordan, Klein, and Bonan. We have tested our method against independent data from 69 meteorological stations throughout the northern hemisphere provided by Mtiller and Zeller. The stations were selected to represent different latitudes, climate zones and elevations. The test showed that the algorithm predicts quite accurately seasonal patterns of solar radiation from the subpolar region down to the tropics and thus can be used with ecological simulation models (e.g. gap-phase succession models or regional landscape models).

INTRODUCTION

Solar radiation is the energy source for important ecological processes such as evapotranspiration and photosynthesis which control to a great extent the distribution, type and physiognomy of terrestrial vegetation. Many recent ecological models include solar radiation as a critical environmental driver (e.g. Running, 1984; Woodward, 1987; Bonan, 1988; Running Correspondence to: N.T. Nikolov, US Forest Service, Rocky Mountain Forest and Range Experiment Station, 240 West Prospect, Fort Collins, CO 80526, USA. 0304-3800/92/$05.00 © 1992 - Elsevier Science Publishers B.V. All rights reserved

150

N.T. NIKOLOV AND K.F. ZELLER

and Coughlan, 1988; Leemans and Prentice, 1989; Bonan, 1991). Most of these models implement algorithms to estimate incident radiation since measured data are rarely available for the site of interest. The radiation algorithms we know (e.g. Garnier and Ohmura, 1968; Swift, 1976; Klein, 1977; Bristow and Campbell, 1984; Running et al., 1987; Bonan, 1988) require site-specific calibrations in order to account for the radiation attenuation effects of the earth's atmosphere at a particular location. The calibrations usually employ actual data from a nearby meteorological station for the incoming solar radiation on a horizontal surface a n d / o r for cloudiness. Since atmospheric effects on radiation vary greatly with latitude, elevation and climatic conditions, and because most meteorological stations do not typically (or continuously) record solar radiation a n d / o r cloudiness (Miiller, 1982), these algorithms have limited applicability or generality for the ecological models they are used with. We describe in this paper a general method for estimating monthly average solar radiation on mountain slopes that uses worldwide commonly recorded topographic and climatic data for input, i.e. latitude, elevation, slope, aspect, mean monthly ambient temperature, mean monthly ambient relative humidity, and total monthly precipitation. We have tested our method against independent data from 69 meteorological stations throughout the northern hemisphere. The close agreement between estimated and measured data called forth our decision to publish this algorithm to be used with ecological simulation models. The Turbo Pascal code of the algorithm is provided in the Appendix. ALGORITHM Our algorithm is based on the methodology developed by Lui and Jordan (1960, 1962, 1963) and Klein (1977) and applied by Bonan (1988). In their approach the calculated potential solar radiation on a horizontal surface outside the atmosphere is attenuated by atmospheric effects to produce the total radiation received on a horizontal surface at the earth's surface. This radiation is then decomposed into its direct and diffuse components which are subsequently adjusted using various tilt factors to components of the surface of interest. We have extended the above methodology by functionally relating the atmospheric transmittance to latitude, elevation and climatic conditions of the site. We have also developed an approach for estimating monthly average cloud cover from mean monthly temperature, mean monthly relative humidity and total monthly precipitation. By substituting our relationship for the cloudiness data required by the original methodology for input, we have actually expanded the algorithm's applicability as temperature,

SOLAR R A D I A T I O N A L G O R I T H M FOR ECOSYSTEM D Y N A M I C MODELS

15 |

relative humidity and precipitation are the most commonly measured climatic parameters throughout the world.

Mean monthly solar radiation on the top of the atmosphere The algorithm estimates average monthly solar radiation received on a horizontal surface outside the atmosphere by dividing the monthly sum of daily potential radiation by the number of days in each month. Daily potential radiation is a function of latitude, solar declination, s u n r i s e / sunset hour angle, and day of year (Klein, 1977): R 0 = S c 458.3711 + 0.033 cos(360J/365)] [cos(L) cos(Ds) sin(hs) + (hs/57.296) sin(L) sin(Ds) ] where R 0 is the solar radiation received on a horizontal plane on the top of the earth's atmosphere (cal cm : day-l); Sc is the solar constant (2.0 cal cm -2 min 1); j the Julian day of year; L the latitude (degree); Ds the solar declination (degree); hs the s u n r i s e / s u n s e t hour angle (degree). We estimate solar declination in radians (Swift, 1976): D S= arcsin{0.39785 sin[4.868961 + 0.017203J +0.033446 sin(6.224111 + 0.017202J)]} Solar s u n r i s e / s u n s e t hour angle is computed by (Keith and Kreider, 1978): h~ = a r c c o s [ - t a n ( L ) tan(Ds)]

Mean monthly solar radiation on the earth's surface The amount of solar radiation received on the earth's surface is less than that at the top of the atmosphere because of the scattering and absorbing effects of air molecules, water vapor, dust, ozone, and carbon dioxide. These attenuation effects can be accounted for on a monthly basis at a particular site as a linear function of undepleted solar radiation and the average portion of the sky covered by clouds (Bonan, 1988): R = Rom(/3 - ~ C )

-a

(1)

where R is the mean monthly solar radiation received on a horizontal surface at the earth's surface (cal cm 2 day-l); Ro m the mean monthly solar radiation incoming to a horizontal plane on the top of the atmosphere (cal cm -z day-1); C the mean monthly cloudiness (tenths); and a,/3 and oare empirical parameters.

152

N.T. N I K O L O V A N D K.F. Z E L L E R

Our analysis showed that mean monthly cloudiness is related to the ratio of mean monthly surface water vapor pressure in the air and total monthly precipitation. Non-linear regression of climate data from Bulgaria (Bulgarian Climate Information System, BULCIS, Nikolov and Bojinov, personal communication, 1991) where cloudiness has been continuously recorded and complete data records are available yielded the following relationship: C = 1 0 - 1.155(ev/P) °5

(2)

where e v is mean monthly water vapor pressure in the atmosphere (Pa) and P is total monthly precipitation (mm). Mean monthly absolute humidity (e v) is estimated from monthly averages of observed surface temperature and relative humidity by applying Dilley's equation (Dilley, 1968) for saturation air water vapor pressure: e v = H r 6.1078 exp[17.269T/(T + 237.3)] where H r is mean monthly relative humidity (%) and T is mean monthly temperature (°C). The empirical parameters of equation (1) can be interpreted in terms of atmospheric properties. /3 is clear-sky atmospheric transmittance, o- accounts for the light absorption effects of clouds, and a is related to the diffuse solar radiation. As solar elevation decreases with latitude so does the atmospheric transparency since lowering solar altitude angle causes solar rays to pass through thicker air masses that attenuate more solar energy. A thicker air layer, however, increases the proportion of diffuse radiation (Rayleigh scattering) since a greater part of the solar beam is scattered by air molecules. Therefore, one could expect that both /3 and a will decrease with latitude. This conclusion is supported by analysis of climatic and radiation data from Mfiller (1982) for the northern hemisphere. We used data from two meteorological stations (i.e. Lansing, Michigan, USA, 42°44'N, 84°29'W, 261 m above sea level; and Jerbogacon, USSR, 61°16'N, 108°01'E, 287 m above sea level) of similar altitude and different latitudes to derive expressions for a, /3 and o- with respect to latitude: a = 32.9835 - 64.88411 - 1.3614 cos(L)]

(3)

/3 =0.715 -0.318311 - 1.3614 cos(L)[

(4)

or = 0.03259

(5)

The analysis indicates that ~r tends to remain constant.

SOLAR RADIATION ALGORITHM FOR ECOSYSTEM DYNAMlC MODELS

153

Solar radiation as a function of site elecation

Equations (1) through (5) account for the depletion of mean monthly solar radiation received on a horizontal surface outside the atmosphere to the radiation received at 274 m above sea level [i.e. (261 + 287)/2] on the earth's surface. Since irradiance normally increases with elevation, a correction for the actual altitude of the site of interest is required. We have derived an elevation correction formula that is based on two assumptions: (i) light attenuation through the atmosphere is exponential with depth and can be described by Beer's law; and (ii) radiation extinction coefficient is a function of solar elevation: R a = R + [(Rom + 1 ) - R ] ( 1 - e x p [ - k / s i n ( E ~ m ) ( Z - 2 7 4 ) / 2 7 4 ] }

where R a is the mean monthly solar radiation received on a horizontal surface at the altitude of the site of interest (cal cm -2 day 1); R the mean monthly solar radiation received on a horizontal surface 274 m above sea level (cal c m - 2 day-l); k the basic atmospheric radiation extinction coefficient (0.014 cal c m - 2 m 1); E s m the average monthly solar altitude angle (degree); and Z the site altitude (m above sea level). The basic light extinction coefficient, k, was estimated by non-linear regression of calculated and measured radiation data and calculated solar elevation data for Brooklyn Lakes, a site within the US Forest Service Glacier Lakes Ecosystem Experiment Site (GLEES) in Northern Rocky Mountains, Wyoming, USA. The site is located at 41022 ' North latitude and has an altitude of ca. 3170 m (Musselman, in press). Mean monthly solar altitude angle (Esm) is calculated by dividing the monthly integral of hourly estimates of solar elevation by the number of hours in a month when the sun is above the horizon. Solar elevation (in degrees) for a particular hour and day (E~) is computed by (Keith and Kreider, 1978): Es = arcsin[sin(L) sin(Ds) + cos( L ) cos(Ds) cos(h)]

(6)

and h = ( t - 12)15 °

(7)

where h is the solar hour angle (degrees) and t is the time (hours) from midnight. Mean monthly solar radiation on a tilted surface

This part of the algorithm follows the approach of Keith and Kreider (1978) as applied by Bonan (1988). The estimated mean monthly solar

154

N.T. NIKOLOV

A N D K.F. Z E L L E R

radiation coming in on a horizontal surface (R a) is decomposed into its direct and diffuse components w h i c h are then adjusted by using corresponding tilt factors to components of the surface of interest. A tilt factor is the ratio of solar radiation received on a tilted surface to that received on a horizontal plane. The fraction of the incoming radiation that is diffuse radiation is calculated based on the fraction ( F r) of solar radiation transmited through the atmosphere: F r = Ra//Rom and Ro/R

a

= 1.0045 + 0.0435F r - 3.5227Fr 2 + 2.6313Fr 3

w h e r e R d is diffuse radiation (cal cm -2 day -1) received on a horizontal surface at the surface of the site of interest. If F r > 0.75, then R o / R a = 0.166. For the direct radiation, the instantaneous tilt factor T b is: T b = coS(/s.)/sin(Es) where I s is the solar incidence angle and E s the solar altitude angle [equation (6)]. The long-term direct radiation tilt factor is calculated by dividing the daily average of cos(I s) by the daily average of sin(E s) for the period in which the sun is both above the horizon ( E s > 0) and in front of the surface ( I s < 90) for all days of the month. In our method, monthly averages of cos(I s) and sin(E s) are obtained from hourly estimates of solar altitude angle [equations (6) and (7)] and solar incidence angle. The solar incidence angle is a function of latitude, solar declination, solar hour angle and surface orientation. For horizontal surfaces: coS(Is) = sin(Es) For surfaces facing due south with a slope inclination angle S: I s = arccos[sin( L - S) sin(Ds) + cos(L - S) cos(Ds) cos(h)] For tilted surfaces facing in other directions: I s = arccos[cos(A s - A w ) cos(Es) sin(S) + sin(Es) cos(S)] where A s is the solar azimuth angle (degrees) and A w the slope aspect (degrees). The azimuth angle of the sun from the south is calculated by: A s = arcsin[cos(Ds) s i n ( h ) / c o s ( E s ) ] Both A s and A w are positive east of south and negative west of south, w h e r e due south is zero.

SOLAR RADIATION ALGORITHM FOR ECOSYSTEM DYNAMICMODELS

155

TABLE 1 Solar radiation m o d e l test for N o r t h America. Climate zones are after C. Troll and K.H. Paffen (MiJller, 1982) Station name

Latitude north

Altitude (m)

C o l d - t e m p e r a t e boreal G. Lakes 41022 ' 3292 Cool-temperate Seattle 47°36 ' 4 Portland 45032 ' 47 Spokane 47047 ' 718 S.S. Marie 46°28 ' 220 Boston 42013 ' 192 Bismark 46046 ' 511 D o d g e City 37036 ' 791 Oklahoma 35029 ' 382 R a p i d city 44002 ' 993 Boise 43°34 ' 866 Lander 42°49 ' 1696 S.L. City 40046 ' 1286 Elko 40050 ' 1547 Chicago 41047 ' 185 Omaha 41018 ' 337 A l b u q u e r q u e 35003 ' 1620 Cleveland 41024 ' 237 New York 40047 ' 96 Indianopolis 39044 ' 241 Washington 38°51 ' 4 Nashville 36°07 ' 176 Fresno 36046 ' 100 W a r m - t e m p e r a t e subtropical Greensboro 36005 ' 273 Little Rock 34044 ' 78 Atlanta 33039 ' 308 Charleston 32054 ' 12 Dallas 32°51' 146 New O r l e a n s 29°57 ' 3 Apalachikola 29044 ' 4 San A n t o n i o 29°32 ' 241 Tampa 27057 ' 11 Los A n g e l e s 34003 ' 103 Brownsville 25054 ' 5 Phoenix 33°26 ' 340 E1 Paso 31°48 ' 1 194 Miami 25048 ' 2

Climate zone

r a

s b

C o n t i n e n t a l boreal

0.987

1.00

20.1

Oceanic Oceanic Sub-continental Sub-continental P e r m e n e n t l y humid H u m i d steppe Humid steppe Humid steppe Dry s t e p p e Dry steppe Dry s t e p p e Dry s t e p p e Dry s t e p p e Humid-and-warm summer Humid-and-warm summer Humid-and-warm summer P e r m a n e n t l y humid Permanently humid P e r m a n e n t l y humid P e r m a n e n t l y humid Permanently humid Dry-summer steppe

0.994 0.984 0.996 0.987 0.996 0.996 0.996 0.996 0.988 0.995 0.943 0.986 0.986 0.998 0.996 0.992 0.963 0.995 0.994 0.991 0.969 0.991

1.23 1.19 0.94 0.89 1.12 0.91 0.90 1.02 0.97 0.98 0.93 0.98 1.11 1.04 0.97 0.98 0.86 1.18 0.93 0.99 1.03 1.00

- 18.4 -30.0 13.8 4.6 -60.1 2.2 12.2 -40.5 6.8 13.1 23.3 20.1 - 75.1 48.3 -20.1 -9.4 53.6 - 22.2 15.6 -4.1 14.6 7.6

Permanently humid P e r m a n e n t l y humid P e r m a n e n t l y humid Permanently humid P e r m a n e n t l y humid Permanently humid P e r m a n e n t l y humid P e r m a n e n t l y humid P e r m a n e n t l y humid Dry-summer steppe Steppe S e m i - d e s e r t and desert S e m i - d e s e r t and desert Tropical h u m i d - s u m m e r

0.997 0.992 0.990 0.992 0.990 0.976 0.992 0.985 0.983 0.974 0.919 0.985 0.957 0.996

1.09 1.09 1.08 1.06 1.07 1.29 1.05 0.99 1.01 1.37 0.89 1.09 1.13 1.01

-42.2 -37.7 -33.9 -29.2 - 53.6 -39.7 -44.3 25.8 - 14.4 - 143.7 123.4 -52.3 -49.0 - 20.5

a P e a r s o n ' s correlation coefficient b e t w e e n e s t i m a t e d and m e a s u r e d data. b Slope of the regression line of e s t i m a t e d versus m e a s u r e d data. c Line intercept.

i "

156

N.T. N I K O L O V A N D K.F. Z E L L E R

TABLE 2 Solar radiation m o d e l test for Eurasia and Africa. Climate zones are after C. Troll and K.H. Paffen (Miiller, 1982) Station name

Latitude north

Altitude (m)

Polar and subpolar Velen 66°10 ' 7 C o l d - t e m p e r a t e boreal Louchi 66005 ' 94 A r c h a n g e l s k 64030 ' 4 Sverdlovsk 56°44 ' 282 Novosibirsk 55002 ' 162 Irkutsk 52016 ' 468 Turuchansk 65047 ' 45 Ojmakon 63016 ' 740 Kirensk 57046 ' 256 Cool-temperate Trondheim 63°25 ' 133 Bergen 60012 ' 45 Paris 48°58 ' 52 Dresden 51007 ' 246 Lyon 45043 ' 200 Bucarest 44°25 ' 82 Scopje 42000 ' 245 Stockholm 59021 ' 44 Moscow 55045 ' 156 Kaunas 54°53 ' 75 Kishinev 47001 ' 95 Omsk 54°26 ' 105 Tbilisi 41014 ' 490 Vladivostok 43o07 ' 138 Fergana 40°23 ' 578 Dushanbe 38035 ' 824 W a r m - t e m p e r a t e subtropical Lisabon 38°~'V 77 Palermo 38°07 71 Athinai 37o58 ' 107 Sevilla 37024 ' 30 Tropical Madras 13°04 ' 16 Bur Sudan 19o34 ' 5

Climate zone

r

Subarctic t u n d r a

0.974

0.95

- 1.5

C o n t i n e n t a l boreal C o n t i n e n t a l boreal C o n t i n e n t a l boreal Continental boreal C o n t i n e n t a l boreal Highly cont. boreal Highly cont. boreal Highly cont. boreal

0.994 0.990 0.987 0.993 0.992 0.976 0.980 0.990

0.94 0.88 0.93 1.04 0.94 0.76 0.92 1.08

-6.9 5.6 14.0 - 18.4 28.1 3.8 1.8 -7.3

Oceanic Oceanic Oceanic Sub-oceanic Sub-oceanic Sub-oceanic Sub-oceanic Sub-continental Sub-continental Sub-continental Sub-continental Continental Continental Permanently humid Humid-summer steppe H u m i d and warm s u m m e r

0.997 0.992 0.999 0.998 0.999 0.993 0.994 0.989 0.984 .0.987 0.990 0.987 0.988 0.919 0.987 0.987

0.93 0.93 1.07 1.06 0.91 0.90 1.04 1.06 0.94 1.04 0.98 0.95 0.95 1.57 1.10 1.09

- 10.0 - 12.3 17.5 11.1 43.2 28.1 3.1 -4.9 21.8 21.6 46.2 11.8 70.5 - 125.2 56.9 26.5

Mediterranean Mediterranean Mediterranean Mediterranean

0.992 0.982 0.997 0.980

1.05 1.51 1.23 1.27

-47.1 -82.8 -2.1 - 48.8

Tropical h u m i d - s u m m e r Tropical s e m i - d e s e r t

0.933 0.967

1.37 0.99

-87.3 113.6

a

a P e a r s o n ' s ' correlation coefficient b e t w e e n estimated and m e a s u r e d data. b Slope o f the regression line of e s t i m a t e d versus m e a s u r e d data. c Line intercept.

Sb

i e

SOLAR RADIATION ALGORITHM FOR ECOSYSTEM DYNAMIC MODELS

157

For diffuse solar radiation, the long-term tilt factor (Td) is estimated by assuming that the sky is an isotropic source of diffuse radiation and thus the instantaneous and long-term tilt factors are equal: Tu

cos2(S/2)

=

Finally, the direct and diffuse radiation tilt factors are combined to produce the ratio (Ts) of total radiation received on a tilted surface to that coming in on a horizontal plane:

+ TdRd/Ra

h = Tb(1 - - R d / R a )

The mean monthly solar radiation received on the tilted surface (Rs) is: Rs = R a T ~ ALGORITHM TEST

We have tested our method by comparing the algorithm's estimates of mean monthly solar radiation received on a horizontal surface with re-

Cloudiness (tenths)

S o l a r r a d i a t i o n ( c a l crn-2 d a y - t )

800~_~___~ 700 -

/F 8 /

',,

7

\\\\\

600 F

x\

x

/

5oo L

15

400

300~-

!3

s001

//

0

\\ I

Jan

Feb

I

I

I

Mar Apr May Jun

Jul

i*

I Aug

Sep

Oct

Nov

Dec

Month --

Measured

radiatior~

Estimated

cloudiness

~

Estimated

radiation

Fig. 1. Seasonal patterns of estimated and measured solar radiation and estimated cloudiness for Velen, USSR: 66°10'N, 169°50'E, 7 m above sea level, subarctic tundra climate.

158

N.T. NIKOLOV

A N D K.F. Z E L L E R

ported radiation data for 69 meteorological stations throughout the northern hemisphere (i.e. 38 in the USA, and 30 in Europe and Asia and one in Africa). The climate and radiation data were provided by Miiller (1982) and Zeller (unpublished). Miiller's data were long-term monthly averages. Meteorological stations were selected to represent different climate zones, latitudes and elevations (i.e. from tundra to tropical climate, and from sea level up to 3300 m). We have used three statistical criteria to evaluate the agreement between predicted and reported seasonal patterns of solar radiation: Pearson's correlation coefficient (r), the slope of the regression line (s) of measured versus calculated data, and the line intercept (i). As indicated by r and s values in Tables 1 and 2, model predictions are highly correlated with reported data in almost all cases. Greater deviations from zero of i values typically reflect different lengths of the averaging period of measurement (see Miiller, 1982) for climatic parameters (i.e. temperature, relative humidity and precipitation) versus solar radiation. The algorithm's estimates of the accumulated annual solar radiation deviate from the reported ones on the average by 5.8% for North America (Table 1), and by 11% for Eurasia and Africa (Table 2). Figures 1 through Total r a d i a t i o n

(cal crn-2 day-l)

Cloudiness

(tenths)

800

700

600

5

500

400

300

200

100 L Jan

Feb

I Mar

I Apr

I May

E

I

Jun

Jul

I

~

L

J

A u g Sep Oct N o v D e c

Month

--

Measured

radiation

-~-

Estimated

cloudiness

~

Estimated

radiation

Fig. 2. Seasonal patterns of estimated and measured solar radiation and estimated cloudiness for Irkutsk, USSR: 52°16'N, 104°19'E, 468 m above sea level, continental boreal climate.

SOLAR

RADIATION

ALGORITHM

Solar radiation

FOR ECOSYSTEM

(cal cm-2 day-l)

DYNAMIC

Cloudiness

MODELS

(tenths)

800~

700

159

8

~

600 i

- 7

_

~

\

,-00

~6

\

/

: \i,-

\\

100-

/"

\

,,' '~

// ,~

Aug

Sep

\ 0 I Jan

,

I

Feb

Mar

Apt

I May

~ Jun

Jul

~t I Oct

Nov

iI I0 Dec

Month

--

Measured

radiation

Estimated

cloudiness

~

Estimated

radiation

Fig. 3. Seasonal patterns of estimated and measured solar radiation and estimated cloudiness for Boise (Idaho), USA: 43°34'N, 116°13'W, 866 m above sea level, cool-temperate dry steppe climate.

4 demonstrate estimated and actual seasonal patterns of solar radiation and estimated mean monthly cloudiness for a range of climate zones in the northern hemisphere. CONCLUSION

The verification test showed that our algorithm produces estimates of incoming solar radiation within an accuracy sufficient for ecosystem models with monthly time steps of computation, e.g. gap-phase succession models (Pastor and Post, 1985; Bonan, 1988; Urban, 1990). In fact we have developed this algorithm as a part of a similar model, currently under development at FS Rocky Mountain Forest and Range Experiment Station in Fort Collins, Colorado. The ability of the algorithm to realisticly predict incident solar energy for a range of latitudes and elevations by using only basic topographic and climatic information for input makes it also applicable to ecological studies on landscape and regional scales.

160

N.T. NIKOLOVAND K.F. ZELLER Solar radiation (cal cm-2 day-l)

Cloudiness (tenths)

800

8

700

7

500

I5

J

400

/

~

~

4

300 ~

3

l

200 I

i2

I00 I O~

II I

i

i

J

I

I

i

i

i

L

0

J a n Feb Mar Apr May Jun Jul Aug S e p Oct N o v D e c Monlh --

Measured radiation ~

Estimated radiation

Estlrnaled c l o u d i n e s s

Fig. 4. Seasonal patterns of estimated and measured solar radiation and estimated cloudiness for Miami (Florida), USA: 25°48'N, 80°16'W, 2 m above sea level, tropical humidsummer climate. APPENDIX Listed below is the Turbo Pascal (ver. 4.0 and higher) code of the solar radiation algorithm described in the paper. The code is written as a Turbo Pascal unit that can be easily incorporated in ecological models. For maximum computation efficiency the code has been organized into two procedures. InitSolRad initializes algorithm's parameters that are bare functions of site location (i.e. latitude, elevation, slope, and aspect) and thus need to be computed only once in the beginning of the program using this unit. The procedure SolRad uses these parameters to calculate incoming m e a n monthly solar radiation (array IRad) as a function of varying climate conditions (i.e. m e a n monthly temperature, m e a n monthly relative humidity, and total monthly precipitation) for each year of simulation. The program that uses this unit must do the following in its initialization section: (1) Assign real values to Lat (site latitude, degrees), Alt (site altitude, m above sea level), Incl (slope inclination, degrees), and Aspect (slope aspect, degrees). The format of all angular data is degree.minutes. Latitude is

SOLAR R A D I A T I O N A L G O R I T H M FOR ECOSYSTEM D Y N A M I C M O D E L S

161

positive north of the equator and negative south of the equator. Slope inclination is always positive and in the range 0-90. Slope aspect is also positive and increases clockwise from 0 to 360 where 0 is due North, 180 is due South, etc. (2) Call procedure InitSolRad to compute mean monthly values for potential solar radiation (array PotHorRad), site elevation correction factor (array AltCor), day length (array DayL, h), and direct radiation tilt factor (array Tb) as well as values for the diffuse radiation tilt factor (Td), and atmospheric light attenuation parameters (Arad and Brad). Then, to calculate the incoming mean monthly solar radiation for every year of simulation the program must (1) assign real values to the arrays MTemp (mean monthly temperature, °C), MHum (mean monthly relative humidity, %), and MPrec (total monthly precipitation, ram), and (2) call procedure SolRad. The result is stored in array IRad (incident radiation, cal cm 2 day 1). If necessary, mean monthly solar radiation can be converted into units of W m -2 by using mean monthly day length estimates stored in the array DayL. unit SRad; interface

type monthly - array[l..12]

var Lat,Alt,lncl,Aspect, month

of real;

Arad,Brad,Td

: real;

: byte;

MTemp,MHum,MPrec, PotHorRad,AltCor,Tb,IRad,DayL

procedure InitSolRad; procedure SolRad;

implementation

: monthly;

162

N . T . N I K O L O V A N D K.F. Z E L L E R

function var

ars

arsin

(s

: real)

: real;

: real;

begin if

(s < l) a n d ars

-1)

then begin

:- a r c t a n ( s / S q r t ( l - S q r ( s ) ) ) ;

if s < 0 t h e n end else else

(s >

arsin

:- pi

- ars else

if s - i then a r s i n

if s -

-i t h e n a r s i n

arsin

:- ars;

areos

:- arc;

:- p i * 0 . 5

:- 1.5*pi;

end;

function var

arc

arcos

(c

: real)

: real;

: real;

begin if c - 0 t h e n a r c o s

:- p i * 0 . 5

else

begin arc

:- a r c t a n ( S q r t ( l - S q r ( c ) ) / e ) ;

if c < 0 t h e n arcos

:- pi + are else

end; end;

procedure eonst

InitSolRad;

Sc - 2.0;

Far

sinD,eosD,Ampl,SINa,SINUSa,COSi,L0,1s,Aw,Reor Num,Day

: word;

LastDay

: byte;

: real;

163

SOLAR RADIATION ALGORITHM FOR ECOSYSTEM DYNAMIC MODELS

procedure

ConvertTopoData;

var Sign

: real;

begin L0

:- Lat;

Is

:- Incl;

Aw

:= Aspect;

if LO >- 0 then Sign

:- i else Sign

:- -i;

L0

:- ( I n t ( L O ) + S i g n * F r a c ( L 0 ) * O . 1 6 6 6 6 6 7 ) * p i / 1 8 0 ;

Is

:- ( I n t ( I s ) + F r a c ( I s ) * O . 1 6 6 6 6 6 7 ) * p i / 1 8 0 ;

Aw

:- I n t ( A w ) + F r a c ( A w ) * 0 . 1 6 6 6 6 6 7 ;

Aw

:- (180

- Aw)*pi/180;

if (Is > pi*0.5) writeln('Slope

or (Is < 0) then b e g i n inclination

angle out of range. ');

Halt; end; end;

function

HorRad

var COShs

(var J

: word)

: real;

: real;

begin sinD

:- 0 . 3 9 7 8 5 * s i n ( A . 8 6 8 9 6 1 + O . 0 1 7 2 0 3 * J + O . 0 3 3 4 4 6 * s i n ( 6 . 2 2 4 1 1 1 + O . 0 1 7 2 0 2 * J ) )

cosD

:- S q r t ( l - S q r ( s i n D ) ) ;

COShs HorRad

~- - s i n ( L 0 ) * s i n D / ( c o s D * c o s ( L 0 ) ) ; :- 4 5 8 . 3 7 * S c * ( l + O . 0 3 3 * c o s ( 2 * p i * O . 0 0 2 7 3 9 7 2 6 * J ) ) Sqrt(l-Sqr(COShs))

end;

* (cos(L0)*cosD*

+ arcos(COShs)*lSO/(pi*57.296)*sin(L0)*sinD)

164

N . T . N I K O L O V A N D K.F. Z E L L E R

function

NumOfDays

(var month

: byte)

: word;

begin case month of 2 : NumOfDays

:- 28;

1,3,5,7,8,10,12 4,6,9,11

: NumOfDays

: NumOfDays

:- 31;

:- 30;

end;

end;

function var dl

DayLengthOf

(var Day: word)

: real;

: real;

begin dl

:- Ampl*sin((Day-79)*0.01721)

+ 12;

if dl < 0 then dl :- 0 else if dl > 24 then dl :- 24; DayLengthOf

:- dl;

end;

procedure

SolarAl;

var k : byte; As,h,SNa,CSa,CSi

: real;

begin for k :- i to 24 do begin h

:- (k-12)*15*pi/180;

SNa

:- sin(L0)*sinD

+ cos(L0)*cosD*cos(h);

if SNa > 0 then begin SINa

:- SINa + SNa;

S O L AR R A D I A T I O N A L G O R I T H M FOR E C O S Y S T E M D Y N A M I C M O D E L S

CSa

:I S q r t ( l

- sqr(SNa));

if A w I 0 t h e n CSi else

sin(LO-Is)*sinD+cos(LO-Is)*cosD*eos(h)

:-

begin

As

:= a r s i n ( c o s D * s i n ( - h ) / C S a ) ;

if As > pi

CSi

t h e n As

:= p i

- As;

:- c o s ( A s - A w ) * C S a * s i n ( I s )

end; if CSi > 0 then b e g i n COSi

:- C O S i + C S i ;

SIN-USa

:= S I N U S a

+ SNa;

end; end

else

if k >- 12 t h e n Exit;

end; end;

begin

{"InitSolRad"}

ConvertTopoData; Ampl

:- e x p ( 7 . A 2 + O . 0 4 5 * L 0 * 1 8 0 / p i ) / 3 6 0 0 ;

Rcor

:~ i - 1 . 3 6 1 4 * c o s ( A b s ( L 0 ) ) ;

Arad

:- 3 2 . 9 8 3 5

Brad

:- 0 . 7 1 5

Td Day

- 64.88A*Rcor;

- 0.31831*Rcor~

:- S q r ( c o s ( I s * 0 . 5 ) ) ; := O;

for m o n t h

:- i to 12 do b e g i n

+ SNa*cos(Is);

165

166

N.T. N I K O L O V A N D K.F. Z E L L E R

DayL[month]

:- 0;

PotHorRad[month] SINa

:~ 0;

SINUSa COSi

:~ 0;

:= 0;

:~ 0;

if Is > 0 then Tb[month]

:- 0 else Tb[month]

LastDay

:- NumOfDays(month);

for Num

:- 1 to LastDay do begin

:- i;

Inc(Day); DayL[month]

:- DayL[month]

PotHorRad[month]

+ DayLengthOf(Day);

:- PotHorRad[month]

+ HorRad(Day);

SolarAI; end; SINa

:- SINa/DayL[month];

if SINa > 0 then AltCor[month]

:- (l-exp(-0.014*(Alt-274)/(SINa*274)))

else AltCor[month] DayL[month]

:- 0;

:- DayL[month]/LastDay;

PotHorRad[month]

:- PotHorRad[month]/LastDay;

if Is > 0 then if COSi > 0 then Tb[month]

:- COSi/SINUSa;

end; end;

{ End of "InitSolRad"

function

Cloudiness

var Ev,CI

: real;

(var Temp,Hum,Precip

begin if Preeip

}

> 0 then begin

: real)

: real;

SOLAR RADIATION A L G O R I T H M FOR ECOSYSTEM DYNAMIC MODELS

Ev

:=

CI

:= i0.0

167

Hum*6.1078*exp(17.269*Temp/(Temp+237.3)); - 1.155*Sqrt(Ev/Precip);

if CI < 0 then CI Cloudiness e n d else

:- 0;

:= CI;

Cloudiness

:- 0;

end;

procedure

SolRad;

var ActHorRad,Kt,DiffToHorRatio,Rt,Clouds

: real;

begin for m o n t h Clouds

:- i to 12 do b e g i n :- C l o u d i n e s s

ActHorRad

:- -Arad + P o t H o r R a d [ m o n t h ]

if A c t H o r R a d ActHorRad

< 0 then A c t H o r R a d

:- A c t H o r R a d

if A c t H o r R a d Kt

(MTemp[month],M1~um[month],MPrec[month]);

+ ((PotHorRad[month]+l)-ActHorRad)

then

else D i f f T o H o r R a t i o

:- 1 . 0 0 4 5 + 0 . 0 4 3 5 * K t - 3 . 5 2 2 7 * S q r ( K t ) + 2 . 6 3 1 3 * S q r ( K t ) * K t :- 0.166;

:= ( l - D i f f T o H o r R a t i o ) * T b [ m o n t h ]

IRad[month] e n d else

:- R~*ActHorRad;

IRad[month]

:- 0;

end;

end.

* AltCor[month];

:- ActHorRad/Por_HorRad[month];

DiffToHorRatio

end

:- 0;

> 0 then b e g i n

if Kt < - 0.75

Rt

* (Brad - 0 . 0 3 2 5 9 * C l o u d s ) ;

{ E n d of "SolRad"};

+ Td*DiffToHorRatio;

168

N . T . N I K O L O V A N D K.F. Z E L L E R

REFERENCES Bonan, G.B., 1988. A Simulation Model of Environmental Processes and Vegetation Patterns in Boreal Forests: Test Case Fairbanks, Alaska. Working Paper WP-88-63, International Institute for Applied Systems Analysis, Laxenburg Austria, 63 pp. Bonan, G.B., 1991. Atmosphere-biosphere exchange of carbon dioxide in boreal forests. J. Geophys. Res., 4: 7301-7312. Bristow, K.L. and Campbell, G.S., 1984. On the relationship between incoming solar radiation and daily maximum and minimum temperature. Agric. For. Meteorol., 31: 159-166. Dilley, A.C., 1968. On the computer calculation of vapor pressure and specific humidity gradients from psychometric data. J. Appl. Meteorol., 7: 717-719. Gamier, B.G. and Ohmura, A., 1968. A method for calculating the direct shortwave radiation income of slopes. J. Appl. Meteorol., 7(5): 996-800. Keith, F. and Kreider, J.F., 1978. Principles of Solar Engineering. Hemisphere Publishing Corporation, Washington, DC, 778 pp. Klein, S.A., 1977. Calculation of monthly average insolation on tilted surfaces. Sol. Energy, 19: 325-329. Leemans, R. and Prentice, I.C., 1989. FORSKA, A General Forest Succession Model. Institute of Ecological Botany, Uppsala, Sweden, 45 pp. Lui, B.Y.H. and Jordan, R.C., 1960. The interrelationship and characteristic distribution of direct, diffuse and total solar radiation. Sol. Energy, 4: 1-19. Lui, B.Y.H. and Jordan, R.C., 1962. Daily insolation on surfaces tilted towards the equator. Trans. Am. Soc. Heating, Refrigerating Air Conditioning Eng., 67: 526-541. Lui, B.Y.H. and Jordan, R.C., 1963. The long-term average performance of flat-plate solar engineering collectors. Sol. Energy, 7: 53-74. Miiller, M.J., 1982. Selected Climatic Data for a Global Set of Standard Stations for Vegetation Science. Dr. W. Junk Publishers, The Hague, The Netherlands, 306 pp. Musselman, R.C. (Editor), in press. The Glacier Lakes Ecosystem Experiments Site - GLEES. A Site Description. General Technical Report, USDA/FS Rocky Mountain Forest and Range Experiment Station, Fort Collins, Colorado, USA, 150 pp. Pastor, J. and Post, W.M., 1985. Development of a linked forest productivity-soil process model. ORNL/TM9519, Oak Ridge National Laboratory, Oak Ridge, Tennessee, USA, 140 pp. Running, S.W., 1984. Documentation and Preliminary Validation of H2OTRANS and DAYTRANS, Two Models for Predicting Transpiration and Water Stress in Western Coniferous Forests. Research Paper RM-252, USDA/FS Rocky Mountain Forest and Range Experimental Station, Fort Collins, Colorado, USA, 45 pp. Running, S.W. and Coughlan, J.C., 1988. A general model of forest ecosystem proccesses for regional applications. I. Hydrological balance, canopy gas exchange and primary production processes. Ecol. Modelling, 42: 125-154. Running, S.W., Nemani, R.R. and Hungerford, R.D., 1987. Extrapolation of synoptic meteorological data in mountainous terrain and its use for simulating forest evapotranspiration and photosynthesis. Can. J. For. Res., 17: 472-483. Swift, L.W., 1976. Algorithm for solar radiation on mountain slopes. Water Resour. Res., 12(1): 108-112. Urban, D.L., 1990. A Versatile Model to Simulate Forest Pattern. A User's Guide to ZELIG version 1.0. Environmental Sciences Department, The University of Virginia, Charlottesville, Virginia 22903, 107 pp. Woodward, F.I., 1987. Climate and Plant Distribution. Cambridge University Press, Cambridge, 174 pp.