Introduction
In a previous paper by one of the authors (Reference BarkstromBarkstrom, 1972, hereafter referred to as paper I) it was shown that multiple scattering in snow could be modelled with considerable success by solving the phenomenological equation of transfer with isotropic scattering. However, as Reference Middleton and MungallMiddleton and Mungall (1952) and Reference Salomonson and MarlattSalomonson and Marlatt (1968) have shown, the light “reflected” from snow is directed in the forward scattering direction, particularly for low solar altitudes. This phenomenon cannot be explained by isotropic scattering. The purpose of this paper is to refine the results in I by modelling the optical properties of snow including the effects of anisotropic scattering in finite layers. In this model the transport equation is solved by an algorithm that gives the optical properties of a multiply scattering layer composed of two stacked layers of known properties. With this algorithm it is possible to infer a phase function that will result in a layer whose surface albedo and angular reflectance is in good accord with experimental measurements. The modelling procedure is sufficiently general that vertically inhomogeneous layers can be constructed which closely simulate stratified layers with different scattering properties. Moreover, the reflectance properties of the “ground” underlying the “snow” layer can be realistically included.
Numerical results of the modelling procedure are given that compare well with Reference Middleton and MungallMiddleton and Mungall’s (1952) measurements on “new snow fallen in calm”. The computations indicate that inclusion of anisotropy causes the flux divergence to maximize deeper in the snow than for isotropic scattering (in terms of optical depth) and that layers of optical depth less than five will exhibit a flux “extinction coefficient” three to five times that observed deep within “semi-infinite” layers (for the same extinction coefficient). The success of the anisotropically scattering model in predicting layer albedo and angular surface reflectance for the above snow suggests that it would be both possible and useful to characterize other types of snow layers by “effective” phase functions. Although no simple analytic procedures are found for deducing these phase functions, they are readily formed by comparing observations of surface albedo and angular reflectance with the results of calculations. A physical theory for the phase function will be found elsewhere (Reference Bohren and BarkstromBohren and Barkstrom, 1974).
Formulation of the Problem
In paper I the equation of transfer, equation (I-11) was solved with an isotropic phase function, equation (I-14), to obtain an analytic expression for specific intensity (spectral radiance) I(τ, μ, φ) within and at the surface of a vertically semi-infinite medium. In this paper we use an alternate, but equivalent, method of solving the transfer equation to obtain the specific intensities within and emerging from the top and bottom of a laterally infinite, but vertically finite, plane-parallel slab containing anisotropic scatterers. The method of solution is sufficiently general that a wide variety of scatterers can be incorporated in the slab, so that realistic, vertically inhomogeneous slabs such as snow packs can be modelled fairly easily. The solution algorithm is well known in the study of stellar and planetary atmospheres, but because it is less well known elsewhere and because it appears to be a useful tool in characterizing the optical properties of snow, the algorithm is discussed in some detail in the text and in the Appendix.
As in paper I, I(τ, μ, φ) is the specific intensity of radiation in a laterally infinite, plane-parallel slab at an optical depth τ travelling in the direction specified by the azimuth φ and the cosine of the polar angle μ = cos θ where θ is measured from the downward normal to the slab. The optical depth is the dimensionless product τ = χd of the volume extinction coefficient χ and the vertical distance d. It is measured downward from the top of the slab. It is convenient to set μ = |cos θ| and explicitly prefix μ with the proper sign so that +μ denotes radiation travelling downward and –μ denotes radiation travelling upward. It is usual (but not necessary) to assume that no sources are present within the slab so that the only radiation sources are the intensities I(0, +μ, φ) and Ι(τ 0, —μ, φ) respectively incident at the top and bottom of a slab of total optical depth τ 0.
The direct problem is to find the intensities I(τ,±μ, φ) within the slab and the emergent intensities I(0, —μ, φ) and I(τ 0, +μ, φ) given the properties of the scatterers and the incident intensities. The algorithm used to find the internal and emergent intensities is variously known as the method of addition of layers (Reference SobolevSobolev, 1963), and doubling method (Reference Van de HulstVan de Hulst, 1963), and the star-product algorithm (Redheffer, 1963; Reference Grant and HuntGrant and Hunt, 1969[a],[b]), but it is formally identical to that given by Reference StokesStokes (1862) for the transmission and reflection from a pile of glass plates. The essence of the method is that if the reflection and transmission properties of one or more scattering layers are known, then the same properties of a stack of two such layers can be found in terms of those of the individual layers. Thus it is possible to compute the properties of a multiply scattering layer by repeated doubling of an initial, optically thin, single-scattering layer of known properties. It is also possible to combine successively layers of dissimilar properties to obtain a vertically inhomogeneous stack with known properties. The algorithm also makes it possible to add arbitrary reflecting boundaries within and at the surface of the slab.
The related inverse problem is to infer the properties of the scatterers given the incident and emergent intensities. The solution to the inverse problem is not as straightforward as that of the direct problem. At present there is no well established procedure for inferring scattering properties when intensities incident upon and emergent from a scattering slab are given. The numerical example given later used an iterative procedure which proceeded from initial estimates of the single scattering properties of the slab with isotropic scattering to a model which produced emergent intensities which compared well with experimental measurements. Fortunately the procedure converged quickly to the results given below.
Method of solution: The Addition Algorithm
The addition algorithm assumes that there exist reflection and transmission operators which linearly transform intensities incident from above and below into intensities emerging from the slab. Assume that the slab is of optical depth τ, and that it is illuminated by either collimated or diffuse radiation from above and below. For brevity let the intensity incident on the top of the slab I(0, +μ, φ) be given by I 0↓ where the down arrow denotes the direction of travel and the subscript denotes the level τ = 0. Let the intensity incident from below be I 1↑ where the subscript denotes τ = τ I,. The emergent intensities are
where R, T, R ∗ and T∗ are the reflection and transmission operators. When the slab is homogeneous or symmetric about the median plane, R = R * and T = T ∗. When the slab is inhomogeneous or when polarization properties are included, R ≠ R ∗ and T = T ∗. In the text which follows, but not in the appendix, it will be assumed that R ∗ = R and T ∗ = T.
Suppose that reflection and transmission operators R 1, and T 1, are known for a slab of optical depth τ 1 and that the corresponding operators R 2 and T 2 are known for a slab of optical depth τ 2. If the two slabs are stacked with the slab of optical depth τ, on top, there will be operators R and T for the stack. In terms of these operators the intensities emerging from the stack are
where the subscript 2 denotes the bottom of the stack. These intensities can also be related to the intensities at the (fictitious) interface between the top and bottom slabs
The intensities at the interface are related to one another and to I 0↓ by
The equations may be solved for explicit expressions for the interface intensities
The inverse operator [I —R 1 R 2] is assumed to exis T, and, as is made clear in the appendix, it may be represented by a simple matrix inversion. Substitution of Equations (9) and (10) into Equations (5) and (6) and comparison with Equations (3) and (4) gives the expressions for R and T:
The R and T operators also may be defined in terms of analogous scattering and transmission functions S + and T + given by Reference ChandrasekharChandrasekhar (1950). For diffuse incident radiation they are defined by
These operators give only the diffuse emergent radiation. Accordingly a term exp (—τ 0/μ) I(0, μ ’, φ ’) δ μ μ , δ φ φ , must be added to Equation (14) to give that portion of I(0, μ, φ') already travelling in the (+μ, φ) direction which is transmitted and attenuated by the slab. The Kronecker deltas δ μ μ , δ φ , serve as a reminder that the term is added only for the μ = μ ’, φ = φ ’ component of I(0, μ ’, φ ’). The definitions for R and T thus become
The analogy between Equations (15) and (16) and Equations (3) and (4) becomes somewhat clearer when the incident radiation is collimated. In that case the incident intensity in the direction (+μ 0, φ 0) is
where φ 0 is the net flux (spectral irradiance) and δ(μ—μ 0) and δ(φ—φ 0) are Dirac delta functions. The emergent intensities are
and the expression for R and T, Equations (15) and (16) remain the same without the integrals.
Although the fundamental procedure for finding R and T from given R 1, R 2, T 1, and T 2 has been given, explicit equations for the operators R 1, and T, or R 2 and T 2 have not been given. Initial values for R and T operators may be obtained in either of two ways. Since the R and T operators are simply related to the single-scattering phase function in the limit τ → 0. optically thin layers (τ ≈ 10-6), can be doubled repetitively to obtain R and T operators for thicker layers. Alternatively, integro-differential equations given by Reference ChandrasekharChandrasekhar (1950) can be solved numerically with reasonable ease to reach depths of τ≈ 10−3 in one integration step with doubling or addition used thereafter. The method used is partly a matter of taste, but should more properly be dictated by the computational economics of the problem to be solved.
Starting values for R and T are doubled from the phase function p using Reference ChandrasekharChandrasekhar’s (1950) small τ limits
These equations reduce to
for small τ 0. The factor Ѓ F(τ 0) is
The procedure for converting S + and T + functions for small τ 0 to the corresponding R and T is evident from Equations (15) and (16), but it is discussed in detail in the appendix. The discussion of the initialization via integro-differential equations is also given in the appendix.
The effects of the albedo of the ground underlying a slab of snow can also be accommodated in the addition algorithm by treating the reflecting layer as if it were a scattering layer. For example, the scattering function S + for a Lambert surface (a Lambert surface has a diffuse reflectance which is equal in all directions) with albedo λ is 4μμ 0 λ so that the reflection operator is μ 0 λ/π. Fresnel reflections a T, for example, an ice-snow interface can be included by using the usual Fresnel reflection factors with the appropriate change in the sign of μ, i.e. μ → — μ upon reflection. Since reflections are not included in the numerical example no further mention of them is made.
Reflectance of a snow cover
Reference Middleton and MungallMiddleton and Mungall (1952) have given measurements of the angular reflectance for several types of snow layers. The measurements share the property that the layer reflectance increasingly diverges from that of a Lambert surface with decreasing solar elevation angles. This is consistent with the anisotropy present in scattering of light by large particles. We suggest that this behavior can be successfully modelled by multiple scattering in a plane parallel layer containing scatterers with an effective (or fictitious) single-scattering albedo
The single-scattering albedo is the fraction of photons incident on a particle which are scattered rather than absorbed. The angular dependence of the scattering is given by the phase function p (cos θ) which explicitly depends on the scattering angle θ. When the phase function has the normalization
the product mp(cos Θ) dΩ/4π is the probability that a photon interacting with the particle will be scattered through an angle θ into an element of solid angle dΩ = d(cos Θ) dø. The anisotropy of the scattering is characterized by the asymmetry parameter g given by
It is convenient to expand the phase function in a Legendre series
whose coefficients are
in which the Pi(cos Θ) are Legendre polynomials. Equations (25)-(28) give f0 ≡ 1 and g ≡ f1/3.
We have attempted to model the angular reflectance data of Reference Middleton and MungallMiddleton and Mungall (1952) shown replotted in Figure 1 for “new snow fallen in calm”. The single-scattering albedo largely determines the level of the curves and is expected (see paper I) to be about 0.99 for isotropic scattering. The shape of the curves is determined by the phase function, with g playing the major role. The angular reflectance for near grazing incidence and reflection (μ, μ 0 ≈ 0) is sensitive to the higher terms. However, the properties important for the energy balance, the albedo and the flux extinction, are extremely insensitive to f 2, f 3., ... (Bohren and Barkstrom, in press; Reference Van de HulstVan de Huls T, 1970). Accordingly, we have searched for
values of
Fig. 1. Plot of the angular dependent of the reflected intensity I from snow new-fallen in calm weather. Data are taken from Reference Middleton and MungallMiddleton and Mungall (1952). The azimuthal angle is denoted by ø, with ø = 0° being the direction in which the solar radiation is travelling, ø = 180° being the direction in which light back-scattered from the surface is travelling, μ is the sine oftht altitude at which Ike light is travelling or, equivalently, the cosine of the zenith angle, μ 0 denotes the sine of the solar altitude.
Our early attempts to reproduce Figure 1 used a Henyey-Greenstein phase function (see for example Reference Van de HulstVan de Huls T, 1970)
which has an infinite Legendre expansion with coefficients
which depend only on l and g. For values of g greater than 0.4 or 0.5 these coefficients decrease very slowly with increasing l and it is necessary to retain many terms in the Legendre expansion of Equation (29) to avoid truncation problems. This makes it necessary to use large matrices for R and T (see the Appendix) to preserve numerical stability and accuracy. It is also possible to ignore truncation problems and use too small an order of quadrature to estimate an appropriate value of g. This procedure suggested that
Fig. 2. Plot of the phase function p for single scattering, as a function of cos θ, where θ is the angle through which the light is scattered.
The alternate phase function, shown in Figure 2, has g = 0.506 and was found by a Mie calculation for a polydispersion of spheres (Reference Querfeld, C.W., M. and J. P.Querfeld and others, 1972) distributed according to a zeroth order logarithmic distribution (Reference KerkerKerker, 1969, p. 356) with modal diameter 0.245 μm' width parameter σ0 = 0.07, and real refractive index n = 1.20 at an incident wavelength of 0.62 μm. The significant Legendre coefficients for this phase function are g.iven in Table I. This phase function with
Table I. Legenfire expansion coefficients of the phase function
The emergent upward intensities in Figure 3 can be integrated (see the Appendix) to obtain a surface albedo A(μ 0). Figure 4 shows the albedo A(μ α) computed with the phase function in Table I and Figure 2 using
Fig. 4. Plot of clear-sky albedo A versus the sine of the solar altitude μ 0. Data are from Antarctic expeditions (Reference LiljequistLiljequis T, 1956; Reference RusinRusin, 1961), as noted in paper I. Solid line is the albedo from a very thick layer of anisotropic scatterer with the phase function shown in Figure 2; dashed line is the albedo from a semi-infinite layer of isotropic scattertr.
Behavior of the Radiation within the snow Layer
It was shown in paper I that the flux divergence in an isotropic scattering atmosphere has a maximum beneath the surface for certain angles of incidence. The flux divergence in an “atmosphere” with a forward-scattering phase function, such as snow, exhibits a similar behavior, except that the maximum of the flux divergence lies deeper in the snow than the maximum in the isotropic scattering case when the distance is measured in photon mean free paths (i.e. in terms of optical depth τ = χz). This behavior is shown in Figure 5. Again, we see that for small solar altitudes, the Sun heats the top few centimeters of the snow with respect to the snow further down in the layer. However, the flux divergence is much flatter than was the flux divergence for isotropic scattering.
Fig. 5. Plot of the flux divergence —d{ΦΦ0}/dτ versus the optical depth τ for various solar altitudes, using the anisotropic scatterer with a phase function shown in Figure 2
We also observe that the non-exponential behavior of the flux extends deeper into the snow than it did in the case of isotropic scattering. Using Figure 5, we find that the flux divergence becomes exponential below τ≈ 5, Thereafter
where we measure va ’ to be 13.1 from the g.raph. This in turn implies that the extinction coefficient is now x= 190 m-1 so that τ = 1 corresponds to about 0.526 cm. Since τ = 1 corresponds roughly to the optical depth at which one can distinguish objects through a layer, the new value of the extinction coefficient g.ives a more accurate idea of the thickness a layer of snow must have to be “opaque” to an underlying surface. The estimate of the physical depth corresponding to τ = 1 was incorrectly g.iven in paper I as 0.851 cm, and should have been x−1 = (0.851 cm−1)−1 = 1.175 cm. This would require a rather thick layer to blanket the g.round so that it could not be seen. The new value of χ appears to bring the theoretical predictions of dΦ/dz into better agreement with the observations of Reference Ambach and MockerAmbach and Mocker (1959). The volume melting rate g.iven in Equation (42) of the first paper is changed as a result of anisotropic scattering to 0.569 kg cm−1 S −1. However, the melting rate per unit area is unchanged. We have also calculated the asymptotic ratio of upward to downward flux deep in the medium. We find φ↑as|φ↓as = 0.815 for the phase function g.iven by Table I.
Fig. 6. Plot of the angular dependence of the intensity of light I transmitted through slabs of various optical depths τ for neat normal incidence (μ 0 = 0.995) on the top of the slab.
Transmission of Light through a Finite slab of snow
The transmission of light through a finite layer of snow is of considerable importance because it is one of the ways in which the “extinction coefficient” of snow is determined. Figure 6 shows the intensity transmitted by a thin layer of snow under nearly normal incidence for the phase function g.iven by Table I. It is clear that thin layers of snow (0.025 m thick or less) may g.ive an entirely erroneous impression of the character of the flux extinction. Such an effect was reported by Reference Giddings and LaChapelleGiddings and LaChapelle (1961), and indeed, was predicted by them on the basis of a diffusion theory for radiation. The steep portion of the curve of extinction coefficient versus optical depth is caused by the extinction of light transmitted directly through the slab without being scattered. It is clear that as the slab becomes thicker, the relative contrast between the directly transmitted light and the diffusely transmitted light decreases, until by τ = 7, only the diffuse light is coming through the slab. The strong directionality of the transmitted light for slabs up to about 3 cm thickness means that one may have to make special care in measuring the flux extinction, since detectors are often strongly direction dependent (Reference GillhamGillham, 1970).
Conclusions
We have extended the computations done in paper I to include anisotropic scattering by finite layers of snow. In doing so, the single scattering albedo and the phase function were adjusted until the calculated angular reflectance of a deep snow layer g.ives reasonable agreement with the observed reflectances. It has been shown that the flux divergence may be expected to peak slightly further beneath the surface than was the case for isotropic scattering (as measured in photon mean free paths). Finally, it has been suggested that measurements of the flux extinction coefficient using the light transmitted through finite layers of snow (under 2.5 cm thick) may lead to serious errors if extrapolated to effectively semi-infinite layers below the boundary layer in the latter.
Acknowledgements
The authors wish to thank Dr Lewis House and Dr Thomas Schlatter of NCAR for a critical reading of this paper. One of us (B.R.B.) would also like to g.ive especial thanks to Dr Thomas Schlatter of NCAR, Dr Kevin Pang of LASP, and Dr Roger Berry of INSTAR for their help in suggesting useful and pertinent references.
MS. received 20 August 1973 and in revised form 2 July 1974
Appendix
In the text the addition algorithm and its initialization were discussed in relatively g.eneral terms, although enough information was included for the preparation of numerical codes. In this appendix we discuss the algorithm and initialization more fully and g.ive the complete reduction of the integral reflection and transmission operators to their equivalent matrix forms. In addition we include procedures for calculating intensities and flux divergences within the slab, as well as both the surface and the Bond albedo of the slab. As in the tex T, polarization is omitted, but that omission should be of little importance in this context.
In the text the term specific intensity or spectral radiance is used without explicitly mentioning that the intensity (units: Wm−2 S −1 Hz−1, for example) is measured with respect to a unit area normal to the direction of propagation. In optics it is frequently customary to measure the radiance of a surface with respect to a unit area parallel to the surface (we denote such a radiance by I). This is the reason for the statement in the text that a Lambert surface radiates isotropically. Often a Lambert surface is thought to show a cosine intensity dependence
The use of the reflection operator for μ 0/π for a Lambert surface illuminated by intensity I 0 (per unit area normal to the incoming direction of propagation) g.ives
because of the different definition of the intensity. If I denotes an intensity measured per unit area parallel to a surface and I per unit area normal to itself, the two are related by
Similar problems arise in the use of the term net flux. This is the spectral irradiance or “exitance” found in optics and is related to the specific intensity by
where μ=1 is the normal to the area under discussion. This definition and that of the specific intensity is g.iven in the hope that it will reduce confusion of the discussion of flux divergence and surface albedo which follows.
In the text the intensities emerging from a scattering slab illuminated from above and below were g.iven as Equations (1) and (2)
As was stated in the tex T, these equations are merely an abbreviated form of the full equations
It is the intension of this section to show that the two forms are numerically equivalent if we interpret the intensities I 0↑, I 0↓ and I 1↑, I 1↓ as column vectors (here written as row vectors, but with curly brackets to denote column vectors) of the form
where the μ 1, μ 2, …, μ n and φ 0, φ 1,…,φ n form a discrete g.rid μ and φ. As m and n = ∞, we expect to recover exactly Equations (A.7)-(A.8). It might appear that I 0 should be interpreted as a rectangular matrix with rows indexed on μand columns indexed on φ, but it will become evident that this can be avoided. The operators R, R∗, T, and T ∗ are evidently square matrices whose rows are indexed on the exit angle μ and whose columns are indexed on the incident angle μ ’ or μ ’’.
As the first step in the reduction of Equations (A.7) and (A.8) to Equations (A.5) and (A.6) we recall from Equations (22)-(24) that the single-scattering (small τ) forms for Chandresekhar’s S + and T + functions were proportional to the phase function ρ(μ, φ; μ 0, φ 0)• In Equation (27) the phase function was g.iven as a Legendre series
The scattering angle θ is related to the slab coordinates through the relation
The addition theorem of spherical harmonics (Reference HobsonHobson, 1955) g.ives the replacement
when the associated Legendre functions Pi m(μ) have the normalization
The numerical implication of this normalization will be discussed below. The substitution of Equation (A-12) into Equation (A, 10) and the interchange of the two summations g.ives
in which
The superscript on ρm(μ,μ 0) should not be taken as an exponent. Equations (A.14)-(A.17) provide the means for computing small τ values of R and T in terms of the single-scattering phase function and suggest that R and T also be expressed as a Fourier series.
The Fourier expansions for the Reflection and transmission functions are
and
with analogous equations for R ∗ and T ∗ when they are needed. The utility of these expansions becomes more apparent if we Restate Equations (II) and (12) which g.ive the reflection and transmission functions
These expressions contain products of the form T 1 ∗ R 2 and R 1 ∗ R 1 ∗ which in g.eneral appear as follows:
The result of the Fourier expansion is that azimuthal integrals in the products can be evaluated explicitly and moreover that cross-products of Fourier terms vanish because of the orthogonality of the cosines. Products of operators thus have the same form as the operators themselves. This also is true of the inverse operator since it can be shown both formally and physically that the inverse operator can be expanded as
with an infinite number of terms. The conventional doubling procedure (Reference HansenHansen, 1969; Reference Van de HulstVan de Huls T, 1963) evaluates the right-hand side of Equation (A.27) explicitly by calculating powers of R 1 ∗ until the ratio of successive terms becomes constant and then the remainder of the terms are summed as a g.eometric series. The star-product algorithm calculates the inverse directly and thereby achieves g.reat computational economy.
The μ ’ ntegral in Equation (A.26) must be evaluated numerically as part of the calculation. This is customarily done by using G.auss-Legendre quadrature which permits the exact replacement of the integral by a sum provided that the i ntegrand meets certain requirements. The quadrature formula (Reference Abramowitz and StegunAbramowitz and Stegun,. 1964) is
with
The Wi are quadrature weights and the xi, are abscissas which are tabulated (Reference Abramowitz and StegunAbramowitz and Stegun, 1964, for example). The yi are shifted abscissas which lie in the interval (a, b) rather than in the interval (–1, 1) in which the xi are g.iven. The sum will be exact as long as the integrand (y) can be represented as a polynomial of degree less than 2L. Computationally this is fulfilled by requiring that the order of quadrature l be g.reater than 2N where n is the number of terms in the phase function (A. 10) and in the operators R and T,
The use of the quadrature formula transforms Equation {A.26) into
It should be apparent that the sum in parentheses is nothing more than a matrix product since both incoming and outgoing values of μ in T ∗ and R assume only discrete values which are determined by the quadrature abscissas. If we label incident angles with a j subscript and exit angles with an i subscript then the (i,j)th element of the matrix T 1 ∗ R 2 will be
It should also be evident that the decoupling of the Fourier terms in Equation (A.26) and therefore in Equations (A.20) and (A.21) allows the computation of each Fourier component separately. Thus if we omit the common cosine factors in Equations (A.20) and (A.21) we have the following expressions for Equations (A.31), (A.20) and (A.21)
The corresponding matrix equations for ths mth Fourier matrices for R∗ and T ∗ are
The inverse operators in Equations (A.33)-(A.36) may be obtained by any matrix inversion method.
The initialization of the R and T functions (or matrices) can be accomplished using Equations (22)-(24) for S + and T + and the Fourier expansion for the phase function, Equations (A.14)-(A,16). The expressions for the mall τ 0 (τ 0 ≈ 220) values of the R and T Fourier coefficients are
where we have used the replacement σi = 1/μ i. As usual δ is the Kronecker delta so that the exponential term occurs only on the main diagonal of the T matrix. The factor (- 1)l+m results from p1 m(-μ i) = (-1)i+m) Pl m(μ) which occurs in Equation (22). Matrix elements for larger values of τ 0 are obtained by doubling or addition using Equation (A.33)-(A.36). The Legendre functions in Equations (A.37) and (A.38) are computed by forward recursion from the Legendre polynomials (m = o)
and the recursion relations
and
These recursion relations differ from those usually found for associated Legendre functions because of the normalization (A. 13).
An alternate method for finding initial values of R and T is to integrate integro-differential equations for S + and T + g.iven by Reference ChandrasekharChandrasekhar (1950, p. 169). The equations for the derivatives are
These equations are easily integrated through one or more steps ▵ τ ≤ 1/2 min μ 0 from τ 0 = о using the initial conditions Sij +m = T +m = 0 and the usual fourth-order Runge-Kutta integration scheme. The constraint on ▵ τ insures numerical stability. Although this initialization may appear more cumbersome than doubling thin layers, it may well be faster since depths of τ 0 = 0.0025 (for N=16) are reached in one step. The R and T functions are related to the S + and T + functions through the equations
The reflection and transmission functions used in Equations (A.7) and (A.8) are the Fourier sums (A. 18) and (A.19) whose coefficients are the R ij m and T ij m. These coefficients contain multiplicative factors and an additive term which requires some modification prior to constructing the Fourier sum. The diagonal term exp (—τ 0σi)δij in the T ij m should be subtracted before the summation. The term is needed as a source in the calculation of terms with m ≠ 0, but is included in the Fourier sum only in the m = 0 term. Indeed the term may be omitted from the m = 0 term when the direc T, attenuated, beam is of no interest. After subtraction of the diagonal term (from all m) the R ij m and T ij m coefficients should be multiplied by (E0 m/πWj to obtain properly normalized coefficients which when summed in Equations (A. 18) and (A. 19) will be normalized properly for use in the integrals (A.7) and (A.8), but not necessarily so for use in the matrix form (A.5) and (A.6). When the incident beam is collimated, the R and T functions so obtained are properly normalized, but when the incident beam is diffuse with a distribution both in μ 0and φ 0, ,the integrals must be evaluated with some suitable (e.g. G.auss-Legendre) quadrature scheme and the appropriate weights and normalizing factors from Equations (A.28) and (A.29) must be reinserted in the R and T functions. When G.auss Legendre quadrature is used for either the μ or Φ integration or both, the problems revert to the discrete matrix form. When the φ dependence of the incident beam is suitable, the φ integration can be done analytically using the Fourier expansion (A.18) and (A.19) and the integration becomes the more usual product of R and T matrices indexed on incoming and outgoing μs and the incident beam in column vector form.
The reflection function can be used to compute a surface angular albedo A(μ 0) which g.ives the fraction of light reflected from a surface illuminated by a collimated incident beam I(o, +μ 0,φ 0). Let the incident beam be Φ0δ(μ ’ — μ 0) so that its net flux (measured per unit area parallel to the surface) is
The net flux reflected from the surface is
After performing the integrations over the incident beam, the surface angular albedo is
The Fourier expansion (A.18) enables an explicit integration over ф,
Gauss-Legendre quadrature is used to evaluate the μ integral
The albedo is thus just a weighted sum over the azimuth independent Fourier coefficient R ij 0. The Bond albedo A, the ratio of the net flux reflected by a spherical surface to that incident in a collimated beam, is simply related to the surface angular albedo by
The Bond albedo is also the albedo under isotropic incident radiation.
In the text the run of flux divergence and thus of heating rate within the slab is g.iven. Reference ChandrasekharChandrasekhar (1950) defines the flux divergence in and normal to a plane parallel layer as
The flux divergence can be calculated within a layer of depth τ 2 at a depth τ 1, by computing reflection and transmission functions for layers of depth τ 1, and τ 2—τ 1 the intensities I 1↑ and I 1↓ as in Equations (9) and (10), and the integral in (A.54) using the value of
It is hoped that this Appendix is sufficiently complete that working numerical routines can be prepared directly from it. The routines used here will be included in the Radiative Transfer Library now being assembled at the High Altitude Observatory for use in the computer at the National Center for Atmospheric Research (NCAR), Although the library is designed and optimized for the NGAR computer, some of the routines, including those described here, are reasonably portable. Correspondence regarding availability of numerical codes should be directed to one of us (C.W.Q.).
