- Email: [email protected]
A note on the variational method of Stokes and DeMarcus for radiative transfer in planetary atmospheres
ICARUS S9, 177--187 (1984)
A Note on the Variational Method of Stokes and DeMarcus for Radiative Transfer in Planetary Atmospheres 1 K E N N E T H S. K. CHOW, A. JAMES FRIEDSON, AND YUK L. YUNG 2 Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, California 91125 Received June 15, 1983; revised April 16, 1984 A generalized functional which yields the Milne integral equation on variation and whose extremum value is proportional to the reflectivity at arbitrary emergent angle is proposed. A similar functional exists for computing the transmissivity at arbitrary emergent angle. This work is a generalization of the variational method of Stokes and DeMarcus (1971, Icarus 14, 307) based on the principle of reciprocity. In the special case of trial functions that are linear in the undetermined parameters, the calculation is greatly simplified. The computational value of our variational principle is demonstrated.
INTRODUCTION
A number of authors (Huang, 1952a,b; Stokes and DeMarcus, 1971; Sze, 1976) have employed the variational method for solving the planetary problem in radiative transfer. Huang worked on conservative scattering in a homogeneous isotropic atmosphere of infinite depth. He constructed a functional which, upon variation, yields the integral equation first derived by Milne. A linear combination of five simple functions was used as the trial source function; and an accuracy of around a few parts per million was achieved for the specific intensity. Stokes and DeMarcus (1971) first pointed out that the reflectivity in the backward direction is proportional to the extremum value of a functional. The reflectivity is then accurate to second order. For isotropic scattering in a homogeneous atmosphere, they used a second-order polynomial trial source function and achieved an accuracy of a few parts in ten thousand for the reflectivity in the backward direction. Sze (1976) used a hybrid of variational and i Contribution 3894 of the Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, Calif. 91125. 2 To whom all correspondence should be addressed.
iterative methods to calculate the reflectivity of a finite atmosphere at arbitrary emergent angles. Without the full benefit of a functional, his trial source functions have to be much more sophisticated in order to match Stokes and DeMarcus' accuracy. In this paper we propose a generalization of the variational principle of Stokes and DeMarcus (1971) by the construction of two bifunctionals. The equation of radiative transfer is obtained by imposing the stationary condition on such functionals and the reflectivity and transmissivity at arbitrary angle are given by their extremum values. Further simplification is possible when one uses trial functions that are linear in the variational parameters, in which case the source function derived from the method of Stokes and DeMarcus can give the reflectivity at arbitrary emergent angle to second order. Finally, we demonstrate the usefulness and explore the limitations of our variational principles with extensive numerical results. THEORY
Following the notation of Chandrasekhar (1960), the equation of transfer in a planeparallel atmosphere with an isotropic scatterer is 177 0019-1035/84 $3.00 Copyright© 1984by AcademicPress. Inc. All rightsof reproductionin any formreserved.
CHOW, FRIEDSON,
178 /J $ (r,P) = ‘(r,P)
-
AND YUNG
emergent zenith angle and the second is the WO(')J('),
(I)
solar zenith angle).
Now let (y(t) and P(t) be two arbitrary integrable functions. Consider the functional,
where J(r)
= ; 1, Qr,~u) dp I
I(T,-l-d
W(a$)
1
OS/L5
normal optical depth where
cosine of zenith angle
c’=
wo(7) = single scattering
albedo.
S(s,wc)=
J(T) is the mean intensity.
It is directly proportional to the source function for an isotropically scattering atmosphere. J(T) satisfies the Milne integral equation F J(T) = $
r T'~(b
wo(s)J(s)E,l~ - s/ ds
(2)
G(t,s) = ;
Z/wo(s)wo(t) E,lr - sl. (5)
This functional is a generalization of the one written by Stokes and DeMarcus (1971). Imposing the condition that W be stationary with iespect to variation in a(t) and /3(t) yields two equations,
where PO = cosine of the solar zenith angle TF = incident solar flux normal to the atmosphere (note that TF in Chandrasekhar (1960) is defined to be normal to the direction of propagation) En(x) =
G(T,s)(Y(s) ds (6a)
(y(T)= s(T,Po) + and P(T)
=
S(T,~)
+
1;’G(T,s)@s)
ds, (6b)
where the functions S and G are as given in Eq. (5). CYand /I can be identified as proportional to the mean intensity for sunlight incident at angles ~0 and CL,respectively. In
The solution of Eq. (I) is related to J(T) fact upon putting >
(Y(T)= V/WO(T)J(T)
(3a) p(T) = aK(T) and substituting the expressions for (Y,p, S, and G into (6a) and (6b) we get +
00(s) J(s)r-l~-.‘“w
$.
(3b)
J(T) = &
In particular we are interested in calculating the reflCCtiVity R(@,j_q~) = I+(T = 0,p,~~)/F and the transmissivity T(p,po) = 1-(~,,p,p~)/F (where the first angle is the
P "w
+;
and
" wo(s)J(s)E,I~
- s/ ds
(6~)
RADIATIVE TRANSFER IN ATMOSPHERES
179
'/
I'(r =0,~)
F
r(~') = ~
e -~/"
I
+
1
1*(r : O,,u. o)
1
f~' too(s)K(s)E, Ir - sl as.
(6d)
J(~') is then readily identified as the mean intensity when sunlight is incident on the atmosphere at an angle/z0 and K(r) is the corresponding mean intensity when the incident angle is tz. 7rF is the flux normal to the atmosphere. See Fig. 1 for details. The extremum value of W is obtained by substituting (6c) into (4); we have Wex =
~,~F
/~
I~F/F.
~°/ F / u °
r
:o
~o(r)
T=
TI
e-S/t'oMs)J(s) ds
F2 = ~- R(~,~o).
(7)
The principle of reciprocity (Minnaert, 1941; van de Hulst, 1980) requires R(~,/.q~) = R(/~o,/~) 4 F = ~ fo' e-S/~otoo(s)K(s) ds 4/x0
l'(rr~o)
l-(ri,M. )
FIG. 1. Schematic diagram showing the geometry of the radiative transfer problem. The incident solar fluxes given are normal to the direction of propagation (see text).
(8)
provided that the incident normal flux component is the same. Reciprocity has been incorporated into the functional by virtue of its symmetry in tt and it0. Equation (8) can be derived from the functional by substituting (6d) into (4) and comparing it with (7). On setting tz = /x0 in W(a,fl), our results reduce to those of Stokes and DeMarcus (1971). In principle, the mean intensity derived from their functional can also be used to calculate the reflectivity at any zenith angle using Eq. (3a). The advantage of the functional W(a,t~) is that its extremum value gives the reflectivity to second order in the error of the trial a and ft. We surmise that as a direct consequence of the principle of reciprocity, similar variational methods exist for transfer problems with arbitrary phase functions. Unfortunately, we are at present unable to ascertain the validity of our conjecture in this work. A similar functional can be written which
yields the transmissivity T(g,~o) on variation,
W'(J,K) = fo~ OJo(s)K(s) F
4tt---ooe-'/"° ds + _ff__Fe_,q_~)/t, ds + f~t
4tx
1
too(s)J(s) J(t)K(s)
tOo(S)wo(t)Ellt - s I ds dt - fo' tOo(S)J(s)K(s) ds.
Since we are solving the same problem, this functional must give the same J(t) and K(t). A little reflection reveals that the same reasoning leading to (7) and (8) also applies here, giving us T(/x,/~o) = T(/xo,/X) 4 = ~
W'x.
180
CHOW, FRIEDSON, AND YUNG ANALYSIS
OF ACCURACY
Analytic expressions for J0") and KiT) can be obtained by choosing some trial functions Jt('r) = Jt(ai, a2 . . . . .
a,z, "r)
Kt(bl, b2 . . . . .
Kt(7") =
J
on W, Oai
OW -
--
Obi
-
0.
ai~J(T)
Jt("r) = Z
b,+, r)
with {ai} = al . . . . . a,, and {bi} - b, . . . . . b,, as unknown p a r a m e t e r s to be determined by imposing the e x t r e m u m condition
OW
above. The validity of these predictions remains to be established by computations. In actual calculations, it is convenient to choose trial functions that are linear in the undetermined parameters:
(9)
With Jt(~') and Kt0-) determined, R(IX,IX0) can be calculated by the three ways:
Kt("r) = ~', bj~qj(r)
where {'1'~i} is s o m e convenient set of functions a n d j runs f r o m 1 to a small integer. In this case there exists further simplification. The variational method for obtaining the {ai} and {hi} then gives F 2
~-- R3(/+,/-+o) = W ( ~ o J t , ds ~ Kt(s)Ogo(s)e .'/,~o- IX~)
I f T,
Rt(IX,IXo) = ~
(10a)
(12)
i
~ff-~oK,)
F2
OW
= ~ R2(~,ixo) + ~ b~ (9hi
(~'~ Jt(s)wo(s')e '/~' ds R2(IX,IXo) = ~1 L /~
(lOb)
F2 = ~ - R2(tx,/Xo)
4 W(X/-~oJ,, ~ o K , ) . R)(~,IX0) = ~..
([0c)
F2 OW = --4-Rffix,lXo) + ~ aj Oaj
Intuitively, one expects R3 to be the most accurate since R~ and R2 are accurate to first order, while R3 is accurate to second order as shown in 1
( " 8K(s)wo(s)e '/~ -ds 8R, = F Jo
ilia)
Ixo
l
682 = F [~' Jo 8J(s)ooo(s)e .,/u -ds tx
6R3-
(13)
where (9) has been used. R~, R2, and R3 are therefore all equal. Further,
OW OW Oai (ai'bD = ~ (bD (lib)
oW OW Obi (aj,bk) = ~ / ( a ~ )
4 F2 f~ 6J(s)6K(s)wo(S) ds 7-~ r I
X/-~o(s)6K(s) ds dt
F2 = ~ - Rff/z,/Xo)
(1 lc)
where ~ J ( s ) = Jt(s) -- Je(s) Jt = trial J J~ = exact J.
N o t e that 8R3 is indeed of second order because f0'r I f0T I ]G(t,s)[ ds dt is bounded from
(14)
i.e., the set of n equations OW/Oai = 0 do not involve the ai's and the other n(OW/Obi = 0) do not involve the bi's. This simplification arises f r o m the linearity of W in J and K and the linearity of J (and K) in the ai's (and the b~'s). If one has chosen some nonlinear trial functional form for J or K, such as the Pade a p p r o x i m a t i o n (Gragg, 1972), then none of the a b o v e simplification will occur. The decoupling of the two sets of equations enables us to solve only one of t h e m for {ai} (or {bi}) in order to calculate R2 (or R,). The
RADIATIVE TRANSFER IN ATMOSPHERES R(/.~,/z0) obtained is then accurate to second order. Let us define T1, T2, and T3 in a manner similar to R1, R2, and R3. The derivation of (14) can be carried through for T exactly as before and we have, W'(J,,Kt)
F~ = -~- T3(/x,/Zo) F2 = - T T2(z,tZo) F 2
= - 4 T,(~,tz0).
In particular, Tz(/Z,/z0) can be calculated as in (10b) by substituting e -s/~' with e -(~ -.,.v~, There is a fundamental limitation to the accuracy of reflectivity and transmissivity evaluated by variational methods. At grazing angle, it can be shown (Chandrasekhar, 1960) that FR(#
where J(0) and J('/'l) are the mean intensities evaluated at the upper and lower boundary, respectively. H e r e the accuracy of reflectivity and transmissivity is obviously of first order. The reason that R and T can be computed to higher accuracy than the source function is that they can be expressed as integrals of the source function. As shown by Sze (1976), the enhanced accuracy for R is due to the cancellation of errors in the integration. At grazing angle, there is no cancellation as the integral reduces to the function at the boundary. Thus we expect the accuracy of R and T t o divide into approximately two regions; a first-order region for emergent angles 0 -
= O,l~o) = 1+(7. = O , t z = O,txo)
RESULTS AND DISCUSSION
= oJ0(0) J ( 0 )
and FT(ix
181
= O,/zo) = 1-(7" = "rl,tZ = O,Izo) = O,~0(7"1)J(,/-1)
A large n u m b e r of simple cases have been studied to illustrate the ideas and formulae developed above. Table I is a list of the cases and their methods of derivation. TABLE I
LISTING OF THE PHYSICAL PARAMETERS OF CASES STUDIED, AND THE METHODS USED" Table
Case
~-~
COo
/Zo
Method
2 2 2 2 2 2 2 3 3 4
A B C D E F G H I J
1 1 1 1 1 1 4 4 4 ~
1 0.9 0.6 0.2 1 0.9 1 0.9 1 1
1 1 1 1 0.1 0.1 1 0.5 0.3 !
5
K
1
1
I
5
L
1
1
0.1
Variational three-term polynomial trial function Variational three-term polynomial trial function Variational three-term polynomial trial function Variational three-term polynomial trial function Variational three-term polynomial trial function Variational three-term polynomial trial function Variational six-term polynomial trial function Variational, six-term polynomial trial function Variational, six-term polynomial trial function Variational, five-term exponential and exponential integral functions as trial function Variational, with three-step functions as trial functions, interactive i m p r o v e m e n t , result fitted to five-term polynomial Variational, with three-step functions as trial functions, interactive i m p r o v e m e n t , result fitted to five-term polynomial
" T h e results are s u m m a r i z e d in Tables I1-V.
182
CHOW,
FRIEDSON,
AND YUNG
TABLE 11 REFLECTION (R) A N D T R A N S M I S S I O N ( i t ) F U N C T I O N S FOR A H O M O G E N E O U S A I M O S P H E R E W I T H 71 = I A N D V A R I O U S V A L U E S OF (oil A N D /~ll (CASES A - F )
/,
0
0.1
0.3
0,5
11.7
R Ro A T T, A
0.463 0.439 0.023 0.243 0.241 0.001
17 34 83 37 51 86
0.467 0.464 0.003 0.292 0.293 0.001
91 88 03 06 40 34
Case A r t I, ¢o~ 0.436 71 0.380 42 0.436 66 0.380 51 0.000 1)5 0.000 119 0.336 24 0.3211 93 0.336 36 0.320 91 -0.000 12 0.000 02
I, /x~ - I 0,328 82 1/.328 91 -0.000 09 0.290 32 0.290 30 0.000 02
R Ro A T 7~ A
0.377 0.359 0.017 0.187 0.186 0.001
24 37 87 29 14 15
0.374 0.372 0.002 0.223 0.224 0.001
99 67 32 01 03 02
Case B ~-i I, 0.344 31 0.298 0.344 24 0.298 0.000 07 0.000 0.257 48 0.246 0.257 58 0.246 -0.0011 IO 0.000
0.9, #~ 0.256 0.256 0.000 0.223 0.223 0.000
R Ro A T 7]~ A
0.200 0.194 0.006 0.087 0.086 0.000
54 11 43 39 91 48
0.191 0.191 0.000 0.101 0.101 -0.000
96 12 84 43 76 33
Case C r~ I. ¢o~ 0.169 38 0.144 42 0.169 33 0.14443 0.000 (15 -0.000 01 0.117 66 0.113 80 0.117 70 0.113 80 0.000 114 0.000 00
R Ro A T 7], A
0.054 0.053 0.000 0.021 0.020 0.000
10 70 40 08 88 20
0.050 0.050 0.000 0.023 0.023 0.000
03 00 03 61 59 02
Case D TT I, 0.042 48 I).035 0.042 47 0.035 0.000 01 0.000 0.027 38 0.026 0.027 38 0.026 0.000 00 0.000
R Ro A l" l~ A
2.132 2.986 -0.853 0.435 0.138 0.296
60 22 62 83 98 85
1.629 1.779 -0.150 0.272 0.177 0.094
58 64 06 07 54 53
Case E 7i I, (o~ 1.054 87 0.776 19 1.058 44 0.771 82 0.003 57 11.004 37 0.279 32 11.316 92 0.273 71 0.318 63 0.005 61 -0.001 71
R Ro A T To A
1.826 2.614 --0.787 0.360 0.091 0.268
56 00 44 17 56 63
1.377 1.516 --0. 139 0.204 0.116 0.087
19 58 39 00 90 10
¢o~ l0 15 05 68 68 00
0.9
F-
I).286 0.287 -0.000 0.260 0.260 0.000
96 03 07 20 18 02
0.269 32 1/.269 39 -0.000 07 0.246 54 0.246 53 0.0110 01
0.341 0.341 0.000 0.658 0.658 0.000
37 33 04 63 67 04
0.223 0.223 0.000 0.200 0.200 0.000
84 88 04 71 70 01
0.209 0.210 -0.000 0.190 0.190 0.000
96 00 04 27 26 01
0.267 0.267 0.000 0.591 0.591 0.000
46 41 05 59 63 04
0.100 41 0.10042 -0.000 01 0.088 68 0.088 68 0.000 00
(l. 129 0.129 1t.000 0.471 0.471 0.000
54 51 03 36 37 01
I 92 97 05 65 63 02
0.6,/*o I 0.123 57 0.12358 I).000 01 0.103 75 0.103 75 0.000 00
(oo = 0.2, /*o 67 0.030 67 0.030 00 0.000 73 0.024 73 0.024 00 0.000
1.0
I).107 21 0.10722 -0.000 01 0.093 43 1/.093 43 0.000 O0
I 29 29 00 51 51 00
I,/*o - (J. I 0.614 20 11.610 02 0.004 18 0.316 99 0.318 98 0.001 99
Case F rl I, oJo - 0.9, /*o = 0.1 0.87t) 65 0.631 78 0.495 83 0.875 57 0.628 83 0.492 72 --0.004 92 0.002 95 0.1103 11 0.200 39 0.234 86 0.238 97 0.194 25 0.235 79 0.240 40 0.006 14 0.000 93 0.001 43
0.026 1/.026 0.000 0.022 0.022 0.000
17 17 00 15 15 00
0.024 0.024 0.000 0.021 0.021 0.000
47 47 00 05 05 00
0.031 0.031 0.000 0.392 0.392 0.000
98 98 00 26 26 00
0.508 0.504 0.003 0.301 0.303 --0.001
25 86 39 68 24 56
0.467 0.464 I).003 0.292 0.293 -0.001
91 88 03 06 40 34
0.695 0.697 --0.002 0.304 0.302 0.002
36 66 30 64 34 30
1/.408 13 lL405 55 0.002 58 0.229 61 11.230 78 -0.001 17
0.374 0.372 0.002 0.223 0.224 0.001
98 67 31 01 03 02
0.566 0.568 -0.002 0.228 0.225 0.002
011 99 99 34 75 59
N o t e . The variation principle with three-term polynomial is used. The " e x a c t " values Ro and To are taken from van de Hulst 119801. The symbol A denotes the difference between "'exact" and approximate solutions.
RADIATIVE TRANSFER IN ATMOSPHERES All integrals are p e r f o r m e d numerically by Gaussian quadrature. The column F ÷,- in the tables serves as a further c h e c k on our calculations. F + (entered in row R) is defined as
Ez(s)J(s) ds,
F + = 2too
which is the u p w a r d flux. F - (entered in row T) is defined as F - = 2too fo' E2(71 - s)J(s) ds + Fe-~l/~O, which is a n o t h e r expression for the d o w n w a r d flux. Their " e x a c t " values are taken f r o m van de Hulst (1980). Table II shows the results of calculations using a s e c o n d - o r d e r polynomial trial function for an a t m o s p h e r e of optical thickness rl = I. too was chosen to be constant for
183
convenience. Typically the error stays small (a few parts per ten thousand) until/z goes below 0.3. In case E and F, the solar zenith angle /z0 is small and the error is significant as expected. Table III contains results for an atmosphere of larger optical depth, r~ = 4. A fifthorder polynomial trial function was used to achieve c o m p a r a b l e a c c u r a c y as in Table II. Again, the error increases f o r / z -< 0.3. F o r an infinite a t m o s p h e r e , we took the results of H u a n g (1952b) and reconstructed his calculations to eight significant figures in Table IV. We used double precision on an I B M 370/3032 and achieved higher numerical a c c u r a c y than did Huang. H u a n g ' s choice of trial function gave him an accuracy of a few parts in ten million. N o t e also that A grows as /z goes to zero. The fluctuating sign of A for /z --< 0.3 has little significance b e c a u s e the last significant fi-
T A B L E 111 REFLECTION ( R ) AND TRANSMISSION ( T ) FUNCTIONS FOR CASES G - I W I T H 7"I = 4. USING SIx-TERM POLYNOMIAL
0
0.1
0.3
0.5
0.7
0.9
1.0
F ~,-
R R0 A T To A
0.623 0.593 0.030 0.142 0.132 0.010
61 10 51 76 28 48
0.661 0.657 0.004 0.167 0.164 0.002
81 37 44 11 76 35
CaseGrt =4, 0.698 63 0.705 0.698 62 0.705 0.000 01 --0.000 0.216 l0 0.262 0.216 05 0.262 0.000 05 -0.000
to0 = l,/z0 = 1 77 0.697 05 98 0.697 23 21 --0.000 18 18 0.302 22 30 0.302 34 12 -0.000 12
0.679 0.679 --0.000 0.333 0.333 -0.000
31 46 15 42 52 10
0.668 0.668 --0.000 0.345 0.345 -0.000
41 55 14 48 58 10
0.690 0.690 --0.000 0.309 0.309 0.000
92 93 01 07 07 00
R R0 A T To A
0.732 0.698 0.043 0.032 0.027 0.005
97 77 20 44 12 32
0.686 0.682 0.004 0.035 0.033 0.001
99 18 81 27 62 65
Case H ~'1 = 4, w0 = 1,/x0 = 0.5 0.606 59 0.541 56 0.488 84 0.606 52 0.541 72 0.488 99 0.000 07 --0.000 16 -0.000 15 0.045 28 0.058 43 0.073 66 0.045 18 0.058 49 0.073 72 0.000 10 --0.000 06 -0.000 06
0.445 0.445 -0.000 0.088 0.088 --0.000
34 45 11 66 72 06
0.426 0.426 -0.000 0.095 0.095 -0.000
31 41 10 47 52 05
0.504 0.504 0.000 0.071 0.071 0.000
14 09 05 84 82 02
R R0 A T To A
1.327 1.292 0.034 0.071 0.075 -0.004
14 99 15 63 88 25
1.186 1.185 0.000 0.092 0.094 -0.001
15 98 17 75 68 93
0.727 0.727 0.000 0.205 0.205 0.000
60 59 01 67 60 07
0.698 0.698 0.000 0.216 0.216 0.000
63 62 01 10 05 05
0.825 0.825 -0.000 0.174 0.174 -0.000
01 02 01 97 98 01
Note.
Case I r1 = 4, too = 1,/x0 = 0.999 49 0.880 57 0.794 0.999 64 0.880 57 0.794 - 0 . 0 0 0 15 0.000 00 0.000 0.124 73 0.153 55 0.181 0.124 80 0.153 46 0.181 - 0 . 0 0 0 07 0.000 09 0.000
Values of R0 and To are taken from van de Hulst (1980).
0.3 85 84 01 22 14 08
184
CHOW, FRIEDSON, A N D YUNG T A B L E IV
REFLECT|VITY (R) AND TRANSMISSION (T) FOR CONSERVATIVE I N F I N I T E ATMOSPHERE (CASE J)
#
R
0 0.1 0.2 0.3 0.4 0.5 0.6 0,7 0.8 0.9 1.0
0.726 0.824 0.878 0.918 0.949 0.975 0.996 1.015 1.030 1.044 1,056
Ro
9193 3332 6156 4890 8546 4632 8943 1573 9418 7410 9203
0.726 0.824 0.878 0.918 0.949 0.975 0.996 1.015 1.030 1.044 1.056
A
9527 3314 6140 4891 8548 4633 8942 1573 9418 7410 9203
0.000 0.000 0.000 0.000 0.000 0,000 0.000 0.000 0.000 0.000 0,000
(1334 0018 0016 0001 0002 0001 0001 0000 0000 0000 0000
Note. T h e v a l u e s for R are c o m p u t e d by us using the p a r a m e t e r s o b t a i n e d by H u a n g (1952b); t h o s e for R0 are c o m p u t e d u s i n g the results o f P l a c z e k (19471.
-i
IOt
,
I
'
O-IL
I
,
I
I
cose A tOo: I
r I:
I
I
gure is probably s w a m p e d by round-off error. Table V s h o w s the results o f an alternative use o f the functional. Instead o f deriving the source function directly from the functional, S z e ' s m e t h o d (Sze, 19761 o f a N - s t e p approximation ( N = 3 in our case) and then iteration w a s used. The iterated J and K were finally fitted to a fifth-order polynomial for numerical integration. Since the J(t) and K(t) used are not directly derived from the function, R~, T1 and R2, T2 are no longer equal to R3, /2~ (to s e c o n d order). The data demonstrated that the error o f the functional m e t h o d is significantly smaller than that o f direct convolution. The variation o f error with /z for a polynomial trial function is s h o w n graphically in Fig. 2 (for ~'j = 11 and Fig. 3 (for 7l = 4). The dip in the curve at around /z - 0.3 is not a true m i n i m u m , but a zero crossing w h e n 2~ changes sign as p, is decreased. Figures 4 and 5 s h o w the i m p r o v e m e n t o f the functional o v e r direct convolution. We
I p.O= I '~-°l
'
I
1(52"-~ 3-
~
\\,
'
I
I
I
'
cose G OJo=l T,=4
I
'
/~o=1
IR-%i
Izxl \
5
\
x
x
_x
~x~
3
"* ~ f IT_ToI
io-5
IAI
\~\\
IR -%1
~ ,
./'-~
I T # ; [- . . . . . .
10-6 3 10-7
0
J
I
,
0.2
FIG. 2. Variation
I
,
0.4
[
,
0.6
I 0.8
5
i L0
/z of
18RI and 18TI with/z
for case A of
Table I. Both ISRIand 18Tl g o t h r o u g h a c h a n g e in sign at a r o u n d p, = 0.3 a n d t h e r e f o r e the dip in the graph. N u m e r i c a l v a l u e s o f R0 a n d To are t a k e n f r o m v a n de H u l s t (1980) for o d d v a l u e s o f p. a n d f r o m Carlstedt and Mullikin (1966) for e v e n v a l u e s o f p,.
id s 5
id7 0
J
I O.2
,
I
J
0.4
1 0.6
~
I 0.8
/z FIG. 3. S a m e as Fig. 2 e x c e p t 7, = 4.0.
i I.O
RADIATIVE
TRANSFER
IN ATMOSPHERES
185
TABLE V THREE DIFFERENT EXPRESSIONS FOR EVALUATING THE REELECTION AND TRANSMISSION FUNCTIONS AS DEFINED aV EQ. (10) ARE COMPARED 0.1
0.3
0.5
0.483 0.018 0.462 -0.001 0.465 0.000
28 40 95 93 21 33
0.437 0.001 0.433 -0.003 0.436 0.000
84 18 58 08 66 00
TI AI T2 A2 7"3 A3
0.316 0.023 0.296 0.002 0.293 0.000
88 48 34 94 90 50
0.340 0.004 0.336 --0.000 0.336 0.000
53 17 11 25 45 09
0.321 0.000 0.319 --0.001 0.320 0.000
Ri A1
1.728 --0.051 1.728 --0.051 1.769 -0.009
16 48 16 48 66 98
1.042 --0.015 1.067 0.008 1.056 -0.001
77 67 39 95 98 46
Case L T1 = l, 0.764 --0.007 0.790 0.018 0.771 -0.000
0.209 0.032 0.209 0.032 0.175 -0.002
81 27 81 27 50 04
0.286 0.013 0.309 0.036 0.273 0.000
76 05 72 01 85 14
A2 R3
A3 TI A1 /'2 A2
T3 A3
0.9
Case K r~ = 1, too = 1,/-I-o = 1 0.379 03 0.326 98 -0.001 48 -0.001 93 0.377 79 0.326 59 --0.002 72 --0.002 32 0.380 48 0.328 87 - 0 . 0 0 0 03 - 0 . 0 0 0 04
Rt A1 RE A2 R3 A3
R2
0.7
0.326 0.007 0.351 0.032 0.319 0.000
1.0
0.285 -0.001 0.285 -0.001 0.286 -0.000
12 91 04 99 99 04
0.267 --0.001 0.267 -0.001 0.269 --0.000
54 85 54 85 35 04
49 81 07 23 29 01
0.259 -0.001 0.258 0.001 0.260 -0.000
04 14 95 23 17 01
0.245 -0.001 0.245 --0.001 0.246 -0.000
32 21 32 21 51 02
~o = l, t~o = 0. l 31 0.605 88 51 --0.004 14 53 0.629 91 71 0.019 89 66 0.610 23 16 0.000 21
0.502 -0.002 0.523 0.019 0.505 0.000
40 46 90 04 18 32
0.462 -0.001 0.483 0.018 0.465 0.000
95 93 28 40 21 33
0.306 0.003 0.328 0.025 0.303 0.000
68 44 28 04 75 51
0.296 0.002 0.316 0.023 0.293 0.000
34 94 88 48 89 49
20 29 86 05 93 02
0.289 -0.000 0.289 --0.001 0.290 -0.000
19 56 48 85 13 50
0.323 0.004 0.347 0.028 0.319 0.000
90 92 65 67 52 54
Note. T h e trial functions are obtained in three steps. First the variation principle using three-step functions is used. T h e result is refined by one iteration (Sze, 1976) using the integral equation, and finally fitted to a fifth-order polynomial. started with a three-step approximation and the numerical integration was done using a fifth-order polynomial as described earlier. A stands for the difference as in Table II, its improvement ranges from a factor of 2 to two orders of magnitude. Another case o f interest is that of a Lambert surface (with albedo h) at the lower boundary of the atmosphere. In this case, the functions S ( s , / z ) a n d G(t,s) i n E q . (5) are modified
S (s ,~ ) =
to
~x/-L~o~ s)
[ F e-s/~' + --~ hF e-'l#'OEz('rl - s) ]
x/co0(t)o~0(s)
G(t,s) =
2
[Ellt Otherwise, through
s I + 2XE2(rj the
as before.
similar degree
entire
-
t)Ez(rl -
calculation
We obtained
of accuracy
s)]. goes
results with
as in the case h
=0.
CONCLUSION
We have provided a generalization of the method of Stokes and DeMarcus (1971) to calculate reflectivity and transmissivity at arbitrary emergent angles. For the special
186
CHOW, FRIEDSON, AND YUNG I0'
'
I
'
I
'
I
'
I
3
o 10-2
\
x
case K ~o =I
L =I /~o =I
\%
3
i 0-~
~'---o-""
+
la,-al
1
3
IAI 10- 4 +
+
+~+~+
3
3 Id 6
3
i0-7 0
1
0.2
,
I
0.4
,
I
0.6
J
0.8
i
1.0
FIG. 4. Comparison of accuracy of the three methods of evaluating the reflectivity using the threestep approximate J as in Table V. Note that the functional method can be more accurate by one to two orders of magnitude. case of trial functions linear in the variational p a r a m e t e r s there is a degeneracy between the Stokes and D e M a r c u s method and ours. Error analysis has d e m o n s t r a t e d the difference b e t w e e n the two independent uses of the functional. If one uses trial functions of the form given by (12) and derives the m e a n intensity J from the functional, then the reflectivity and transmissivity can be c o m p u t e d to second order by direct convolution. It is u n n e c e s s a r y to evaluate K or the full functional. I f one uses a J (and K) not directly derived from the functional then the functional gives better results than direct convolution. There is a natural limitation on the a c c u r a c y of such methods in the sense that errors b e c o m e significant for # <~ 0.2. This is to be expected since small # is equivalent to large effective optical depth and the simple trial functions used in the sample calculations are insufficient.
The construction of our functional was motivated by Levine and Schwinger's (1950) variational principles of electromagnetic w a v e diffraction. The theory of radiative transfer is closely related to the theory of neutron diffusion. Therefore it is not surprising to detect striking similarity of our functional with the one given by Davison (1957). The relation between the functional and simple physical quantities, first discovered by Stokes and DeMarcus (1971), has been extended by us. Although this p a p e r does not significantly advance the computational aspects of the planetary problem beyond Sze (1976), it has brought new i n s i g h t - - t h e connection between the principle of reciprocity and the variational m e t h o d - - t o the problem of radiative transfer. I n d e p e n d e n t work by C h e y n e y and Arking (1976) and Shia and Yung (personal communication, 1984) indicates that generalization to problems with nontrivial phase function and non-plane-parallel g e o m e t r y can be made.
id' '
3
1
'
elK
'~Oo::
T,'=I I/.,o='1
]
~52 3
"',
i0-~
"%
..... I#_~:#J.....
I
5
IAI
10-4 3
3 1~56 3 i0 7
0
I
]
0.2
I
I
0.4
L
I
0.6
~
I
0.8
FIG. 5. Same as Fig. 4, for transmissivity.
1,0
RADIATIVE TRANSFER ACKNOWLEDGMENTS One of us (Y.L.Y.) acknowledges helpful discussions with N. D. Sze and R. M. Goody. This research was supported by NASA Grant NSG 7376 to California Institute of Technology.
REFERENCES
CARLSTEDT, J. L., AND T. W. MULLIKIN (1966). Chandrasekhar's X- and Y-functions. Astrophys. J.
Suppl. Ser. 113(12), 449. CHANDRASEKHAR, S. (1960). Radiative Transfer. Dover, New York. CHEYNEY, H., III, AND A. ARKING (1976). A new formulation for anisotropic radiative transfer problems: I. Solution with a variational technique. Astrophys. J. 207, 808. DAVlSON, B. (1957). Neutron Transport Theory. Oxford Univ. Press, London/New York. GRAGG, W. B. (1972). The Pade Table and its relation to certain algorithms of numerical analysis. SLAM. Review, 14, 1. VAN DE HULST, H. C. (1980). Multiple Light Scattering, Vol. I. Academic Press, New York.
IN ATMOSPHERES HUANG, blems with a HUANG, blems
187
S. S. (1952a). The variational method for proof radiative transfer: I. Isotropic scattering constant net flux. Astrophys. J. 117, 211. S. S. (1952b). The variational method for proof radiative transfer: III. Reflection effect. Astrophys. J. 117, 221. LEVINE,H., AND J. SCHWINGER (1950). On the theory of electromagnetic wave diffraction by an aperture in an infinite plane conducting screen. Commun. Pure Appl. Math. 3, 355. MINNAERT, M. (1941). The reciprocity principle in lunar photometry. Astrophys. J. 93, 403. PLACZEK, G. (1947). The angular distribution of neutrons emerging from a plane surface. Phys. Rev. 72, 556. STOKES, R. A., AND W. C. DEMARCUS (1971). A simple technique for calculating line profiles of inhomogeneous planetary atmospheres. Icarus 14, 307. SZE, N. D. (1976). Variational methods in radiative transfer problems. J. Quant. Spectrosc. Radiat. Transfer 16, 763. YUNG, Y. L. (1976). A numerical method for calculating the mean intensity in an inhomogeneous Rayleigh scattering atmosphere. J. Quant. Spectrosc. Radiat. Transfer 16, 755.
A Note on the Variational Method of Stokes and DeMarcus for Radiative Transfer in Planetary Atmospheres 1 K E N N E T H S. K. CHOW, A. JAMES FRIEDSON, AND YUK L. YUNG 2 Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, California 91125 Received June 15, 1983; revised April 16, 1984 A generalized functional which yields the Milne integral equation on variation and whose extremum value is proportional to the reflectivity at arbitrary emergent angle is proposed. A similar functional exists for computing the transmissivity at arbitrary emergent angle. This work is a generalization of the variational method of Stokes and DeMarcus (1971, Icarus 14, 307) based on the principle of reciprocity. In the special case of trial functions that are linear in the undetermined parameters, the calculation is greatly simplified. The computational value of our variational principle is demonstrated.
INTRODUCTION
A number of authors (Huang, 1952a,b; Stokes and DeMarcus, 1971; Sze, 1976) have employed the variational method for solving the planetary problem in radiative transfer. Huang worked on conservative scattering in a homogeneous isotropic atmosphere of infinite depth. He constructed a functional which, upon variation, yields the integral equation first derived by Milne. A linear combination of five simple functions was used as the trial source function; and an accuracy of around a few parts per million was achieved for the specific intensity. Stokes and DeMarcus (1971) first pointed out that the reflectivity in the backward direction is proportional to the extremum value of a functional. The reflectivity is then accurate to second order. For isotropic scattering in a homogeneous atmosphere, they used a second-order polynomial trial source function and achieved an accuracy of a few parts in ten thousand for the reflectivity in the backward direction. Sze (1976) used a hybrid of variational and i Contribution 3894 of the Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, Calif. 91125. 2 To whom all correspondence should be addressed.
iterative methods to calculate the reflectivity of a finite atmosphere at arbitrary emergent angles. Without the full benefit of a functional, his trial source functions have to be much more sophisticated in order to match Stokes and DeMarcus' accuracy. In this paper we propose a generalization of the variational principle of Stokes and DeMarcus (1971) by the construction of two bifunctionals. The equation of radiative transfer is obtained by imposing the stationary condition on such functionals and the reflectivity and transmissivity at arbitrary angle are given by their extremum values. Further simplification is possible when one uses trial functions that are linear in the variational parameters, in which case the source function derived from the method of Stokes and DeMarcus can give the reflectivity at arbitrary emergent angle to second order. Finally, we demonstrate the usefulness and explore the limitations of our variational principles with extensive numerical results. THEORY
Following the notation of Chandrasekhar (1960), the equation of transfer in a planeparallel atmosphere with an isotropic scatterer is 177 0019-1035/84 $3.00 Copyright© 1984by AcademicPress. Inc. All rightsof reproductionin any formreserved.
CHOW, FRIEDSON,
178 /J $ (r,P) = ‘(r,P)
-
AND YUNG
emergent zenith angle and the second is the WO(')J('),
(I)
solar zenith angle).
Now let (y(t) and P(t) be two arbitrary integrable functions. Consider the functional,
where J(r)
= ; 1, Qr,~u) dp I
I(T,-l-d
W(a$)
1
OS/L5
normal optical depth where
cosine of zenith angle
c’=
wo(7) = single scattering
albedo.
S(s,wc)=
J(T) is the mean intensity.
It is directly proportional to the source function for an isotropically scattering atmosphere. J(T) satisfies the Milne integral equation F J(T) = $
r T'~(b
wo(s)J(s)E,l~ - s/ ds
(2)
G(t,s) = ;
Z/wo(s)wo(t) E,lr - sl. (5)
This functional is a generalization of the one written by Stokes and DeMarcus (1971). Imposing the condition that W be stationary with iespect to variation in a(t) and /3(t) yields two equations,
where PO = cosine of the solar zenith angle TF = incident solar flux normal to the atmosphere (note that TF in Chandrasekhar (1960) is defined to be normal to the direction of propagation) En(x) =
G(T,s)(Y(s) ds (6a)
(y(T)= s(T,Po) + and P(T)
=
S(T,~)
+
1;’G(T,s)@s)
ds, (6b)
where the functions S and G are as given in Eq. (5). CYand /I can be identified as proportional to the mean intensity for sunlight incident at angles ~0 and CL,respectively. In
The solution of Eq. (I) is related to J(T) fact upon putting >
(Y(T)= V/WO(T)J(T)
(3a) p(T) = aK(T) and substituting the expressions for (Y,p, S, and G into (6a) and (6b) we get +
00(s) J(s)r-l~-.‘“w
$.
(3b)
J(T) = &
In particular we are interested in calculating the reflCCtiVity R(@,j_q~) = I+(T = 0,p,~~)/F and the transmissivity T(p,po) = 1-(~,,p,p~)/F (where the first angle is the
P "w
+;
and
" wo(s)J(s)E,I~
- s/ ds
(6~)
RADIATIVE TRANSFER IN ATMOSPHERES
179
'/
I'(r =0,~)
F
r(~') = ~
e -~/"
I
+
1
1*(r : O,,u. o)
1
f~' too(s)K(s)E, Ir - sl as.
(6d)
J(~') is then readily identified as the mean intensity when sunlight is incident on the atmosphere at an angle/z0 and K(r) is the corresponding mean intensity when the incident angle is tz. 7rF is the flux normal to the atmosphere. See Fig. 1 for details. The extremum value of W is obtained by substituting (6c) into (4); we have Wex =
~,~F
/~
I~F/F.
~°/ F / u °
r
:o
~o(r)
T=
TI
e-S/t'oMs)J(s) ds
F2 = ~- R(~,~o).
(7)
The principle of reciprocity (Minnaert, 1941; van de Hulst, 1980) requires R(~,/.q~) = R(/~o,/~) 4 F = ~ fo' e-S/~otoo(s)K(s) ds 4/x0
l'(rr~o)
l-(ri,M. )
FIG. 1. Schematic diagram showing the geometry of the radiative transfer problem. The incident solar fluxes given are normal to the direction of propagation (see text).
(8)
provided that the incident normal flux component is the same. Reciprocity has been incorporated into the functional by virtue of its symmetry in tt and it0. Equation (8) can be derived from the functional by substituting (6d) into (4) and comparing it with (7). On setting tz = /x0 in W(a,fl), our results reduce to those of Stokes and DeMarcus (1971). In principle, the mean intensity derived from their functional can also be used to calculate the reflectivity at any zenith angle using Eq. (3a). The advantage of the functional W(a,t~) is that its extremum value gives the reflectivity to second order in the error of the trial a and ft. We surmise that as a direct consequence of the principle of reciprocity, similar variational methods exist for transfer problems with arbitrary phase functions. Unfortunately, we are at present unable to ascertain the validity of our conjecture in this work. A similar functional can be written which
yields the transmissivity T(g,~o) on variation,
W'(J,K) = fo~ OJo(s)K(s) F
4tt---ooe-'/"° ds + _ff__Fe_,q_~)/t, ds + f~t
4tx
1
too(s)J(s) J(t)K(s)
tOo(S)wo(t)Ellt - s I ds dt - fo' tOo(S)J(s)K(s) ds.
Since we are solving the same problem, this functional must give the same J(t) and K(t). A little reflection reveals that the same reasoning leading to (7) and (8) also applies here, giving us T(/x,/~o) = T(/xo,/X) 4 = ~
W'x.
180
CHOW, FRIEDSON, AND YUNG ANALYSIS
OF ACCURACY
Analytic expressions for J0") and KiT) can be obtained by choosing some trial functions Jt('r) = Jt(ai, a2 . . . . .
a,z, "r)
Kt(bl, b2 . . . . .
Kt(7") =
J
on W, Oai
OW -
--
Obi
-
0.
ai~J(T)
Jt("r) = Z
b,+, r)
with {ai} = al . . . . . a,, and {bi} - b, . . . . . b,, as unknown p a r a m e t e r s to be determined by imposing the e x t r e m u m condition
OW
above. The validity of these predictions remains to be established by computations. In actual calculations, it is convenient to choose trial functions that are linear in the undetermined parameters:
(9)
With Jt(~') and Kt0-) determined, R(IX,IX0) can be calculated by the three ways:
Kt("r) = ~', bj~qj(r)
where {'1'~i} is s o m e convenient set of functions a n d j runs f r o m 1 to a small integer. In this case there exists further simplification. The variational method for obtaining the {ai} and {hi} then gives F 2
~-- R3(/+,/-+o) = W ( ~ o J t , ds ~ Kt(s)Ogo(s)e .'/,~o- IX~)
I f T,
Rt(IX,IXo) = ~
(10a)
(12)
i
~ff-~oK,)
F2
OW
= ~ R2(~,ixo) + ~ b~ (9hi
(~'~ Jt(s)wo(s')e '/~' ds R2(IX,IXo) = ~1 L /~
(lOb)
F2 = ~ - R2(tx,/Xo)
4 W(X/-~oJ,, ~ o K , ) . R)(~,IX0) = ~..
([0c)
F2 OW = --4-Rffix,lXo) + ~ aj Oaj
Intuitively, one expects R3 to be the most accurate since R~ and R2 are accurate to first order, while R3 is accurate to second order as shown in 1
( " 8K(s)wo(s)e '/~ -ds 8R, = F Jo
ilia)
Ixo
l
682 = F [~' Jo 8J(s)ooo(s)e .,/u -ds tx
6R3-
(13)
where (9) has been used. R~, R2, and R3 are therefore all equal. Further,
OW OW Oai (ai'bD = ~ (bD (lib)
oW OW Obi (aj,bk) = ~ / ( a ~ )
4 F2 f~ 6J(s)6K(s)wo(S) ds 7-~ r I
X/-~o(s)6K(s) ds dt
F2 = ~ - Rff/z,/Xo)
(1 lc)
where ~ J ( s ) = Jt(s) -- Je(s) Jt = trial J J~ = exact J.
N o t e that 8R3 is indeed of second order because f0'r I f0T I ]G(t,s)[ ds dt is bounded from
(14)
i.e., the set of n equations OW/Oai = 0 do not involve the ai's and the other n(OW/Obi = 0) do not involve the bi's. This simplification arises f r o m the linearity of W in J and K and the linearity of J (and K) in the ai's (and the b~'s). If one has chosen some nonlinear trial functional form for J or K, such as the Pade a p p r o x i m a t i o n (Gragg, 1972), then none of the a b o v e simplification will occur. The decoupling of the two sets of equations enables us to solve only one of t h e m for {ai} (or {bi}) in order to calculate R2 (or R,). The
RADIATIVE TRANSFER IN ATMOSPHERES R(/.~,/z0) obtained is then accurate to second order. Let us define T1, T2, and T3 in a manner similar to R1, R2, and R3. The derivation of (14) can be carried through for T exactly as before and we have, W'(J,,Kt)
F~ = -~- T3(/x,/Zo) F2 = - T T2(z,tZo) F 2
= - 4 T,(~,tz0).
In particular, Tz(/Z,/z0) can be calculated as in (10b) by substituting e -s/~' with e -(~ -.,.v~, There is a fundamental limitation to the accuracy of reflectivity and transmissivity evaluated by variational methods. At grazing angle, it can be shown (Chandrasekhar, 1960) that FR(#
where J(0) and J('/'l) are the mean intensities evaluated at the upper and lower boundary, respectively. H e r e the accuracy of reflectivity and transmissivity is obviously of first order. The reason that R and T can be computed to higher accuracy than the source function is that they can be expressed as integrals of the source function. As shown by Sze (1976), the enhanced accuracy for R is due to the cancellation of errors in the integration. At grazing angle, there is no cancellation as the integral reduces to the function at the boundary. Thus we expect the accuracy of R and T t o divide into approximately two regions; a first-order region for emergent angles 0 -
= O,l~o) = 1+(7. = O , t z = O,txo)
RESULTS AND DISCUSSION
= oJ0(0) J ( 0 )
and FT(ix
181
= O,/zo) = 1-(7" = "rl,tZ = O,Izo) = O,~0(7"1)J(,/-1)
A large n u m b e r of simple cases have been studied to illustrate the ideas and formulae developed above. Table I is a list of the cases and their methods of derivation. TABLE I
LISTING OF THE PHYSICAL PARAMETERS OF CASES STUDIED, AND THE METHODS USED" Table
Case
~-~
COo
/Zo
Method
2 2 2 2 2 2 2 3 3 4
A B C D E F G H I J
1 1 1 1 1 1 4 4 4 ~
1 0.9 0.6 0.2 1 0.9 1 0.9 1 1
1 1 1 1 0.1 0.1 1 0.5 0.3 !
5
K
1
1
I
5
L
1
1
0.1
Variational three-term polynomial trial function Variational three-term polynomial trial function Variational three-term polynomial trial function Variational three-term polynomial trial function Variational three-term polynomial trial function Variational three-term polynomial trial function Variational six-term polynomial trial function Variational, six-term polynomial trial function Variational, six-term polynomial trial function Variational, five-term exponential and exponential integral functions as trial function Variational, with three-step functions as trial functions, interactive i m p r o v e m e n t , result fitted to five-term polynomial Variational, with three-step functions as trial functions, interactive i m p r o v e m e n t , result fitted to five-term polynomial
" T h e results are s u m m a r i z e d in Tables I1-V.
182
CHOW,
FRIEDSON,
AND YUNG
TABLE 11 REFLECTION (R) A N D T R A N S M I S S I O N ( i t ) F U N C T I O N S FOR A H O M O G E N E O U S A I M O S P H E R E W I T H 71 = I A N D V A R I O U S V A L U E S OF (oil A N D /~ll (CASES A - F )
/,
0
0.1
0.3
0,5
11.7
R Ro A T T, A
0.463 0.439 0.023 0.243 0.241 0.001
17 34 83 37 51 86
0.467 0.464 0.003 0.292 0.293 0.001
91 88 03 06 40 34
Case A r t I, ¢o~ 0.436 71 0.380 42 0.436 66 0.380 51 0.000 1)5 0.000 119 0.336 24 0.3211 93 0.336 36 0.320 91 -0.000 12 0.000 02
I, /x~ - I 0,328 82 1/.328 91 -0.000 09 0.290 32 0.290 30 0.000 02
R Ro A T 7~ A
0.377 0.359 0.017 0.187 0.186 0.001
24 37 87 29 14 15
0.374 0.372 0.002 0.223 0.224 0.001
99 67 32 01 03 02
Case B ~-i I, 0.344 31 0.298 0.344 24 0.298 0.000 07 0.000 0.257 48 0.246 0.257 58 0.246 -0.0011 IO 0.000
0.9, #~ 0.256 0.256 0.000 0.223 0.223 0.000
R Ro A T 7]~ A
0.200 0.194 0.006 0.087 0.086 0.000
54 11 43 39 91 48
0.191 0.191 0.000 0.101 0.101 -0.000
96 12 84 43 76 33
Case C r~ I. ¢o~ 0.169 38 0.144 42 0.169 33 0.14443 0.000 (15 -0.000 01 0.117 66 0.113 80 0.117 70 0.113 80 0.000 114 0.000 00
R Ro A T 7], A
0.054 0.053 0.000 0.021 0.020 0.000
10 70 40 08 88 20
0.050 0.050 0.000 0.023 0.023 0.000
03 00 03 61 59 02
Case D TT I, 0.042 48 I).035 0.042 47 0.035 0.000 01 0.000 0.027 38 0.026 0.027 38 0.026 0.000 00 0.000
R Ro A l" l~ A
2.132 2.986 -0.853 0.435 0.138 0.296
60 22 62 83 98 85
1.629 1.779 -0.150 0.272 0.177 0.094
58 64 06 07 54 53
Case E 7i I, (o~ 1.054 87 0.776 19 1.058 44 0.771 82 0.003 57 11.004 37 0.279 32 11.316 92 0.273 71 0.318 63 0.005 61 -0.001 71
R Ro A T To A
1.826 2.614 --0.787 0.360 0.091 0.268
56 00 44 17 56 63
1.377 1.516 --0. 139 0.204 0.116 0.087
19 58 39 00 90 10
¢o~ l0 15 05 68 68 00
0.9
F-
I).286 0.287 -0.000 0.260 0.260 0.000
96 03 07 20 18 02
0.269 32 1/.269 39 -0.000 07 0.246 54 0.246 53 0.0110 01
0.341 0.341 0.000 0.658 0.658 0.000
37 33 04 63 67 04
0.223 0.223 0.000 0.200 0.200 0.000
84 88 04 71 70 01
0.209 0.210 -0.000 0.190 0.190 0.000
96 00 04 27 26 01
0.267 0.267 0.000 0.591 0.591 0.000
46 41 05 59 63 04
0.100 41 0.10042 -0.000 01 0.088 68 0.088 68 0.000 00
(l. 129 0.129 1t.000 0.471 0.471 0.000
54 51 03 36 37 01
I 92 97 05 65 63 02
0.6,/*o I 0.123 57 0.12358 I).000 01 0.103 75 0.103 75 0.000 00
(oo = 0.2, /*o 67 0.030 67 0.030 00 0.000 73 0.024 73 0.024 00 0.000
1.0
I).107 21 0.10722 -0.000 01 0.093 43 1/.093 43 0.000 O0
I 29 29 00 51 51 00
I,/*o - (J. I 0.614 20 11.610 02 0.004 18 0.316 99 0.318 98 0.001 99
Case F rl I, oJo - 0.9, /*o = 0.1 0.87t) 65 0.631 78 0.495 83 0.875 57 0.628 83 0.492 72 --0.004 92 0.002 95 0.1103 11 0.200 39 0.234 86 0.238 97 0.194 25 0.235 79 0.240 40 0.006 14 0.000 93 0.001 43
0.026 1/.026 0.000 0.022 0.022 0.000
17 17 00 15 15 00
0.024 0.024 0.000 0.021 0.021 0.000
47 47 00 05 05 00
0.031 0.031 0.000 0.392 0.392 0.000
98 98 00 26 26 00
0.508 0.504 0.003 0.301 0.303 --0.001
25 86 39 68 24 56
0.467 0.464 I).003 0.292 0.293 -0.001
91 88 03 06 40 34
0.695 0.697 --0.002 0.304 0.302 0.002
36 66 30 64 34 30
1/.408 13 lL405 55 0.002 58 0.229 61 11.230 78 -0.001 17
0.374 0.372 0.002 0.223 0.224 0.001
98 67 31 01 03 02
0.566 0.568 -0.002 0.228 0.225 0.002
011 99 99 34 75 59
N o t e . The variation principle with three-term polynomial is used. The " e x a c t " values Ro and To are taken from van de Hulst 119801. The symbol A denotes the difference between "'exact" and approximate solutions.
RADIATIVE TRANSFER IN ATMOSPHERES All integrals are p e r f o r m e d numerically by Gaussian quadrature. The column F ÷,- in the tables serves as a further c h e c k on our calculations. F + (entered in row R) is defined as
Ez(s)J(s) ds,
F + = 2too
which is the u p w a r d flux. F - (entered in row T) is defined as F - = 2too fo' E2(71 - s)J(s) ds + Fe-~l/~O, which is a n o t h e r expression for the d o w n w a r d flux. Their " e x a c t " values are taken f r o m van de Hulst (1980). Table II shows the results of calculations using a s e c o n d - o r d e r polynomial trial function for an a t m o s p h e r e of optical thickness rl = I. too was chosen to be constant for
183
convenience. Typically the error stays small (a few parts per ten thousand) until/z goes below 0.3. In case E and F, the solar zenith angle /z0 is small and the error is significant as expected. Table III contains results for an atmosphere of larger optical depth, r~ = 4. A fifthorder polynomial trial function was used to achieve c o m p a r a b l e a c c u r a c y as in Table II. Again, the error increases f o r / z -< 0.3. F o r an infinite a t m o s p h e r e , we took the results of H u a n g (1952b) and reconstructed his calculations to eight significant figures in Table IV. We used double precision on an I B M 370/3032 and achieved higher numerical a c c u r a c y than did Huang. H u a n g ' s choice of trial function gave him an accuracy of a few parts in ten million. N o t e also that A grows as /z goes to zero. The fluctuating sign of A for /z --< 0.3 has little significance b e c a u s e the last significant fi-
T A B L E 111 REFLECTION ( R ) AND TRANSMISSION ( T ) FUNCTIONS FOR CASES G - I W I T H 7"I = 4. USING SIx-TERM POLYNOMIAL
0
0.1
0.3
0.5
0.7
0.9
1.0
F ~,-
R R0 A T To A
0.623 0.593 0.030 0.142 0.132 0.010
61 10 51 76 28 48
0.661 0.657 0.004 0.167 0.164 0.002
81 37 44 11 76 35
CaseGrt =4, 0.698 63 0.705 0.698 62 0.705 0.000 01 --0.000 0.216 l0 0.262 0.216 05 0.262 0.000 05 -0.000
to0 = l,/z0 = 1 77 0.697 05 98 0.697 23 21 --0.000 18 18 0.302 22 30 0.302 34 12 -0.000 12
0.679 0.679 --0.000 0.333 0.333 -0.000
31 46 15 42 52 10
0.668 0.668 --0.000 0.345 0.345 -0.000
41 55 14 48 58 10
0.690 0.690 --0.000 0.309 0.309 0.000
92 93 01 07 07 00
R R0 A T To A
0.732 0.698 0.043 0.032 0.027 0.005
97 77 20 44 12 32
0.686 0.682 0.004 0.035 0.033 0.001
99 18 81 27 62 65
Case H ~'1 = 4, w0 = 1,/x0 = 0.5 0.606 59 0.541 56 0.488 84 0.606 52 0.541 72 0.488 99 0.000 07 --0.000 16 -0.000 15 0.045 28 0.058 43 0.073 66 0.045 18 0.058 49 0.073 72 0.000 10 --0.000 06 -0.000 06
0.445 0.445 -0.000 0.088 0.088 --0.000
34 45 11 66 72 06
0.426 0.426 -0.000 0.095 0.095 -0.000
31 41 10 47 52 05
0.504 0.504 0.000 0.071 0.071 0.000
14 09 05 84 82 02
R R0 A T To A
1.327 1.292 0.034 0.071 0.075 -0.004
14 99 15 63 88 25
1.186 1.185 0.000 0.092 0.094 -0.001
15 98 17 75 68 93
0.727 0.727 0.000 0.205 0.205 0.000
60 59 01 67 60 07
0.698 0.698 0.000 0.216 0.216 0.000
63 62 01 10 05 05
0.825 0.825 -0.000 0.174 0.174 -0.000
01 02 01 97 98 01
Note.
Case I r1 = 4, too = 1,/x0 = 0.999 49 0.880 57 0.794 0.999 64 0.880 57 0.794 - 0 . 0 0 0 15 0.000 00 0.000 0.124 73 0.153 55 0.181 0.124 80 0.153 46 0.181 - 0 . 0 0 0 07 0.000 09 0.000
Values of R0 and To are taken from van de Hulst (1980).
0.3 85 84 01 22 14 08
184
CHOW, FRIEDSON, A N D YUNG T A B L E IV
REFLECT|VITY (R) AND TRANSMISSION (T) FOR CONSERVATIVE I N F I N I T E ATMOSPHERE (CASE J)
#
R
0 0.1 0.2 0.3 0.4 0.5 0.6 0,7 0.8 0.9 1.0
0.726 0.824 0.878 0.918 0.949 0.975 0.996 1.015 1.030 1.044 1,056
Ro
9193 3332 6156 4890 8546 4632 8943 1573 9418 7410 9203
0.726 0.824 0.878 0.918 0.949 0.975 0.996 1.015 1.030 1.044 1.056
A
9527 3314 6140 4891 8548 4633 8942 1573 9418 7410 9203
0.000 0.000 0.000 0.000 0.000 0,000 0.000 0.000 0.000 0.000 0,000
(1334 0018 0016 0001 0002 0001 0001 0000 0000 0000 0000
Note. T h e v a l u e s for R are c o m p u t e d by us using the p a r a m e t e r s o b t a i n e d by H u a n g (1952b); t h o s e for R0 are c o m p u t e d u s i n g the results o f P l a c z e k (19471.
-i
IOt
,
I
'
O-IL
I
,
I
I
cose A tOo: I
r I:
I
I
gure is probably s w a m p e d by round-off error. Table V s h o w s the results o f an alternative use o f the functional. Instead o f deriving the source function directly from the functional, S z e ' s m e t h o d (Sze, 19761 o f a N - s t e p approximation ( N = 3 in our case) and then iteration w a s used. The iterated J and K were finally fitted to a fifth-order polynomial for numerical integration. Since the J(t) and K(t) used are not directly derived from the function, R~, T1 and R2, T2 are no longer equal to R3, /2~ (to s e c o n d order). The data demonstrated that the error o f the functional m e t h o d is significantly smaller than that o f direct convolution. The variation o f error with /z for a polynomial trial function is s h o w n graphically in Fig. 2 (for ~'j = 11 and Fig. 3 (for 7l = 4). The dip in the curve at around /z - 0.3 is not a true m i n i m u m , but a zero crossing w h e n 2~ changes sign as p, is decreased. Figures 4 and 5 s h o w the i m p r o v e m e n t o f the functional o v e r direct convolution. We
I p.O= I '~-°l
'
I
1(52"-~ 3-
~
\\,
'
I
I
I
'
cose G OJo=l T,=4
I
'
/~o=1
IR-%i
Izxl \
5
\
x
x
_x
~x~
3
"* ~ f IT_ToI
io-5
IAI
\~\\
IR -%1
~ ,
./'-~
I T # ; [- . . . . . .
10-6 3 10-7
0
J
I
,
0.2
FIG. 2. Variation
I
,
0.4
[
,
0.6
I 0.8
5
i L0
/z of
18RI and 18TI with/z
for case A of
Table I. Both ISRIand 18Tl g o t h r o u g h a c h a n g e in sign at a r o u n d p, = 0.3 a n d t h e r e f o r e the dip in the graph. N u m e r i c a l v a l u e s o f R0 a n d To are t a k e n f r o m v a n de H u l s t (1980) for o d d v a l u e s o f p. a n d f r o m Carlstedt and Mullikin (1966) for e v e n v a l u e s o f p,.
id s 5
id7 0
J
I O.2
,
I
J
0.4
1 0.6
~
I 0.8
/z FIG. 3. S a m e as Fig. 2 e x c e p t 7, = 4.0.
i I.O
RADIATIVE
TRANSFER
IN ATMOSPHERES
185
TABLE V THREE DIFFERENT EXPRESSIONS FOR EVALUATING THE REELECTION AND TRANSMISSION FUNCTIONS AS DEFINED aV EQ. (10) ARE COMPARED 0.1
0.3
0.5
0.483 0.018 0.462 -0.001 0.465 0.000
28 40 95 93 21 33
0.437 0.001 0.433 -0.003 0.436 0.000
84 18 58 08 66 00
TI AI T2 A2 7"3 A3
0.316 0.023 0.296 0.002 0.293 0.000
88 48 34 94 90 50
0.340 0.004 0.336 --0.000 0.336 0.000
53 17 11 25 45 09
0.321 0.000 0.319 --0.001 0.320 0.000
Ri A1
1.728 --0.051 1.728 --0.051 1.769 -0.009
16 48 16 48 66 98
1.042 --0.015 1.067 0.008 1.056 -0.001
77 67 39 95 98 46
Case L T1 = l, 0.764 --0.007 0.790 0.018 0.771 -0.000
0.209 0.032 0.209 0.032 0.175 -0.002
81 27 81 27 50 04
0.286 0.013 0.309 0.036 0.273 0.000
76 05 72 01 85 14
A2 R3
A3 TI A1 /'2 A2
T3 A3
0.9
Case K r~ = 1, too = 1,/-I-o = 1 0.379 03 0.326 98 -0.001 48 -0.001 93 0.377 79 0.326 59 --0.002 72 --0.002 32 0.380 48 0.328 87 - 0 . 0 0 0 03 - 0 . 0 0 0 04
Rt A1 RE A2 R3 A3
R2
0.7
0.326 0.007 0.351 0.032 0.319 0.000
1.0
0.285 -0.001 0.285 -0.001 0.286 -0.000
12 91 04 99 99 04
0.267 --0.001 0.267 -0.001 0.269 --0.000
54 85 54 85 35 04
49 81 07 23 29 01
0.259 -0.001 0.258 0.001 0.260 -0.000
04 14 95 23 17 01
0.245 -0.001 0.245 --0.001 0.246 -0.000
32 21 32 21 51 02
~o = l, t~o = 0. l 31 0.605 88 51 --0.004 14 53 0.629 91 71 0.019 89 66 0.610 23 16 0.000 21
0.502 -0.002 0.523 0.019 0.505 0.000
40 46 90 04 18 32
0.462 -0.001 0.483 0.018 0.465 0.000
95 93 28 40 21 33
0.306 0.003 0.328 0.025 0.303 0.000
68 44 28 04 75 51
0.296 0.002 0.316 0.023 0.293 0.000
34 94 88 48 89 49
20 29 86 05 93 02
0.289 -0.000 0.289 --0.001 0.290 -0.000
19 56 48 85 13 50
0.323 0.004 0.347 0.028 0.319 0.000
90 92 65 67 52 54
Note. T h e trial functions are obtained in three steps. First the variation principle using three-step functions is used. T h e result is refined by one iteration (Sze, 1976) using the integral equation, and finally fitted to a fifth-order polynomial. started with a three-step approximation and the numerical integration was done using a fifth-order polynomial as described earlier. A stands for the difference as in Table II, its improvement ranges from a factor of 2 to two orders of magnitude. Another case o f interest is that of a Lambert surface (with albedo h) at the lower boundary of the atmosphere. In this case, the functions S ( s , / z ) a n d G(t,s) i n E q . (5) are modified
S (s ,~ ) =
to
~x/-L~o~ s)
[ F e-s/~' + --~ hF e-'l#'OEz('rl - s) ]
x/co0(t)o~0(s)
G(t,s) =
2
[Ellt Otherwise, through
s I + 2XE2(rj the
as before.
similar degree
entire
-
t)Ez(rl -
calculation
We obtained
of accuracy
s)]. goes
results with
as in the case h
=0.
CONCLUSION
We have provided a generalization of the method of Stokes and DeMarcus (1971) to calculate reflectivity and transmissivity at arbitrary emergent angles. For the special
186
CHOW, FRIEDSON, AND YUNG I0'
'
I
'
I
'
I
'
I
3
o 10-2
\
x
case K ~o =I
L =I /~o =I
\%
3
i 0-~
~'---o-""
+
la,-al
1
3
IAI 10- 4 +
+
+~+~+
3
3 Id 6
3
i0-7 0
1
0.2
,
I
0.4
,
I
0.6
J
0.8
i
1.0
FIG. 4. Comparison of accuracy of the three methods of evaluating the reflectivity using the threestep approximate J as in Table V. Note that the functional method can be more accurate by one to two orders of magnitude. case of trial functions linear in the variational p a r a m e t e r s there is a degeneracy between the Stokes and D e M a r c u s method and ours. Error analysis has d e m o n s t r a t e d the difference b e t w e e n the two independent uses of the functional. If one uses trial functions of the form given by (12) and derives the m e a n intensity J from the functional, then the reflectivity and transmissivity can be c o m p u t e d to second order by direct convolution. It is u n n e c e s s a r y to evaluate K or the full functional. I f one uses a J (and K) not directly derived from the functional then the functional gives better results than direct convolution. There is a natural limitation on the a c c u r a c y of such methods in the sense that errors b e c o m e significant for # <~ 0.2. This is to be expected since small # is equivalent to large effective optical depth and the simple trial functions used in the sample calculations are insufficient.
The construction of our functional was motivated by Levine and Schwinger's (1950) variational principles of electromagnetic w a v e diffraction. The theory of radiative transfer is closely related to the theory of neutron diffusion. Therefore it is not surprising to detect striking similarity of our functional with the one given by Davison (1957). The relation between the functional and simple physical quantities, first discovered by Stokes and DeMarcus (1971), has been extended by us. Although this p a p e r does not significantly advance the computational aspects of the planetary problem beyond Sze (1976), it has brought new i n s i g h t - - t h e connection between the principle of reciprocity and the variational m e t h o d - - t o the problem of radiative transfer. I n d e p e n d e n t work by C h e y n e y and Arking (1976) and Shia and Yung (personal communication, 1984) indicates that generalization to problems with nontrivial phase function and non-plane-parallel g e o m e t r y can be made.
id' '
3
1
'
elK
'~Oo::
T,'=I I/.,o='1
]
~52 3
"',
i0-~
"%
..... I#_~:#J.....
I
5
IAI
10-4 3
3 1~56 3 i0 7
0
I
]
0.2
I
I
0.4
L
I
0.6
~
I
0.8
FIG. 5. Same as Fig. 4, for transmissivity.
1,0
RADIATIVE TRANSFER ACKNOWLEDGMENTS One of us (Y.L.Y.) acknowledges helpful discussions with N. D. Sze and R. M. Goody. This research was supported by NASA Grant NSG 7376 to California Institute of Technology.
REFERENCES
CARLSTEDT, J. L., AND T. W. MULLIKIN (1966). Chandrasekhar's X- and Y-functions. Astrophys. J.
Suppl. Ser. 113(12), 449. CHANDRASEKHAR, S. (1960). Radiative Transfer. Dover, New York. CHEYNEY, H., III, AND A. ARKING (1976). A new formulation for anisotropic radiative transfer problems: I. Solution with a variational technique. Astrophys. J. 207, 808. DAVlSON, B. (1957). Neutron Transport Theory. Oxford Univ. Press, London/New York. GRAGG, W. B. (1972). The Pade Table and its relation to certain algorithms of numerical analysis. SLAM. Review, 14, 1. VAN DE HULST, H. C. (1980). Multiple Light Scattering, Vol. I. Academic Press, New York.
IN ATMOSPHERES HUANG, blems with a HUANG, blems
187
S. S. (1952a). The variational method for proof radiative transfer: I. Isotropic scattering constant net flux. Astrophys. J. 117, 211. S. S. (1952b). The variational method for proof radiative transfer: III. Reflection effect. Astrophys. J. 117, 221. LEVINE,H., AND J. SCHWINGER (1950). On the theory of electromagnetic wave diffraction by an aperture in an infinite plane conducting screen. Commun. Pure Appl. Math. 3, 355. MINNAERT, M. (1941). The reciprocity principle in lunar photometry. Astrophys. J. 93, 403. PLACZEK, G. (1947). The angular distribution of neutrons emerging from a plane surface. Phys. Rev. 72, 556. STOKES, R. A., AND W. C. DEMARCUS (1971). A simple technique for calculating line profiles of inhomogeneous planetary atmospheres. Icarus 14, 307. SZE, N. D. (1976). Variational methods in radiative transfer problems. J. Quant. Spectrosc. Radiat. Transfer 16, 763. YUNG, Y. L. (1976). A numerical method for calculating the mean intensity in an inhomogeneous Rayleigh scattering atmosphere. J. Quant. Spectrosc. Radiat. Transfer 16, 755.