Computation of diffuse solar radiation

Solar Ener~,,y Vol. 39, No. 6, pp. 521-532, 1987 Printed in the U,S.A.

0038-092X,'87 $3.0(} + .0~) ,~' 1987 Pergamon Journals I.td.

COMPUTATION OF DIFFUSE SOLAR RADIATION C. R. PRASAD Department of Mechanical Engineering and Centre for Atmospheric Sciences. Indian Institute of Science, Bangalore 560012, India a . K. INAMDAR D e p a r t m e n t o f M e c h a n i c a l E n g i n e e r i n g , I n d i a n Institute of S c i e n c e , B a n g a l o r e 560012, I n d i a

and PRABHA V E N K A I ESH

Department of Mechanical Engineering and Centre for Atmospheric Sciences. Indian Institute of Science, Bangalore 560012, India Abstract--A technique for computing the spectral and angular (both the zenith and azimuthal) distribution of the solar energy reaching the surface of earth and any other plane in the atmosphere has been developed. Here the computer code LOWTRAN is used for getting the atmospheric transmitIances in conjunction with two approximate procedures: one based on the Eddington method and the other on van de Hulst's adding method, for solving the equation of radiative transfer to obtain the diffuse radiation in the cloud-free situation. The aerosol scattering phase functions are approximated by the Hyeney-Greenstein functions. When the equation of radiative transfer is solved using the adding method, the azimuthal and zenith angle dependence of the scattered radiation is evaluated, whereas when the Eddington technique is utilized only the total downward flux of scattered solar radiation is obtained. Results of the diffuse and beam components of solar radiation received on surface of earth compare very well with those computed by other methods such as the more exact calculations using spherical harmonics and when atmospheric conditions corresponding to that prevailing locally in a tropical location (as in India) are used as inputs the computed values agree closely with the measured values.

1. INTRODUCTION It is well known that the sun is the primary source of energy for the earth, supplying more than 99.99% of its energy input. Solar radiation absorbed to different extents by the various constituents of the e a r t h - - a t m o s p h e r e , oceans, cryosphere, biomass, and land, e t c . - - i s responsible for the driving forces that control the climate of the world, circulation of atmosphere and oceans, extent of ice cover, ternperature structures of atmosphere and oceans, and so on. It is thus necessary to have a knowledge of solar radiation received at the surface of the earth and the radiation absorbed by the atmosphere, in studies of earth's radiation budget, climate, etc. There are many technological applications where solar energy is utilized directly such as in solar water and air heaters, dryers, solar ponds, photovoltaic devices, solar refrigerators, furnaces, and the like. For evaluating the performance of these devices, a detailed knowledge of the spectral as well as angular distribution of solar energy incident on these energy collecting devices is essential, Extensive solar radiation data sets are being published frequently[l]. But usually these are not completely adequate in that neither continuous observations at all locations nor spectral and/or angular distributions are available. Also, since data is available at only a few fixed locations, the radiation incident at other locations have to be inferred by some extrapolative or empirical techniques. Several such techniques are in use. If, however, data

on spectral distribution is required, it becomes necessary to use an accurate predictive technique for obtaining this quantity. In this article the details of a computational scheme that has been developed to calculate the solar radiation reaching any altitude in atmosphere and the surface of earth are presented. The scheme needs, as data, the profiles of temperature and concentrations of gases and aerosols, in addition to the extraterrestrial distribution of solar energy. The spectral distributions of the direct as well as diffuse components of solar radiation and the angular distribution of the diffuse radiation can be computed as a function of the solar zenith and observer zenith and azimuthal angles. 2. PROPAGATIONOF SOLARRADIATIONIN THE ATMOSPHERE 2.1 Characterization of solar radiation For most purposes, the sun can be considered as a blackbody source at about 5860 K. While the radiation reaching the top of the earth's atmosphere is characterized as above, the solar radiation reaching the surface of the earth is modified considerably during its passage through the atmosphere. It undergoes attenuations due to (1) absorption by gaseous components in certain wavelengths, (2) scattering by gaseous molecules--Rayleigh scattering (which varies as the inverse fourth power of wavelength M, and (3) absorption and scattering by particulate matter present in a t m o s p h e r e - - M i e 521

C.R. PRASADet al.

522

scattering. The radiation thus scattered is distributed in all directions and is referred to as the diffuse radiation. Hence, the radiation reaching the surface of the earth has both a beam (or direct) and diffuse component, whereas the extraterrestrial radiation consists only of the beam component, Figure 1 shows the spectral distribution[2] of the extraterrestrial[3] and the total (beam + diffuse) radiation reaching the surface. The total path length that the solar radiation has to travel through the atmosphere increases as the solar zenith angle increases, resulting in increased attenuation of the radiation. Thus, at noon when the zenith angle is at its minimum, the attenuation is a minimum, whereas at sunrise and sunset the attenuation is a maximum. 2.2 Earth's atmosphere The earth's atmosphere extends up to about 100 kin. It consists of[4] mixed gases like N2, 02, CO, CO2, CH4, N20, and so on, whose relative concentrations do not vary with height, and unmixed gases like H20, 03, and HNO3, whose concentrations are a function of height. As these unmixed species have strong absorption bands within the solar spectrum, they exert a varying influence on the incoming solar radiation. In a cloud-free situation, aerosols play a very crucial role in modifying the solar radiation reaching the surface of the earth, as they strongly scatter and to some extent absorb the solar radiation. The relative proportion of direct-to-diffuse radiation is primarily controlled by

the aerosols. Aerosols--both natural and man m a d e - - v a r y in a random manner geographically and temporally. Their size and number density distributions vary considerably depending on their origin[5]. When water precipitates like rain and fog are also present the situation is further complicated. In the presence of clouds, the beam component is drastically reduced and depending on the density and extent of clouds it may only be the diffuse radiation that reaches the ground. Here the atmosphere is modeled so that the mean climatological conditions are normally used as the inputs in the default case. However, when actual meteorological conditions are known they can be used in place of these mean condition as inputs. 2.3 Atmospheric model The atmospheric model chosen for our computations is the one used in the LOWTRAN code[6] since they are representative of the average conditions. Prasad and Venkatesh found, in a recent study[7], that the mean temperature and water vapor profiles obtained from averaged aerological data over India for the years 1975-1981 agreed fairly well with the tropical profiles used in LOWTRAN code (Fig. 2). To account for the great variation in the composition, concentration, and sizes of aerosols, the

lOO

\~

\

90

200o

80

\

7c 1600

\/ ~ \\ ~ '

~E 6C

~

EXTRATERRESTRIAL RADIATION

OZONE

~

~ 50 ~:

, 4o

= 0o 60 o

z<

30

800

TEMPE

20 WATER VAPOR

-10

400 t l / "

0 0.25

IIl/I/~

0.85 1.45 WAVELENGTH (~tm)

I 200

2.05

Fig. 1. Extraterrestrial and total (beam + diffuse) radiation for two different incident directions computed using the model developed here: Tropical model with rural aerosols and 23 km visibility.

240 240 2'60 TEMPERATURE ~3K )

X

i & 280

300

[

I

I

I

I

I

I

10- 4

10-3

10 2

10-1

100

101

102

103

I 10-6

PRESSURE( mb ) I I I ] 10- 4 10- 2 100 102 WATER VAPOR (gm/m 3 )

I 104

I

10 10

I 10-8 I

I

I

10-11 10 10 10-9

I

I

I

I

I

10-8

10-7

10-6

10-5

10-4

OZONE(gm/m 3 )

Fig. 2. Profiles for tropical atmospheric model.

I

Diffuse solar radiation atmosphere has been divided into four separate regions: the boundary layer (0-2 km), the troposphere (0-10 km), the stratosphere (10-30 km), and the mesosphere (30-100 km). In the boundary layer itself, there are three different types of aerosols: rural, urban, and maritime. They are defined and presumed to have a bimodal size distribution[5]; their compositions as well as refractive indices are specified depending on where they originate from. The appropriate distributions and number densities of the aerosols are automatically chosen once the boundary layer visibility is specified. The effect of the relative humidity on the particle size and refractive index is also included. Tropospheric aerosols are assumed to have the same composition as the rural aerosol except that large particles (diameters greater than 1 pxm) are absent. Five different situations are considered for the stratospheric aerosols: background, moderate, and high aged volcanic and moderate and high fresh volcanic[6]. 2.4

523 SOLAR RADIATION

\

8 ~ o ~ / /

\

\

/

/

/

~

TOP

~E R E ~)

1 l x

d× | I

SCA~ RADIATION

OPTICAL DEPTH 7- =/oX/~(x) dx

\\ \ \ 8

Equation and radiative transfer

Computation of diffuse solar radiation in the atmosphere has been attempted by many researchers. These attempts range from the direct solution of the equation of radiative transfer describing the propagation of solar energy through the atmosphere to many other simplified procedures. Simplifications become essential in any attempt to compute the diffuse radiation except for a few very simple conditions in view of the complications introduced by single and multiple scattering by particulates, The equation of radiative transfer (ERT) that expresses the radiative energy balance in the medium can be expressed[8] as follows:

Fig. 3. Schematic diagram showing the solar zenith and observer zenith angles.

transmitted intensity) is obtained[8] by integrating (1) as

I(T, dl(T, IX, ¢)/d~ = I(~, Ix, +) -

-ix, 4)) = I(0, -Ix, ¢) e x p ( - "r/tx)

[too/4"rr]~vFo -(too/4w)

x POx, ¢; - DXo,qbo) e x p ( - dlxo) -

[too/4~r]

r 30

d+'

r ' [d~x' f -1

X I(-r, ~x', (b')P(P, Ix', d~, ¢b')]

× (1)

for transport of solar radiation in the atmosphere. l(v, ix, ¢) denotes the intensity at optical depth ~ in the direction specified by the zenith angle O (O = c o s - 1 Ix) and azimuth ~b (see Fig. 3). Here, I refers only to the diffusely transmitted or reflected quantity and the beam component of the solar radiation is found separately by Beer's law from wFo, the solar irradiance falling at the top of the atmosphere (~ = 0) on a unit area normal to the direction ( - Ixo, 60), as (¼) Fo e x p ( - djXo) The sign convention used in this analysis is: a positive Ix for upgoing radiation (0 < 0 < ~r/2) and a negative ix for downcoming radiation (~v/2 < O < v). ~Oo is the single scatter albedo defined as the ratio of the scattering optical depth to the total optical depth. A general expression (for downward

:

f~ dr' f~'~ dqb' f j, +;

tiT':

+')

× exp[--(v -- ~')/IX]/Ix] -- (eoo/4~r) / cT × t'rrF° Jo d'r'P(r'; Ix, ¢b; --P,o, 4)o) × exp(-T'/#o) x exp[-(~ -

"r')llx]llx)

(2)

/

In the above solution (2), the quantities that need to be specified are the optical depth T and the phase function P. The optical depth has to be computed along the path of the ray from a knowledge of the amount of absorbers in the path and their scattering cross sections. When there is no scattering present, solution (2) simplifies to Beer's law. But, except under some special situations like isotropic scattering, single scattering, etc., closed-form solutions of (2) are not possible. In hazy atmospheres, where multiple scattering occurs, the solution of eqn (2) poses formidable difficulties unless some simplifications are made.

524

C.R. PRASADet al.

The phase function, which describes the angular distribution of scattered energy, depends on the scattering properties of the gas molecules and particulates in the atmosphere. When the wavelength of the incident light is much larger than the particle size, as in the case of gas molecules (i.e., Rayleigh scattering) the phase function for incident unpolarized light is given by[8] P(cos 0) = ~(1 + cos 2 ~)

(3)

where + is the scattering angle. For this case, simple solutions to the ERT can be obtained and the scattering characteristics of a pure Rayleigh atmosphere have been studied in detail[9], The computation of phase functions for scattering by particles whose size is of the order of or larger than the wavelength of incident radiation involves the use of the Mie scattering computation codes[10, 11] derived from the classical electromagnetic theory. Instead of undertaking Mie scattering computations for each case, parametric representations of phase functions can be employed, One such is the Hyeney-Greenstein phase function[12]. By using a parametric form of phase function, solutions to the ERT become fairly simple and are of sufficient accuracy in most cases. Throughout this work, the Hyeney-Greenstein function has been used, which contains an asymmetry parameter g, defined as the first moment of the phase function, This phase function is mathematically well behaved and can be made to represent very different scattering conditions quite well by varying the value of g. 2.5 Approximate solutions There are many approximate solutions that give solar radiation integrated over the full angular range of the observer angles. These include the twostream[13] andmultistreamapproximations[14, 15]. F o r the general case of nonisotropic aerosol scattering, generalized phase functions like the Hyeney-Greenstein function ( H - G function) have been used by several authors with good success in estimating the diffuse radiation in the galaxy[12] and in our atmosphere[16]. In the case of Eddington approximation[14], the intensity is expressed as a series wherein only the first two terms are taken: 1(% tx) = 10(7) + txlj('r)

(4)

When this is substituted in (1), terms of order higher than l will vanish. The phase function is also expanded in a series and taking the first two terms P(cos 0) = 1 + oJ~

(5)

In other two-stream approximations, the H - G function also provides a very effective simplification for the analysis of anisotropic scattering in the atmosphere[15]. In a modified two-stream model[13], the

dependence of the reflectivity and transmission of the atmosphere on the angle off the incident radiation and on the angular dependence of the scattering phase function of the medium have been specifically included. This model[13] has been employed tO assess the effects of regional aerosol layers on planetary albedo. But the model is applicable only for optically thin plane parallel atmospheric layers with clear-sky conditions. Although all the above methods are easy to compute numerically, their greatest drawback lies in the fact that they do not yield the detailed angular distribution of diffuse radiation, but only their integrated quantities. Another limitation of the Eddington approximation is that it fails for very small optical depths and large zenith angles. In fact, the two-stream approximation[15] gives even negative values of transmission ['or low total optical depths. The Eddington approximation[14] has been chosen here as one of the methods for solution of ERT since the computations involved are very simple and hence provides a means for rapidly calculating the irradiances within as well as those emerging from the boundaries of the atmosphere. 2.6 Numerical methods" Many techniques have been developed in the recent years for solving the equation of radiative transfer through terrestrial atmosphere for the evaluation of the diffuse radiation. One such technique for handling multiple scattering calculations in a plane parallel atmosphere, is the adding or doubling method outlined by van de Hulst[17] as far back as 1963, based on simple geometrical ray tracing and involving the doubling of very thin layers[18]. Basically, in this method, two individual, very thin layers are taken first, with the reflectance and transmittance of the combined layer found next. In fact, many of the methods known as the matrix operator method, star-product method, etc. are nothing more than variations of this procedure. The adding method has been employed by many authors[17-19] and yields accurate solutions to ERT without the need to fully solve it. Another approach is the iterative technique[20, 21], wherein the atmosphere is divided into a number of layers, each with optical thickness ranging from 0.02-0.05. The single-scatter solution is then formulated and the value of the diffuse intensity at the end of the first layer is found. This intensity distribution is used to calculate the value for the second, third layers, and so on. These values are then substituted into the original full transfer equation (1) to obtain the next set of values and a G a u s s Siedel iterative scheme is used. The integrations over azimuth angle, zenith angle, and optical depth are carried out numerically. Such an iterative technique is too slow to be of great value particularly for large optical depths. In yet another method, for multiple scattering calculations in plane-parallel atmosphere, the scat-

Diffuse solar radiation

525

tering phase function is expanded[22] in terms of a Fourier series as,

an excessively large number of computations are performed.

P[P., Ix', (do - do')] x = ~7; F,,(IX, Ix') cos[(n - 1) (6 - 6')] , 1

2.7 Solution of the equation of transfer Our interest in solving the ERT arises primarily from the desire to compute the spectral and directional characteristics of the solar beam and diffuse radiation reaching the earth. Of the many possible methods, we have chosen van de Hulst's adding method[17] for getting the angular distribution of diffuse radiation, despite the fact that there are other methods that are more rigorous and accurate. This is because of the following reasons: the adding method is simple and needs modest computational resources even when the optical depths are larger and the result is accurate enough especially since the uncertainties in the input quantities like the meteorological parameters, aerosol properties--like number density, size distribution, and optical properties--are themselves large enough as to not to warrant the use of a very accurate computational procedure[27]. The Eddington approximation[14] is also separately included as an alternative method of solution, which can be used whenever only total diffuse radiation is required. These two programs are integrated with the well-known comprehensive program LOWTRAN 5[6], which computes transmittance of a real atmosphere. There are a few other programs reported in the literature for computing the solar spectral direct and diffuse irradiances. One of these is SOLTRAN[28], which uses a Monte Carlo technique for taking care of multiple scattering. In LOWTRAN 6129], a single scattering scheme has been introduced for taking the solar or linear scattered radiation into account. The work detailed in this article was carried out independently[30] and was used for several other solar energy applications[2, 30].

(6)

It is shown that[22], for a desired accuracy, the seties can be terminated with fewer and fewer terms as the direction of incidence (or of scattering) is moved toward the zenith. An iterative scheme is then used for computing the scattered radiation up to the nth approximation, which amounts to considering up to nth order of scattering. The iterations, however, are cumbersome, requiring evaluation of a number of integrals at each step and a large number of iterations are needed for the solution to converge even for a homogeneous atmosphere. This technique has been utilized[23] for analyzing the effect of aerosols on the absorbed, reflected, and transmitted solar energy in cloudless, nonhomogeneous, plane-parallel realistic atmospheric models, In C h a n d r a s e k h a r ' s discrete ordinate method[24], the solution of ERT is explicity derived by employing a finite set of discrete streams, representing emergent angles, in the integral form. The phase fimction is expanded in Legendre polynomials: ,v P(cos +) = ~ mkPk(cos tb) k o

(7)

where mkare a set o f N + 1 constants and~oo = 1. The intensity may also be expanded in the form ,~, I(~; ~x, do) = ~ tn

Im(~, p~) COS rn(cb - doo)

(8)

o

And by letting

3. ANALYSIS OF RADIATIVE TRANSFER IN THE ATMOSPHERE

iv

Im(~, IX) = ~ la(T) Pk(t*), k

0

(9)

the ERT is reduced to a set of simultaneous firstorder differential equations, which can then be solved by conventional methods. Dave[25] employed this method, also known as the spherical harmonics technique, to obtain the diffuse solar flux in the earth's atmosphere. The Monte Carlo techniques[26] for multiple scattering have also been used for solving the ERT. Solutions for a plane parallel inhomogeneous Mie scattering atmosphere without recourse to any approximation to the phase function and/or vertical nonhomogeneity have been obtained[26], This method being statistical gives inevitable fluctuations and yields results of rather poor quality unless

For scattering atmospheres, the diffuse up and downward flux densities for any given v are given, respectively, by

Fd " (t) =

l(x; ~x, 6)Ix dtx ddo

(10)

Jo l(v; -tx, do) tx dtx ddo

(11)

(.2w (- 1

Fd ~, ('r) = )o

where 1(~; Ix, do) is the intensity of diffuse radiation given by eqn (2). For axJmuthally averaged intensities, the azimuth-independent phase function can be defined as

P(IX, ~') = I (>" ~

j.

P(ix, do; tx', do') ddo'

(12)

C.R. PRASADet al.

526

3.1 Single-scatter solution The source function for radiation that has been scattered only once is

The scattering angle ~ is given by cos + = be be' + [(1 - be2) x (1 - be,2)]l/2 cos(+ - + ' )

(13) J(*; be, ~b) = (¢Oo/4W)~rFoP(~x, 6;

The phase function P(be, ~b; ~', 6 ' ) is generally not available in closed form. Based on the spherical geometry involved in the scattering process, it is found convenient to expand the phase function in Legendre polynomials with a finite number of terms, N as in eqn (7). Using the addition theorem for spherical harmonics, this equation can be expanded for the argument shown in (13) as

- ~xo, (bo) exp(-,/beo)

It is assumed that there is no diffuse radiation incident in the downward direction at top nor any upward radiation at the base of the finite atmosphere (i.e., a perfectly absorbing ground), which gives the boundary conditions as: 6) = 0 I(Tj; be, +) = 0

I(0; - ~ ,

P(be' 6; be" 6 ' ) N

= ~,

(18)

(19)

N

~

rn--O k~rrl

to'~nP~(be) P~(be') cos m($ - 6 ' ) (14)

where P;Y are the associated Legendre polynomi-

The reflected and transmitted intensities for a finite layer of optical depth v~ are given by 1(0; be, 6) = [~oobeoFo/4(be + beo)] P(P-, 6;

als[31] and

-beo, Oo) {1 - e x p [ - ,,(l/be + l/beo)]} (20)

to~ = (2 - 8~")(k - m)!/(k + m)!

and

with

I(TI: - be, +) = [~o0beoFo/4(~ - beo)] 8~~ = 1 f o r m = 0,

x P(-be, (b; -beo,~bo) [(exp(-vl/be)

= 0 otherwise

- exp(-T~/beo)]

for ~t # be0

= [o~0*lFo/4~to] P ( - b e , cb; -beo, +o) A great deal of simplification occurs if an analyrical approximate expression is used for the phase function. One such approximation is the one suggested by Hyeney and Greenstein[12].

x exp(-¢l/be0)

for be = beo (21)

Now the reflection and transmission functions can be defined as follows:

PM(COS O) = (1 -- g2)/(1 + g2 _ 2g COS ~)B/2

(15) where g is an asymmetry factor representing the cosine-weighted average of the exact phase function. Subscript M indicates phase function for aerosols (or Mie scattering). The values of g for many atmospheric aerosols have been found from exact Mie theory computations by Shettle and Fenn[5]. The Legendre coefficients (see eqn (7)) for the H e y n e y - G r e e n s t e i n function have a simple form and are given by: cok = (2k + 1)g * For Rayleigh scattering: PR(COS 6) = 43(1 + COSz 6)

(16)

If VS,R denotes the Rayleigh optical depth and • S.M denotes the aerosol optical depth, the combined phase function P(cos tb), can be written as P(cos +) = [*s,uPR(COS +) + Ts,MPM(COS +)]/(*Sm + *S,M)

(17)

R(~, ~b; beo, +o) = I(0, ~, +)/(beoFo)

(22a)

T(~x, qb; beo, +o) = I(v~, -~x, +)/(p~oFo)

(22b)

3.2 Adding method As multiple scattering occurs for many of the situations encountered in the earth's atmosphere, the single scattering analysis will not be valid. Under such circumstances we use the van de Hulst's adding method wherein the atmosphere is divided into thin plane layers such that the optical depth of each of the layers is very small--in each such layers only single scattering occurs. After obraining the reflectivity and transmissivity of each such optically thin layer, composite layers of greater optical thickness are built up by "adding" the properties of the thin layers. Basically, for adding, a straightforward geometric ray-tracing technique is used. Let 1 and 2 denote two adjacent plane parallel layers of the atmosphere, beo and be ( = cos Oo and cos O) the incident and emergent directions, R~ and l'z the reflection and total (direct plus diffuse) transmission functions for the first layer and R2 and T2 for the second layer. If the same layers are illuminated from below, these functions are, in general, different from those due

Diffuse solar radiation

527

to illumination from above and are denoted by R * , TT, R~ and T*. L e t the dombined total transmission and reflection functions at the interface of the layers 1 and 2 be D and U. The set of simul-

sity term is e x p a n d e d as in eqn (4), substituted in eqn (1), and the integrated o v e r + and ix giving

taneous equations governing the diffuse reflection and transmission functions for the two layers can now be written as follows:

- (too/4Fo)(l - 3gtxp~o) exp(-~/txo) (27)

(1o + txl~) - ~oo(1o + gpd~)

U p o n integrating eqn (27) o v e r ix, and taking the second m o m e n t of(27), obtained by multiplying the

D = T~ + ST~ + S exp(-T1/tXo) U = R 2 D + R2 e x p ( - ~ / t X o )

R12 = RI + e x p ( - ~ l / I x ) U + T * U

d u o + Ixl~)/d~ -

same by ~x and integrating, a pair of ordinary differential equations is obtained as follows: (23) d I i / d ~ = 311 - too('r)] Io

T12 = e x p ( - r2/Ix)D

- 4aCOo(r)Fo exp(-~/txo)

(28a)

+ T2 e x p ( - ' r j / t x o ) + T z D dlo/d'r = [1 - tOo('r)g("r)] 11

w h e r e S = R T R 2 ( 1 - R ~ R 2 ) -1 By analogy, expressions for R*2 and T*2 for illumination from below can be written down simply by inspection as:

+ ~o~o(T)g(v)p~oFoexp(-~/txo)

(28b)

S = RzR~

These are solved subject to the radiation boundary conditions given by eqn (19). F o r an atmosphere consisting of N layers, the solutions for the ith layer

D -- R* exp(-~2/tXo) + R * U

can be written as:

U = T* + S e x p ( - ' r z / g o ) + S T *

(24)

lo.,.('r) = C~,,. e x p ( - k i ' r ) + C2.~ exp(kiT)

R*2 = R* + e x p ( - ' r 2 / ~ ) D + T2D

- e~ exp(-'ri/ix~)

T~2 = e x p ( - v l / i x ) U

(29a)

ll,,(~) = Pi[Ci,i e x p ( - k i ' ~ ) - C-_.i exp

+ T~ exp(-v2/IXo) + T * U

x (k~T)] - f3i exp(-~J~xo)

(32b)

The quantities S, U, and D in (24) are not necessarily the same as those in (23). In the foregoing equations, the product of any two parameters implies integration o v e r the adjoining solid angle. Thus, for the azimuth independent

where

case,

c~i = 3CooiFoix2[1 + gi(1 - tOoi)]/[4(1 - ki2 ~.L~)]' 1

R*Rz

= 2

f0

p, = [3(1 - ~Oo~)/(1 - too~g~)]v2

[3~ = 3too~Folxo[1 + 3g~(1 - too~)tXo]/[4(1 - k~txo)]

~

R~ (ix, Ix')R2(IX , IXo) Ix'dlx'

k~ = [3(1 - too~) (1 - too~gi)] v2

(25) ro = 0,

Azimuthal variation can be incorporated by expanding these functions in a Fourier Series with (+ +o) as the argument -

and

r,~, is the total optical depth.

The constants C~,~ and C,~ are determined from the two boundary conditions for zero flux at top and b o t t o m (see eqn 19) and the continuity of intensities at each layer i (i.e., flux)

T(ix, IXo; + - +o) = To(p, ixo)/2 +

~

Tm(tx, tXo) cos m(cb - (bo)

(26)

F j = 2'rr

(I~ + txl~)lxdix = ~(1o + ~1~)

(30)

rn--I

where T,,(ix, IXo) = (l/~r) fo ~ T(ix, txo; ~5 - +o) cos

where IX < 0 c o r r e s p o n d s t o F j

re(d0' - +o) dqb'.

Thus

A similar relation also holds for Rm. It is seen that Rm( . . . . ) and T,,( . . . . ) also double and add independently like R and T (m = 0 corresponds to the azimuth independent case), 3.3 E d d i n g t o n a p p r o x i m a t i o n As m e n t i o n e d earlier, w h e n one is interested only in the total flux of scattered radiaton in the d o w n w a r d direction and not its angular dependence, a quick and simple m e t h o d of obtaining it is to use Eddington approximation. In the Eddington approximation[14], the inten-

+ (r) and vice versa.

Fj(0) = 7r[lo(0) - ~lr(0)] = 0

(31a)

Fa('r) = v[lo('r) + ~l~('r)] = 0 lo,i('ri) = Io,~+ t(ri),

i = I, 2 . . . . .

(N -

I~,~(~) = I~,~+~(~),

i = 1, 2 . . . . .

(N -

(31b) 1) (31 c) 1) (31d)

3.4 M e t h o d o f c o m p u t a t i o n s As we m e n t i o n e d earlier, the L O W T R A N 5 c o m p u t e r code d e v e l o p e d at the A F G L laboratory

528

C.R. PRASADet al.

is a very comprehensive computational scheme for obtaining the transmittance and thermal radiance (i.e., in the absence of any imposed external radiation field) along any slant path in any specified atmosphere. But when an external radiation field is present, or when the atmosphere is not in local thermodynamic equilibrium, the radiance computation as given in L O W T R A N 5 cannot be used, since the source function in that scheme has been taken to be merely the Planck blackbody function at the local temperature, Here the computational method that has been developed starts with the solution of ERT by either the Eddington approximation or the adding method and then uses the LOWTRAN 5 computer code only to calculate the transmittances that are subsequently utilized in the solutions of the ERT to derive the radiances when solar radiation is present, Thus the computer code evolved here consists of a modified subroutine TRANS of LOWTRAN along with a few additional subroutines such as DUBL, P H A S E , PRAL, SS, ADD, etc. These subroutines accept the values of optical depth, single-scatter albedo, and asymmetry parameter for each layer of the atmosphere generated in TRANS and compute the spectral radiances for any specified incident zenith and observer zenith angles. The data for the extraterrestrial solar spectrum based on the standard NASA/ASTM[3] curve and also the asymmetry factors for the tropical atmosphere[5] have been included in the subroutine TRANS. The additional subroutines developed here compute the reflection (R) and transmission (T) functions (eqns 22) for a multilayered atmosphere starting from the single-scatter solution for very thin layers, The computations are first started for each layer using the single scatter expressions (modified as below) for an optical depth less than 10 -5 (~ 2 ~7). This optimum value has been arrived at after striking a balance between the computational time and accuracy. Any smaller starting optical depth than this does not improve the accuracy more than 1%. When v ~ 1, the eqns (20) and (21) can be rewritten for R and Tfunctions by expanding the exponentials in a power series and retaining terms only up to v2, for example, R(tx, P~0) = [~oo'r/41XlXo][1 - ~o/2 x (l/tx + 1/IX~)]P(Ix, ~o)

(32)

T(u, uo) = [~oo~/4~xp.o][1 - To/2

x (1/tx + l/tx0)] POx, txo) Beginning at top of the atmosphere and depending on the optical thicknesses obtained from the subroutine TRANS, the layers of the L O W T R A N model are thus further subdivided in the subroutine DUBL, until a layer of optical thickness less than 10 -5 is obtained. The single scatter eqn (32) are then evaluated by the subroutine SS, and returned

to DUBL. The azimuthally averaged phase functions used in these equations are obtained from the subroutines P H A S E and PRAL. The H y e n e y Greenstein function is evaluated using the Legendre expansion (eqn 14) up to 40 terms. The reflection and transmission functions thus obtained for each layer by doubling are then added (subroutine ADD), beginning at the topmost (first) layer, adding with the second layer and then adding the combined first and second layer with the third and so on. When computations were carried out tbr a number of wavelengths for the azimuth-dependent case, up to 20 terms were tbund to be sufficient for convergence of the series. For Eddington's approximation, the relations (29) are evaluated using the boundary conditions (31) and a matrix inversion algorithm for the solution of simultaneous equations thus obtained, in subroutine EDD. Here it was found to be sufficient to divide the atmosphere into only four layers: 10010, 10-2, 2-1 and 1-0 km. This has been done to ensure that the optical thickness of each layer lies within the valid range of the Eddington's approximarion. The phase function and single-scatter albedos (weighted with optical depth) are averaged (subroutine EDD4) over these layers. However, it has been found that this does not introduce any significant errors in the consideration of the inhomogeneity of the atmosphere, as proved by the results.

4. COMPUTATIONSAND RESULTS Computations of spectral, diffuse solar radiation in cloud-free situations have been carried out for tropical and midlatitude summer models of the atmosphere using the procedures described in the previous section. The meteorological range or visibility was taken to be 23 km corresponding to a clear atmosphere and 10 km for hazy conditions. A boundary layer with rural aerosols and stratosphere with background aerosols were assumed. Eighty-eight unequally spaced wavelengths in the solar spectrum (between 0.29-2.5 ixm) were chosen for the computations. Values of the extraterrestrial solar flux at these wavelengths were obtained from the ASTM standard solar spectral curve[3]. Figures 4 to 6 show the azimuthally averaged intensities as a function of the observer zenith angle for a few selected wavelengths, for the tropical, rural aerosol model with 23 km visibility. Results presented are for observations at ground level with a perfectly absorbing ground (i.e., ground albedo is zero). In these figures, solar zenith angle is 0 °, 60 °, and 84.261 ° from the zenith. Strong forward lobes due to aerosol scattering result in the intensity peak around the direction of incidence. The strength of this maximum decreases with an increase in wavelength with this peak disappearing altogether for longer wavelengths. This can be attributed to decrease in Rayleigh scattering (k -4)

Diffuse solar radiation

40

480 L

400

~

TROPICALATMOSPHERE RURAL AEROSOLVlS :23km INCIDENTZENITHANGLE:O°

~

/ / TROPICALATMOSPHERE RURAL AEROSOLV I S : 23 krn INCIDENTZENITHANGLE: 84.261°

r. ~ 32

©

-

320

529

\

//

E

:~ ~ 24

/

~ 24C

/

/

5

/

a: 16

4[

~

/~-- 0.615

1 ZENITHANGLEOF OBSERVATION(DEG) 20

40

60

80

Fig. 6. Azimuth averaged diffuse solar radiation reaching

ZENITH ANGLE OF OBSERVATION ( DEC )

the earth shown as a function of observer zenith angle.

Fig. 4. Azimuth averaged diffuse solar radiation reaching the earth shown as a function of observer zenith angle. Solar zenith angle is 0 ° (directly overhead).

Solar zenith angle is 84.261 ° . The peak in the curve for 0.535 and 0.615 ixm are displayed from that for l.l Ixm. the

and also due to a very rapid increase in intensity that occurs at observer angles close to the horizon as the air mass increases sharply. Furthermore, the w a v e l e n g t h s 0.9935, 1.578 p,m (Figs. 4 and 5) lie in

regions

of

atmospheric

absorption,

bands.

MIDLATITUDESUMMERATMOSPHERE RURAL AEROSOL VlS: 50 km 80= ~ = 600

~ ~ io a

weak

whereas the wavelengths 0.93 and 1.1 Ixm (Figs. 5 and 6) are located at the water vapor absorption

~/+

,'~-~\~ i

TROPICALATMOSPHERE RURAL AEROSOLV l S :23 km INCIDENTZENITHANGLE: 80o

'~ E ..=

Ul°

~o

~\ k~,

s ~ .

\\\~ "X.\~

L

LO

\\

10

20 40 80 80 ZENITHANGLEOF OBSERVATION(DEC)

Fig. 5. Azimuth averaged diffuse solar radiation reaching the earth shown as a function of observer zenith angle, Solar zenith angle is 60 °. Note that logarthmic scale has been used on the ordinate,

0

4'0

80 120 AZIMUTH ANGLE (DEC)

160

Fig. 7. Intensity of diffuse solar radiation as a function of the observer azimuth angle reckoned from the solar meridian for midlatitude summer model. Dave's results (dotted lines) shown for comparison.

530

C . R . PRASAD et al.

r| _

MIDLATITUDE SUMMER ATMOSPHERE RURAL AEROSOL VIS: 50 km = 0.615 #m : 45 °

/

/~\

/

\

103 \ /

//

.%,

=90 ° /

\~

10

.~ -/

~

10'

.

- 0.5350um =

\\

..

//

MIDLATITUDE SUMMER ATMOSPHERE RURAL AEROSOL VIS: 15km

~

/

DAVE

--

~ --

~

\

20 310 410 ;o ;o ~o ~o ZENITH ANGLE OF OBSERVATION (DEG)

~'o j

Fig. 8. Diffuse solar radiation intensity shown as a function of o b s e r v e r zenith angle for midlatitude s u m m e r model, The observer azimuth angles are used as a parameter, D a v e ' s results (dotted lines) shown for comparison,

i 10

I I I 1 I 170 180 20 30 40 50 60 ZENITH ANGLE OF OBSERVATION (DEG)

/k/ \ ~

TROPICAL ATMOSPHERE RURAL AEROSOL VPS: 23km ZENITH ANGLE : 0°

\

1500

~t~ ~" 1000

'~ATERRESTRIAL

<:

.EDDiNGTON

\

d£ -~ 500

0-0

0.4

0.8

!

~o

Fig. 9. Diffuse solar radiation intensity as a function of observer zenith angle for midlatitude s u m m e r model s h o w n for one azimuth angle. Visibility is 15 km. Dotted lines show D a v e ' s results.

^

2000

~

0

1.2 1.6 WAVELENGTH (prn)

2.0

2.4

Fig. 10. Graph s h o w s the good agreement between the adding m e t h o d and E d d i n g t o n ' s method.

Diffuse solar radiation Computations have also been carried out for the azimuth dependent intensities by considering the Fourier expansion, eqn (24) of the reflection and transmission functions. Results presented in Figs. 4 to 6 were for the azimuth independent case wherein only the first term of the Fourier series is taken. However, to obtain the azimuth dependence of the intensities, higher-order terms have to be ineluded and it was found that the series converged after including about 15 to 20 terms. The variation of diffuse radiation intensity with azimuth for the sun's and observer's zenith position of (0 = 0 = 60°) is shown in Fig. 7 for three different wavelengths. Figure 8 shows the diffuse intensity as a function of the zenith angle of observation for different azimuth angles (referred to the sun's meridian) for 0.615 pum wavelength. The results shown in Fig. 9 are for 0.535 p~m wavelength for intensities in the sun's meridial plane with sun overhead (0 = 0°). The dotted lines in these figures (7-9) show Dave's[32] results for comparison, who has carried out the computations using the spherical harmonics techniques and exact Mie phase functions. Since the aerosol models used by Dave are slightly different from those used here, the visibilities V were chosen so as to give aerosol optical depths matching with those of Dave[32]. We see from these figures that there is very good agreement between our cornputation and Dave's results. Figure 10 shows the fluxes computed both by adding and Eddington methods when the sun is directly overhead, The CPU time taken for the method varied between 45 to 60 s for one wavelength for the azimuth independent case on DEC 1090 system. For the azimuth dependent intensities, it took about 2.5 min when 20 Fourier coefficients were used[30].

5. CONCLUSIONS In the work presented here, a fairly simple and effective method has been developed to compute the spectral and angular distributions of both the beam and diffuse components of solar radiation reaching the surface of the earth and any other location in the atmosphere for a given set of local meteorological conditions. The climatological mean profiles of temperature, humidity, and ozone are available as default values if the local meteorological conditions are not known. Here the ground has been taken to be a perfect absorber, although it is possible to use any finite value of its reflectivity. Hansen[18, 33] demonstrated that the use of the H y e n e y - G r e e n s t e i n phase function approximation gets to be more and more accurate with increasing optical depths. Here we have used the H y e n e y Greenstein phase function even for moderate optical depths and the results are seen to be quite good. The agreement of our results with those of Dave are very good although some discrepancies can be observed near the peak of diffuse scattering.

531

It may be recalled that the aerosol models that have been used by Dave[32] and in our analysis are not the same and the phase function used by Dave was rigorously computed. An advantage of the adding method is that it is sufficient to calculate the transmission and reflection matrices only once, since the same values can be stored and used over and again for any combination of zenith and azimuth angles. Furthermore, the adding method is computationally very stable and is uniformly valid for any range of optical thickness[33]. The subroutines developed in our work can be directly incorporated into the LOWTRAN[6] scheme to include the effect of an external source function (the sun's radiation in this case). Acknowledgments--The authors acknowledge the sup-

port of the Indian Space Research Organisation for some of the computations presented here.

REFWREnCES 1. A. Mani, Handbook gfRadiation Data for India. AIlied Publishers Pvt. Ltd., New Delhi, India (1980). 2. Prabha Venkatesh and C. R. Prasad, Computation of solar spectral irradiance. Proc. National Solar Energy Convention 1981, Paper # IN:04. Allied Publishers Limited, New Delhi (1982). 3. M. P. Thekaekara, Solar energy outside the earth's atmosphere. Solar Energy 14, 109 (1973). 4. S. L. Valley, Handbook o f Geophysics and SpaceEnvironments. AFCRL, Bedford, MA (1975). 5. E. P. Shettle and R. W. Fenn, Models of the aerosols of the lower atmosphere and the effects of humidity variations on their optical properties. Report AFGLTR-79-0214,Hanscom Air Force Base~ MA 01731 (1979). 6. F. X. Kneizys, E. P. Shettle, W. O. Gallery, J. H. Chetwynd Jr., L. W. Abreu, J. E. A. Selby, R. W. Fenn and R. A. McClatchey, Atmospheric Transmittance/Radiance--Computer code LOWTRAN 5. Report AFGL-TR-80-0067, AFGU, Hanscom Air Force Base, City (1980). 7. C. R. Prasad and Prabha Venkatesh, Final report of ISRO RESPOND project, Model for atmospheric transmittance and radiance, Dept. Mech. Engg., Indian Institute of Science, Bangalore, India (1986). 8. S. Chandrasekhar, Radiative Tran~fer. Dover, New York (1950). 9. K. L. Coulson, J. V. Dave and Z. Sekera, Tables Re-

lated to Radiation Emerging from a Planetmy Atmosphere with Rayleigh Scattering. University of

California Press, Berkeley and Los Angeles (1960). 10. M. Kerker, The Scattering o f Light. Academic Press, New York (1969). 11. J. V. Dave, Subroutines for computing the parameters of electromagnetic radiation scattered by a sphere. IBM Scientific Center, Palo Alto, California, Report No. 320 (May 1968). 12. L. G. Hyeney and J. L. Greenstein, Diffuse radiation in the galaxy. Astrophys. J. 93, 73 (1941). 13. J. A. Coakley, and P. Chylek, The two-stream approximation in radiative transfer--includingthe angle of the incident radiation. J. Arm. Sci. 32, 409 (1975). 14. E. P. Shettle and J. Wienmann, The transfer of solar irradiance through inhomogeneous turbid atmospheres evaluated by Eddington's approximation. J. Atm. Sci. 27, 1048 (1969). 15. K. N. Liou, Analytic two-stream and four-stream so-

532

16. 17. 18. 19. 20. 21. 22.

23.

24.

25.

C . R . PRASADet al. lutions for radiative transfer. J. Atm. Sci. 31, 1473 (1974). W. M. Irvine, Multiple scattering by large particles: II--Optically thick layers. Astrophys. J. 152, 823 (1966). H. C. van de Hulst, A new look at multiple scattering, Tech. Rep. Goddard Instt. for Space studies, NASA, New York (1963). J. E. Hansen, Radiation transfer by doubling very thin layers. Astrophys. J. 155, 565 (1969). A. A. Lacis and J. E. Hansen, A parametrization for the absorption of solar radiation in the earth's atmosphere. J. Atm. Sci. 31, 113 (1974). B. M. Herman and S. R. Browning, A numerical solution to the equation of radiative transfer. J. Atm. Sci. 22, 559 (1965). B. M. Herman and S. R. Browning, A semi-analytic technique to integrate the radiative transfer equation over optical depth. J. Atm. Sci. 37, 1828 (1980). J. V. Dave and J, Gazdag, A modified Fourier transform method for multiple scattering calculations in a plane parallel Mie tatmosphere. Appl. Opt. 9, 1457 (1970). N. Brasslau and J. V. Dave, Effects of aerosols on the transfer of solar energy through realistic model atmospheres, Parts I and II. J. Appl. Meteorol. 30, 601 (1973). K. N. Liou, A numerical experiment on Chandrasekhar's discrete-ordinate method for radiative transfer: Applications to cloudy and hazy atmospheres. J. Atm. Sci. 30, 1303 (1973). J. V. Dave, A direct solution of the spherical har-

26. 27. 28. 29.

30. 31.

32. 33.

monics approximation to the radiative transfer equation for an arbitrary solar elevation, Part I: Theory. J. Atm. Sci. 32, 790, Part II: Results. J. Atm. Sci. 32, 1463 (1975). G. N. Plass and G. W. Kattawar, Calculations of reflected and transmitted radiance for earth's atmosphere. Appl. Opt. 7, 1129 (1968). K. Ranholm, M. B. Baker, and H. Harrison, Radiative transfer through media with uncertain or variable parameters. J. Atm. Sci. 37, 1279 (1980). R. E. Bird, R. L. Hulstrom and L. J. Lewis, Terrestrial solar spectral data sets. Solar Energy 30, 563 (1983). F. X. Kneizys, E. P. Shettle, W. O. Gallery, J. H. Chetwynd Jr., L. W. Abreu, J. E. A. Selby, S. A. Clough and R. W. Fenn, Atmospheric transmittance/ radiance: Computer code LOWTRAN 6. Report AFGL-TR-83-0187, Hanscom Air Force Base, MA 01731 (1983). A. K. Inamdar, Solar spectral diffuse radiance for a specified atmosphere. M. Sc. thesis, Dept. Mech. Engg., Indian Institute of Science, Bangalore (1984). M. Abramowitz and I. A. Stegun, Handbook o f Mathematical Functions, ed. by M. Abramowitz and I . A . Stegun, National Bureau of Standards, Applied Mathematics Series, 55 (1964). J. V. Dave, Extensive data sets of the diffuse radiation inrealistic atmospheric models with aerosls and common absorbing gases. Solar Energy 21, 361 (1978). J. E. Hansen, Exact and aproximate solutions for multiple scattering by cloudy and hazy planetary atmospheres. J. Atm. Sci. 26, 478 (1969).