Department of Earth System Science zender@uci.edu
University of California Voice: (949) 824-2987
Irvine, CA  92697-3100 Fax: (949) 824-3256

Copyright © 1998–2008, Charles S. Zender
Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. The license is available online at http://www.gnu.org/copyleft/fdl.html.

Facts about FACTs: This document is part of the Freely Available Community Text (fact) project. facts are created, reviewed, and continuously maintained and updated by members of the international academic community communicating with eachother through a well-organized project website. facts are intended to standardize and disseminate our fundamental knowledge of Earth System Sciences in a flexible, adaptive, distributed framework which can evolve to fit the changing needs and technology of the geosciences community. Currently available facts and their urls are listed in Table 1.

Because of its international scope and availability to students of all income levels, the fact project may impact more students, and to a greater depth, than imaginable to before the advent of the Internet. If you are interested in learning more about facts and how you might contribute to or benefit from the project please contact zender@uci.edu.

This monograph on Radiative Transfer provides some core and some supplementary reading material for ESS 200B. We will discuss much of the material in the first twenty pages, and the figures at the end. The Index beginning on page 366 is also helpful.

Yada yada yada.

Contents

List of Figures

List of Tables

1 Introduction

This document describes mathematical and computational considerations pertaining to radiative transfer processes and radiative transfer models of the Earth system. Our approach is to present a detailed derivation of the tools of radiative transfer needed to predict the radiative quantities (irradiance, mean intensity, and heating rates) which drive climate. In so doing we begin with discussion of the intensity field which is the quantity most often measured by satellite remote sensing instruments. Our approach owes much to Bohren and Huffman [1983] (particle scattering), Goody and Yung [1989] (band models), and Thomas and Stamnes [1999] (nomenclature, discrete ordinate methods, general approach). The nomenclature follows these authors where possible. These sections will evolve and differentiate from their original sources as the manuscript takes on the flavor of the researchers who contribute to it.

1.1 Planetary Radiative Equilibrium

The important role that radiation plays in the climate system is perhaps best illustrated by a simple example showing that without atmospheric radiative feedbacks (especially, ironically, the greenhouse effect), our planet’s mean temperature would be well below the freezing point of water. Earth is surrounded by the near vacuum of space so the only way to transport energy to or from the planet is via radiative processes. If E is the thermal energy of the planet, and FASR and FOLR are the absorbed solar radiation and emitted longwave radiation, respectively, then

Earth does not cool to space as a perfect blackbody (41a) of a single temperature and emissivity. Nevertheless the spectrum of thermal radiation FOLR which escapes to space and thus cools Earth does resemble blackbody emission with a characteristic temperature. The effective temperature TE of an object is the temperature of the blackbody which would produce the same irradiance. Inverting the Stefan-Boltzmann Law (73) yields

Substituting (3) and (4) into (2)

1.2 Fundamentals

The fundamental quantity describing the electromagnetic spectrum is frequency, ν. Frequency measures the oscillatory speed of a system, counting the number of oscillations (waves) per unit time. Usually ν is expressed in cycles-per-second, or Hertz. Units of Hertz may be abbreviated Hz, hz, cps, or, as we prefer, s-1. Frequency is intrinsic to the oscillator and does not depend on the medium in which the waves are travelling. The energy carried by a photon is proportional to its frequency

A related quantity, the angular frequency ω measures the rate of change of wave phase in radians per second. Wave phase proceeds through 2π radians in a complete cycle. Thus the frequency and angular frequency are simply related

Most radiative transfer literature use wavelength or wavenumber. Wavelength, λ (m), measures the distance between two adjacent peaks or troughs in the wavefield. The universal relation between wavelength and frequency is

The wavenumber m-1, is exactly the inverse of wavelength

There is another, distinct quantity also called wavenumber. This secondary usage of wavenumber in this text is the traditional measure of spatial wave propagation and is denoted by k.

Table 2 summarizes the relationships between the fundamental parameters which describe wave-like phenomena.

---

Table 2: Wave Parameter Conversion Tablea b Variable ν λ ω k τ Units s-1 m cm-1 s-1 m-1 s ν - 2πν               λ -               100c - 200π               ω -               k ck -               τ cτ -

---

2 Radiative Transfer Equation

2.1 Definitions

2.1.1 Intensity

The fundamental quantity defining the radiation field is the specific intensity of radiation. Specific intensity, also known as radiance, measures the flux of radiant energy transported in a given direction per unit cross sectional area orthogonal to the beam per unit time per unit solid angle per unit frequency (or wavelength, or wavenumber). The units of Iλ are Joule meter^-2 second^-1 steradian^-1 meter^-1. In SI dimensional notation, the units condense to J m-2 s-1 sr-1 m-1. The SI unit of power (1 Watt ≡ 1 Joule per second) is preferred, leading to units of W m-2 sr-1 m-1. Often the specific intensity is expressed in terms of spectral frequency Iν with units W m-2 sr-1 Hz-1 or spectral wavenumber (also I) with units W m-2 sr-1 (cm-1)-1.

Consider light travelling in the direction through the point r. Construct an infinitesimal element of surface area dS intersecting r and orthogonal to . The radiant energy dE crossing dS in time dt in the solid angle dΩ in the frequency range [ν,ν + dν] is related to Iν(r,) by

(13)

It is not convenient to measure the radiant flux across surface orthogonal to , as in (13), when we consider properties of radiation fields with preferred directions. If instead, we measure the intensity orthogonal to an arbitrarily oriented surface element dA with surface normal , then we must alter (13) to account for projection of dS onto dA. If the angle between and is θ then

(14)

and the projection of dS onto dA yields

(15)

so that

(16)

The conceptual advantage that (16) has over (13) is that (16) builds in the geometric factor required to convert to any preferred coordinate system defined by dA and its normal . In practice dA is often chosen to be the local horizon.

The radiation field is a seven-dimensional quantity, depending upon three coordinates in space, one in time, two in angle, and one in frequency. We shall usually indicate the dependence of spectral radiance and irradiance on frequency by using ν as a subscript, as in Iν, in favor of the more explicit, but lengthier, notation I(ν). Three of the dimensions are superfluous to climate models and will be discarded: The time dependence of Iν is a function of the atmospheric state and solar zenith angle and will only be discussed further in those terms, so we shall drop the explicitly dependence on t. We reduce the number of spatial dimensions from three to one by assuming a stratified atmosphere which is horizontally homogeneous and in which physical quantities may vary only in the vertical dimension z. Thus we replace r by z. This approximation is also known as a plane-parallel atmosphere, and comes with at least two important caveats: The first is the neglect of horizontal photon transport which can be important in inhomogeneous cloud and surface environments. The second is the neglect of path length effects at large solar zenith angles which can dramatically affect the mean intensity of the radiation field, and thus the atmospheric photochemistry.

With these assumptions, the intensity is a function only of vertical position and of direction, Iν(z,). Often the optical depth τ (defined below), which increases monotonically with z, is used for the vertical coordinate instead of z. The angular direction of the radiation is specified in terms of the polar angle θ and the azithumal angle ϕ. The polar angle θ is the angle between and the normal surface that defines the coordinate system. The specific intensity of radiation traveling at polar angle θ and azimuthal angle ϕ at optical depth level τ in a plane parallel atmosphere is denoted by Iν(τ,θ,ϕ). Specific intensity is also referred to as intensity.

Further simplification of the intensity field is possible if it meets certain criteria. If Iν is not a function of position (τ), then the field is homogeneous. If Iν is not a function of direction (), then the field is isotropic.

2.1.2 Mean Intensity

The mean intensity is an integrated measure of the radiation field at any point r. Mean intensity ν is defined as the mean value of the intensity field integrated over all angles.

The definition of ν (17) demands the radiation field be integrated over all angles, i.e., over all 4π steradians. Evaluating the denominator demonstrates the properties of angular integrals. The denominator of (17) is

It is convenient to return briefly to the definition of isotropic radiation. Isotropic radiation is, by definition, equal intensity in all directions so that the total emitted radiation is simply 4π times the intensity of emission in any direction.

Applying (19) to (17) yields

Let us simplify (21) by introducing the change of variables

The hemispheric intensities or half-range intensities are simply the up- and downwelling components of which the full intensity is composed

Iν(τ,) = Iν(τ,θ,ϕ)	=	(25a)
Iν(τ,) = Iν(τ,u,ϕ)	=	(25b)

2.1.3 Irradiance

The spectral irradiance Fν measures the radiant energy flux transported through a given surface per unit area per unit time per unit wavelength. Although it is somewhat ambiguous, “flux” is used a synonym for irradiance, and has become deeply embedded in the literature [Madronich, 1987]. Consider a surface orthogonal to the ′ direction. All radiant energy travelling parallel to ′ crosses this surface and thus contributes to the irradiance with 100% efficiency. Energy travelling orthogonal to ′ (and thus parallel to the surface), however, never crosses the surface and does not contribute to the irradiance. In general, the intensity Iν() projects onto the surface with an efficiency cos Θ = ⋅′ , thus

Let us simplify (27) by introducing the change of variables u = cos θ, du = - sin θ dθ. This maps θ [0, π] into u [1,-1]:

The irradiance per unit frequency, Fν, is simply related to the irradiance per unit wavelength, Fλ. The total irradiance over any given frequency range, [ν,ν + dν], say, is Fν dν. The irradiance over the same physical range when expressed in wavelength, [λ,λ - dλ], say, is Fλ dλ. The negative sign is introduced since -dλ increases in the same direction as +dν. Equating the total irradiance over the same region of frequency/wavelength, we obtain

2.1.4 Actinic Flux

A quantity of great importance in photochemistry is the total convergence of radiation at a point. This quantity, called the actinic flux, FJ, determines the availability of photons for photochemical reactions. By definition, the intensity passing through a point P in the direction within the solid angle dΩ is Iν dΩ. We have not multiplied by cos θ since we are interested in the energy passing along (i.e., θ = 0). The energy from all directions passing through P is thus

The usefulness of actinic flux FJ ν becomes apparent only in conjunction with additional, species-dependent data describing the probability of photon absorption, or photo-absorption. Photo-absorption is the process of molecules absorbing photons. Each absorption removes energy (a photon) from the actinic radiation field. The amount of photo-absorption per unit volume is proportional to the number concentration of the absorbing species NA [m-3], the actinic radiation field FJ ν , and the efficiency with with each molecule absorbs photons. This efficiency is called the absorption cross-section, molecular cross section, or simply cross-section. The absorption cross-section is denoted by α and has units of [m-2]. In the literature, however, values of α usually appear in CGS units [cm-2]. To make the frequency-dependence of α explicit we write α(ν). If α depends significantly on temperature, too (as is true for ozone), we must consider α(ν,T).

The probability, per unit time, per unit frequency that a single molecule of species A will absorb a photon with frequency in [ν,ν + dν] is proportional to FJ ν (ν)α(ν)² . Thus α(ν) is the effective cross-sectional area of a molecule for absorption. The absorption cross-section is the ratio between the number of photons (or total energy) absorbed by a molecule to the number (or total energy) per unit area convergent on the molecule. Let Fα ν [W m-3] be the energy absorbed per unit time, per unit frequency, per unit volume of air. Then

Photochemists are interested in the probability of absorbed radiation severing molecular bonds, and thus decomposing species AB into constituent species A and B. Notationally this process may be written in any of the equivalent forms

The probability that a photon absorbed by AB will result in the photodissociation of AB, and the completion of (33), is called the quantum yield or quantum efficiency and is represented by ϕ. As a probability, ϕ is dimensionless³ . In addition to its dependence on ν, ϕ depends on temperature T for some important atmospheric reactions (such as ozone photolysis). We explicitly annotate the T -dependence of ϕ only for pertinent reactions. Measurement of ϕ(ν) for all conditions and reactions of atmospheric interest is an ongoing and important laboratory task.

The specific photolysis rate coefficient for the photodissociation of a species A is the number of photodissociations of A occuring per unit time, per unit volume of air, per unit frequency, per molecule of A. In accord with convention we denote the specific photolysis rate coefficient by Jν. The units of Jν are s-1 Hz-1.

The total photolysis rate coefficient J is obtained by integrating (34) over all frequencies which may contribute to photodissociation

Integration errors due to the discretization of (35) are quite common. To compute J with high accuracy, regular grids must have resolotion of ~1 nm in the ultraviolet, [Madronich, 1989]. Much of the difficulty is due to the steep but opposite gradients of FJ ν and ϕ which occur in the ultraviolet. High frequency features in α worsens this problem for some molecules.

The utility of J has motivated researchers to overcome these computational difficulties by brute force techniques and by clever parameterizations and numerical techniques [Chang et al., 1987; Dahlback and Stamnes, 1991; Toon et al., 1989; Stamnes and Tsay, 1990; Petropavlovskikh, 1995; Landgraf and Crutzen, 1998; Wild et al., 2000]. It is common to refer to photolysis rate coeffients as “J-rates”, and to affix the name of the molecule to specify which individual reaction is pertinent. Another description for J is the first order rate coefficient in photochemical reactions. For example, JNO2 is the first order rate coefficient for

Figure 1 shows the spectral distribution of actinic flux in a clear mid-latitude summer atmosphere, and the absorption cross-section and quantum yield of NO2.

Figure 1: (a) Spectral distribution of actinic flux FJ [# m-2 s-1 μm-1] at TOA and at the surface for a mid-latitude summer (MLS) atmosphere with a unit optical depth of dust or sulfate in the lowest kilometer. (b) Absorption cross section of NO2, αNO2 [m2 molecule-1]. (c) Quantum yield of NO2, ϕNO2 (36).

Figure 2 shows the vertical distribution of JNO2 [s-1] for the conditions shown in Figure (1).

Figure 2: Vertical distribution of JNO2 [s-1] (36) for the conditions shown in Figure (1). (a) Absolute rates. (b) Rates normalized by clean sky rates.

In this section we have assumed the quantities FJ ν , α, and ϕ are somehow known and therefore available to use to compute J. Typically, α and ϕ are considered known quantities since they usually do not vary with time or space. Models may store their values in lookup tables or precompute their contributions to (35). The essence of forward radiative problems is to determine I so that quantities such as Jν and Fα may be determined. In inverse radiative transfer problems which are encountered in much of remote sensing, both J and the species concentration are initially unknown and must be determined. We shall continue describing the methods of forward radiative transfer until we have tools at our disposal to solve for Jν. At that point we shall re-visit the inverse problem.

2.1.5 Actinic Flux Enhancement

The actinic flux FJ ν (31) is sensitive to the angular distribution of radiance Iν. Nearly all scattering processes diffuse the radiation field, i.e., convert collimated photons to more isotropic photons. Such diffusion causes actinic flux enhancement. It is instructive to examine how the relationship between downwelling flux and actinic flux changes in the presence of scattering. Four limiting cases may be identified and are summarized in Table 3.

---

Table 3: Actinic Flux Enhancement by Scatteringa Description F- ν ρL(ν) Iν() F + ν ν FJ ν Collimated, non-reflecting F⊙ ν 0 F⊙ ν δ( -⊙) 0 F⊙ ν               Isotropic, non-reflecting F⊙ ν 0 I- ν = F⊙ ν ∕π I+ ν = 0 0 2F⊙ ν               Collimated, reflecting F⊙ ν 1 I- ν = F⊙ ν δ( -⊙) I+ ν = F⊙ ν ∕π F⊙ ν 3F⊙ ν               Isotropic, reflecting F⊙ ν 1 F⊙ ν ∕π F⊙ ν 4F⊙ ν

---

The scenarios differ in the isotropicity of the downwelling radiance (collimated or isotropic) and the reflectance ρL(ν) (0 or 1) of the lower boundary, which is taken to be a Lambertian surface. All scenarios are driven by the same downwelling irradiance, which is taken to be the direct solar beam. The scenarios are arranged in order of increasing actinic flux FJ ν , shown in the final column. FJ ν is seen to increase with both the number and the brightness of reflecting surfaces. The reflectivity of natural surfaces is no more than 90% nor less than 5%⁴ , so that the FJ ν in Table 3 represent bounds on realistic systems.

The Collimated-Non-Reflecting scenario assumes all light travels unidirectionally in a tightly collimated direct beam with irradiance Fν = F⊙ ν . In this case the radiation field is the delta function in the direction of the solar beam Iν() = F⊙ ν δ( -⊙). The closest approximation to this scenario in the natural environment is the pristine atmosphere high above the ocean in daylight. The actinic flux FJ ν = F⊙ ν follows directly from (31). We define the actinic flux enhancement S of a medium as the ratio of the actinic flux to the actinic flux of a collimated beam with the same incident irradiance

The Isotropic-Non-Reflecting scenario assumes isotropic downwelling radiance above a completely black surface. Under these conditions, the actinic flux is twice the incident irradiance because the photons are evenly distributed over the hemisphere rather than collimated. Moderately thick clouds (τ ~ > 3) over a dark surface such as the ocean create a radiance field approximately like this. However, because clouds are efficient at diffusing downwelling irradiance, they are also efficient reflectors and this significantly reduces the incident irradiance F- ν relative to the total extraterrestrial irradiance F⊙ ν . Whether photochemistry is enhanced or diminished beneath real clouds depends on whether the actinic flux enhancement factor FJ ν = 2 compensates the reduced sub-cloud insolation due to photons up-scattering off the cloud and back to space.

The Collimated-Reflecting scenario is a very important limit in nature because the net effect on photochemistry can approach the theoretical photochemical enhancement of S = 3. In this limit, collimated downwelling radiation and diffuse upwelling radiation combine to drive photochemistry from both hemispheres. A clear atmosphere above bright surfaces (clouds, desert, snow) approaches this limit. These conditions describe a large fraction of the atmosphere, which is 50–60% cloud-covered. Moreover, the incident flux is not attenuated by cloud transmission, so F- ν ≈ F⊙ ν . The relatively high frequency of occurance of this scenario, combined with the large photochemical enhancement of S = 3, are unequalled by any other scenario.

We may also define the actinic flux efficiency E of a medium as the actinic flux relative to the actinic flux of an isotropic radiation field with the same incident irradiance

The Isotropic-Reflecting scenario shows that an isotropic radiation field is most efficient for driving photochemistry. The maximum four-fold increase in efficiency relative to the collimated field arises from the radiation field interacting with the particle from all directions, rather than from one direction only. This is geometrically equivalent to multiplying the molecular cross section by a factor of four, the ratio between the surface and cross-sectional areas of a sphere. Note that photochemistry itself is driven by molecular absorption which reduces FJ ν from the values in Table 3.

In summary, we have learned that a reactant molecule with a spherically symmetric field of influence receives photochemical radiation much like Earth receives solar irradiance. In both cases the collimated beam intercepts one fourth of the total area of matter (molecule or Earth), while an equal flux of diffuse (or diurnal average) irradiance impinges on four times as much area. Since, in the geometric limit, absorption probability depends upon area, not direction, collimated beams have one fourth the photochemical potential as isotropic radiation.

2.1.6 Energy Density

Another quantity of interest is the density of radiant energy per unit volume of space. We call this quantity the energy density Uν. The energy density is the number of photons per unit volume in the frequency range [ν,ν + dν] times the energy per photon, hν. Uν is simply related to the actinic flux FJ ν (31) and thus to the mean intensity ν.

2.1.7 Spectral vs. Broadband

Until now we have considered only spectrally dependent quantities such as the spectral radiance Iν, spectral irradiance Fν, spectral actinic flux FJ ν , and spectral energy density Uν. These quantities are called spectral and are given a subscript of ν, λ, or because they are expressed per unit frequency, wavelength, or wavenumber, respectively. Each spectral radiant quantity may be integrated over a frequency range to obtain the corresponding band-integrated radiant quantity. Band-integrated radiant fields are often called narrowband or broadband Depending on the size of the frequency range, broadband radiant are obtained by integrating over all frequencies:

2.1.8 Thermodynamic Equilibria

Temperature plays a fundamental role in radiative transfer because T determines the population of excited atomic states, which in turn determines the potential for thermal emission. Thermal emission occurs as matter at any temperature above absolute zero undergoes quantum state transitions from higher energy to lower energy states. The difference in energy between the higher and lower level states is transferred via the electromagnetic field by photons. Thus an important problem in radiative transfer is quantifying the contribution to the radiation field from all emissive matter in a physical system. For the atmosphere the system of interest includes, e.g., clouds, aerosols, and the surface.

To develop this understanding we must discuss various forms of energetic equilibria in which a physical system may reside. Earth (and the other terrestrial planets, Mercury, Venus, and Mars) are said to be in planetary radiative equilibrium because, on an annual timescale the solar energy absorbed by the Earth system balances the thermal energy emitted to space by Earth. Radiation and matter inside a constant temperature enclosure are said to be in thermodynamic equilibrium, or TE.

Radiation in thermodynamic equilibrium with matter plays a fundamental role in radiation transfer. Such radiation is most commonly known as blackbody radiation. Kirchoff first deduced the properties of blackbody radiation.

Thermodynamic equilibrium (TE) is an idealized state, but, fortunately, the properties of radiation in TE can be shown to apply to a less restrictive equilibrium known as local thermodynamic equilibrium, or LTE.

2.1.9 Planck Function

The Planck function Bν describes the intensity of blackbody radiation as a function of temperature and wavelength

Bν(T,ν)	=	(41a)
Bλ(T,λ)	=	(41b)

The relations (41a) and (41b) predict slightly different quantities. The former predicts the blackbody radiance per unit frequency, while the latter predicts the blackbody radiance per unit wavelength. Of course these quantities are related since the blackbody energy within any given spectral band must be the same regardless of which formula is used to describe it. Expressed mathematically, this constraint means

The Planck function (41) has interesting behavior in both the high and the low energy photon limits. In the high energy limit, known as Wien’s limit, the photon energy greatly exceeds the ambient thermal energy

Bν(T,ν) =	e-hν∕kT	(46a)
Bλ(T,λ) =	e-hc∕λkT	(46b)

In the very low energy limit, known as the Rayleigh-Jeans limit, the photon energy is much less than the ambient thermal energy

Bν(T,ν)	≈	(49a)
Bλ(T,λ)	≈	(49b)

The frequency of extreme emission is obtained by taking the partial derivative of (41a) with respect to frequency with the temperature held constant

A separate relation may be derived for the wavelength of maximum emission λM by an analogous procedure starting from (41b). The result is

It is possible to use (51) to estimate the temperature of remotely sensed surfaces. For example, a satellite-borne tunable spectral radiometer may measure the emission of a newly discovered planet at all wavelengths of interest. Assuming the wavelength of peak measured emission is λM μm. Then a first approximation is that the planetary temperature is close to 2897.8∕λM.

2.1.10 Hemispheric Quantities

In climate studies we are most interested in the irradiance passing upwards or downwards through horizontal surfaces, e.g., the ground or certain layers in the atmosphere. These hemispheric irradiances measure the radiant energy transport in the vertical direction. These hemispheric irradiances depend only on the corresponding hemispheric intensities. Let us assume the intensity field is azimuthally independent, i.e., Iν = Iν(θ) only. Then the azimuthal contribution to (28) is 2π and

μ = | cos θ|	=	(54a)
μ = |u|	=	(54b)

(57)

(58)

The superscripts ⁺ and ^- denote upwelling (towards the upper hemisphere) and downwelling (towards the lower hemisphere) quantities, respectively. The net flux Fν is the difference between the upwelling and downwelling hemispheric fluxes, F + ν and F- ν , which are both positive definite quantities.

The hemispheric flux transport in isotropic radiation fields is worth examining in detail since this condition is often met in practice. It will be seen that isotropy considerably simplifies many of the troublesome integrals encountered. When Iν has no directional dependence (i.e., it is a constant) then F + ν = F- ν (58). Thus the net radiative energy transport is zero in an isotropic radiation field, such as a cavity filled with blackbody radiation.

Let us compute the upward transport of radiation Bν(T) (41a) emitted by a perfect blackbody such as the ocean surface. Since Bν(T) is isotropic, the intensity may be factored out of the definition of the upwelling flux (56a) and we obtain

2.1.11 Stefan-Boltzmann Law

The frequency-integrated hemispheric irradiance emanating from a blackbody of great interest since it describes, e.g., the radiant power of most surfaces on Earth. Although we could integrate the Bν (41a) directly to obtain the broadband intensity B ≡∫ ₀^∞B ν(ν) dν, it is traditional to integrate Bν first over the hemisphere (59). By proceeding in this order, we shall obtain the total hemispheric blackbody irradiance F + B in terms fundamental physical constants and the temperature of the body.

The definite integral in (61) is π⁴∕15. Proving this is a classic problem in mathematical physics which involves the Riemann zeta function (and thus prime number theory), the Gamma function, and contour integration. The procedure used to obtain this result is interesting so we briefly summarize it here. The Riemann zeta function ζ(x) for real x > 1 may be defined as

The integral (62) is analytically solvable for the special case of integers x = n. We may transform the integrand from a rational fraction into a power series using algebraic manipulation:

Recall that the sum of an infinite power series with initial term a₀ and ratio r is

Using the inifinite series representation (67) for the integrand of (62) yields

Contour integration in the complex plane gives analytic closed-form solutions to Σ₁^∞k^-n (70) for positive, even integers n. Our immediate concern is n = 4. Consider the complex function [Carrier et al., 1983, p. 97]

1. Is analytic throughout the complex plane
2. Has first order poles at all integer values on the real axis (except the origin)
3. Has a fifth order pole at the origin
4. Satisfies lim _|R|→∞f(z) = 0 where z = Re^iθ

Therefore (71) obeys the residue theorem for suitably chosen contours. In other words, fxm.

Having shown

Using Γ(4) = 3! = 3 × 2 × 1 = 6 and ζ(4) = π⁴∕90 from (72), Equation (63) becomes

We derived F + B (73) directly so that the Stefan-Boltzmann constant would fall naturally from the derivation. For completeness we now present the broadband blackbody intensity B

2.1.12 Luminosity

The total thermal emission of a body (e.g., star or planet) is called its luminosity, L. The luminosity is thermal irradiance integrated over the surface-area of a volume containing the body. The luminosity of bodies with atmospheres is usually taken to be the total outgoing thermal emission at the top of the atmosphere. If the time-mean thermal radiation field of a body is spherically symmetric, then its luminosity is easily obtained from a time-mean measurement of the thermal irradiance normal to any unit area of the surrounding surface. This technique is used on Earth to determine the intrinsic luminosity of stars whose distance D⊙ is known (e.g., through parallax) including our own.

An independent means of estimating the luminosity of a celestial body is to integrate the surface thermal irradiance over the surface area of the body. For planets without atmospheres, L is simply the integrated surface emission. Assuming the effective temperature (4) of our Sun is TE, the radius r⊙ at which this emission must originate is defined by combining (73) with (76)

2.1.13 Extinction and Emission

Radiation and matter have only two forms of interactions, extinction and emission. Lambert⁵ , first proposed that the extinction (i.e., reduction) of radiation traversing an infinitesimal path ds is linearly proportional to the incident radiation and the amount of interacting matter along the path

(78)

Here k(ν) is the extinction coefficient, a measurable property of the medium, and s is the absorber path length. k(ν) is proportional to the local density of the medium and is positive definite.

The term extinction coefficient and the exact definition of k(ν) are somewhat ambiguous until their physical dimensions are specified. Path length, for example, can be measured in terms of column mass path M [kg m-2], number path of molecules N [# m-2], and geometric distance s [m]. Each path measure has a commensurate extinction coefficient: the mass extinction coefficient ψe [m2 kg-1] (i.e., optical cross-section per unit mass), the number extinction coefficient κe [m2 molecule-1] (i.e., optical cross-section per molecule), the volume extinction coefficient ke [m2 m-3]=[m-1] (i.e., optical cross-section per unit concentration). These coefficients are inter-related by

We will develop the formalism of radiative transfer in this chapter in terms using the geometric path and volume extinction coefficient (s, ke) formalism. Our intent is to be concrete, rather than leaving the choice of units unstated. However, there is nothing fundamental about (s, ke). In Section 4.1.1 we state our preference for working in mass units (M, ψe).

Extinction includes all processes which reduce the radiant intensity. As will be described below, these processes include absorption and scattering, both of which remove photons from the beam. Similarly the radiative emission is also proportional to the amount of matter along the path

(80)

where Sν is known as the source function. The source function plays an important role in radiative transfer theory. We show in §2.2.1 that if Sν is known, then the full radiance field Iν is determined by an integration of Sν with the appropriate boundary conditions. Emission includes all processes which increase the radiant intensity. As will be described below, these processes include thermal emission and scattering which adds photons to the beam. Determination of k(ν), which contains all the information about the electromagnetic properties of the media, is the subject of active theoretical, laboratory and field research.

Extinction and emission are linear processes, and thus additive. Since they are the only two processes which alter the intensity of radiation,

2.1.14 Optical Depth

We define the optical path between points P₁ and P₂ as

The most frequently used form of optical path is the optical depth. The optical depth τ is the vertical component of the optical path , i.e., τ measures extinction between vertical levels. For historical reasons, the optical depth in planetary atmospheres is defined τ = 0 at the top of the atmosphere and τ = τ* at the surface. This convention reflects the astrophysical origin of much of radiative transfer theory. Much like pressure, τ is a positive definite coordinate which increases monotonically from zero at the top of the atmosphere to its surface value. Consider the optical depth between two levels z₂ > z₁, and then allow z₂ →∞

2.1.15 Geometric Derivation of Optical Depth

The optical depth of a column containing spherical particles may be derived by appealing to intuitive geometric arguments. Consider a concentration of N m-3 identical spherical particles of radius r residing in a rectangular chamber measuring one meter in the x and y dimensions and of arbitrary height. If the chamber is uniformly illuminated by a collimated beam of sunlight from one side, how much energy reaches the opposite side?

For this thought experiment, we will neglect the effects of scattering⁶ . Moreover, we will assume that the particles are partially opaque so that the incident radiation which they do not absorb is transmitted without any change in direction. Finally, assume the particles are homogeneously distributed in the horizontal so that the radiative flux F(z) is a function only of height in the chamber. Let us denote the power per unit area of the incident collimated beam of sunlight as F⊙. Our goal is to compute F(z) as a function of N and of r.

Since each particle absorbs incident energy in proportion to its geometric cross-section, the total absorption of radiation per particle is proportional to πr²Q e, where Qe = 1 for perfectly absorbing particles. If Qe < 1, then each particle removes fewer photons than suggested by its geometric size⁷ . Conversely, if Qe > 1 each particle removes more photons than suggested by its geometric size. Each particle encountered removes πr²Q eF(z) W m-2 from the incident beam. The maximum flux which can be removed from the beam is, of course, F⊙.

In a section of height Δz, the collimated beam passing through the chamber will encounter a total of NΔz particles. The number of particles encounted times the flux removed per particle gives the change in the radiative flux of the beam between the entrance and the rear wall

(88)

Finally we define the optical depth τ in terms of the extinction

We are now prepared to solve (87) for F(τ). From the theory of first order differential equations, we know that F must be an exponential function whose solution decays from its initial value with an e-folding constant of τ

(91)

Thus, in the limit of geometrical optics, the optical depth measures the number of e-foldings undergone by the radiative flux of a collimated beam passing through a given medium. This result is one form of the extinction law.

The exact value of Qe(r,λ) depends on the composition of the aerosol. However, there is a limiting value of Qe as particles become large compared to the wavelength of light.

(92)

Thus particles larger than 5–10 μm extinguish twice as much visible light as their geometric cross-section suggests.

To gain more insight into the usefulness of the optical depth, we can express τ (90) in terms of the aerosol mass M, rather than number concentration N. For a monodisperse aerosol of density ρ, the mass concentration is

(93)

If we substitute

into (90) we obtain

(94)

Typical cloud particles have r ~ 10 μm so that for visible solar radiation with λ ~ 0.5 μm we may employ (92) to obtain

(95)

Thus τ increases linearly with M for a given r. Note however, that a given mass M produces an optical depth that is inversely proportional to the radius of the particles!

2.1.16 Stratified Atmosphere

We obtain the radiative transfer equation in terms of optical path by substituting (82) into (81)

d	= u-1 dτ	(97a)
dτ	= u d	(97b)

Using the definitions of the half-range intensities (25) in (98) we obtain

- u	= I- ν - S- ν 0 < θ < π∕2,u = cos θ > 0	(99a)
- u	= I+ ν - S+ ν π∕2 < θ < π,u = cos θ < 0	(99b)

We now change variables from u to μ (54). Replacing u by μ in (99a) is allowed since μ = u > 0 in this hemisphere. In the upwelling hemisphere where π∕2 < θ < π, u < 0 so that μ = -u (54). This negative sign cancels the negative sign on the LHS of (99b), resulting in

- μ	= I- ν - S- ν	(100a)
μ	= I+ ν - S+ ν	(100b)

Equation (99) states that I± depends explicitly only on I± and on S±, but has no explicit dependence on I∓ or on S∓. Thus it may appear that I+ and I- are completely decoupled from eachother. However, we shall see that in problems involving scattering, S± depends explicitly on I∓ because scattering may change the trajectory of photons from upwelling to downwelling and visa versa. By coupling I+ to I-, scattering allows the entire radiance field to affect the radiance field at every point and in every direction (modulo the speed of light, of course). Thus scattering changes the solutions to (100) from being locally-dependent to depending on the global radiation field.

As a special case of (98), consider a stratified, non-scattering atmosphere in thermodynamic equilibrium. Then the source function equals the Planck function Sν = Bν = Bν(T) and we have

(101)

Equation (101) is the basis of radiative transfer in the thermal infrared, where scattering effects are often negligible. The solution to (101) is described in §2.2.2.

2.2 Integral Equations

2.2.1 Formal Solutions

It is useful to write down the formal solution to (98) before making additional assumption about the form of the source function Sν.

The solution for upwelling radiance is obtained by replacing u in (103) by -μ because u < 0 for upwelling intensities.

(105)

The first term on the RHS is the contribution of the boundary (e.g., Earth’s surface) to the upwelling intensity at level τ. This contribution is attenuated by the optical path of the radiation between the ground and level τ. The second term on the RHS is the contribution of the atmosphere to the upwelling intensity at level τ. The net upward emission of each parcel of air between the surface and level τ is S+ ν (τ′), but this internally emitted radiation is attenuated along the slant path between τ′ and τ. The μ^-1 factor in front of the integral accounts for the slant path of the emitting mass in the atmosphere.

The solution for downwelling radiance is obtained by replacing u in (103) by μ because μ = u in the downwelling hemisphere. The resulting expression must be integrated from the upper boundary down to level τ, and a boundary condition applied at the top.

The upwelling and downwelling intensities in a stratified atmosphere are fully described by (105) and (106). Such formal solutions to the equation of radiative transfer are of great heuristic value but limited practical use until the source function is known. Note that we have assumed a source function and boundary conditions which are azimuthally independent, but that the derivation of (105) and (106) does not rely on this assumption. It is straightforward to relax this assumption and replace Iν(τ,μ) by Iν(τ,μ,ϕ) and Sν(τ,μ) by Sν(τ,μ,ϕ) in the above.

2.2.2 Thermal Radiation In A Stratified Atmosphere

Consider a purely absorbing, stratified atmosphere in thermodynamic equilibrium. Then the source function equals the Planck function Sν = Bν = Bν(T) (41) and the radiative transfer equation is given by (101). It is important to remember that Bν is the complete source function only because we are explicitly neglecting all scattering processes. Thus we need only define the boundary conditions in order to use (105) and (106) to fully specify Iν. We assume that the surface emits blackbody radiation into the upper hemisphere

(107)

For brevity we shall define B* ν = Bν[T(τ*)]. At the top of the atmosphere, we assume there is no downwelling thermal radiation.

(108)

The solutions for upwelling and downwelling intensities are then obtained by using Sν(τ′,μ) = B ν(τ′) (the Planck function is isotropic) and substituting (107) and (108) into (105) and (106), respectively

2.2.3 Angular Integration

Once the solutions for the hemispheric intensities are known, it is straightforward to obtain the hemispheric fluxes by performing the angular integrations (56a)–(56b).

(111)

Consider first the downwelling flux in a non-scattering, thermal, isotropic, stratified atmosphere obtained by substituting (110) into (111) and interchanging the order of integration

In terms of exponential integrals defined in Appendix 10.7, the inner integral in parentheses in (112) is E₂(τ - τ′) (603c).

(113)

Similar terms arise when we consider the horizontal upwelling flux obtained by substituting (109) into (56a) and we obtain

2.2.4 Thermal Irradiance

Assume a non-scattering planetary surface at temperature T emits blackbody radiation such that Iν(τ*, +μ) = B ν(T) (107). What is the total upwelling thermal irradiance from the surface? From (56a) we have

We integrate (116) over frequency to obtain the total upwelling thermal irradiance

2.2.5 Grey Atmosphere

Consider an atmosphere transparent to solar radiation and partially opaque to thermal radiation governed by (101). The exact solution, including angular dependence, is given in §2.2.2. The hemispheric fluxes and net flux may only be obtained exactly by accounting for the angular dependence as in §2.2.3. We may eliminate the angular dependence of the net flux by making the simplifying assumption that the hemispheric up and downwelling irradiances equal a constant times the corresponding intensity. A number of methods exist to determine this constant, called the diffusivity factor (e.g., §2.2.15). These methods are all related to the two-stream approximation.

One such method [Salby, 1996, p. 232] is to identify an effective inclination along which all radiation is assumed to travel. With this assumption, the contribution to upwelling irradiance from the lower boundary, the first term on the RHS of 115), is

Using this assumption (118), the upwelling hemispheric irradiance (56a) for blackbody radiation is

Based on (119) and (120), the irradiance structure of the thermal atmosphere may approximated by performing a direct angular integration of (100). With our approximation, radiances I integrate directly to irradiances F modulo the diffusivity factor . (101)

-	= F- ν - πFB ν	(121a)
	= F + ν - πFB ν	(121b)

It is instructive to examine an idealized grey atmosphere, where the fluxes of interest have no spectral dependence. Although this is far from true in Earth’s atmosphere, the solution is straightforward and sheds light on the radiative equilibrium temperature profile and the greenhouse effect. We simplify (121) in two ways. First, we introduce a scaled optical depth d = ^-1dτ = dτ. Second, we drop the frequency dependence, which is equivalent to integrating over a broad range of frequencies. For heuristic purposes, think of this integration as being over the relatively narrow λ = 5–20 μm range where most of Earth’s terrestrial radiative energy resides.

-	= F-- πB	(122a)
	= F+ - πB	(122b)

In the absence of dynamical, chemical, and latent heating, the energy deposition in a parcel of air is entirely radiative. Under these conditions the idealize grey atmosphere described by (122) will adjust to a temperature profile determined by radiative equilibrium. Let the time rate of change of temperature T of a parcel be denoted by h [K s-1], the parcel warming rate⁸ or cooling rate. The warming rate is the rate of net energy deposition divided by the specific heat capacity at constant pressure cp [J kg-1 K-1] times the density ρ [kg m-3]

By definition, the net radiative heating vanishes (dT∕dz = 0) at all levels of an atmosphere in radiative equilibrium. Setting (125) equal to zero and integrating, we obtain

Adding and subtracting (122), we obtain

(F+ - F-)	= F+ + F-- 2πB	(126a)
(F+ + F-)	= F+ - F- = F 0	(126b)

	= ψ - 2πB	(127a)
	= ϕ = F0	(127b)

The upwelling irradiance at the surface is F+( *) = πB(T s) where Ts [K] is the surface skin temperature. The atmospheric temperature just above the surface is given by (131) as B( *) = ( * + 1). Thus there is a temperature discontinuity between the near-surface air and the ground.

Climate models typically express h (124) in terms of the flux gradient with respect to pressure by invoking the hydrostatic equilibrium condition (456)

2.2.6 Scattering

Energy interacting with matter undergoes one of two processes, scattering or absorption. Scattering occurs when a photon reflects off matter without absorption. The direction of the photon after the interaction is usually not the same as the incoming direction. The case where the scattered photons are homogeneously distributed throughout all 4π steradians is called isotropic scattering. In general, the angular dependence of the scattering is described by the phase function of the interaction. The phase function p is closely related to the probability that photons incoming from the direction ′ = (θ′,ϕ′) will (if scattered) scatter into outgoing direction = (θ,ϕ). It is usually assumed that p depends only on the scattering angle Θ between incident and emergent directions.

The Cartesian components of Ω′ and Ω are straigtforward to obtain in spherical polar coordinates.

The scattering angle Θ is simply related to the inner product of ′ and by the cosine law

2.2.7 Phase Function

Accurate treatment of the angular scattering of radiation, i.e., the phase function, is, perhaps, makes rigorous demands of radiative transfer applications. A correspondingly large body of literature is devoted to this topic. Essential references include van de Hulst [1957], Joseph et al. [1976], Wiscombe and Grams [1976], Wiscombe [1977], Wiscombe [1979, edited/revised 1996], and Boucher [1998].

The phase function p(cos Θ) is normalized so that the total probability of scattering is unity

∫ Ωp(cos Θ) dΩ	= 1	(138a)
∫ ϕ=0ϕ=2π ∫ θ=0θ=πp(θ′,ϕ′; θ,ϕ) sin θ dθ dϕ	= 1	(138b)

Scattering may depend on the absolute directions ′ and themselves, rather than just their relative orientations as measured by the angle Θ between them. This might be the case, for example, in a broken sea-ice field. For the time being, however, we shall assume that the phase function depends only on Θ.

In atmospheric problems, the phase function may often be independent of the azimuthal angle ϕ, and depend only on θ. In this case the phase function normalization (138b) simplifies to

2.2.8 Legendre Basis Functions

The phase function (138) specifies the angular distribution of scattering. Solution of any radiative transfer equation involving scattering depends on it. Numerical approaches aim for a suitable approximation or discretization which represents p(Θ) to some desired level of accuracy. The optimal basis functions for representing p(Θ) are the Legendre polynomials.

The Legendre polynomial expansion of the phase function is

The Legendre polynomials are orthonormal on the interval [-1, 1]. The factor of 2n + 1 that appears in the numerator of the Legendre expansion (146) also appears in the denominator of the Legendre polynomial orthonormality property:

∫ u=-1u=1P n(u)Pk(u) du	= δn,k	(141a)
∫ u=-1u=1P n(cos Θ)Pk(cos Θ) d(cos Θ)	= δn,k	(141b)

The expansion coefficients χn are defined by the projection of the corresponding Legendre polynomial onto the phase function.

χn	= ∫ u=-1u=1P n(u)p(u) du	(142a)
	= ∫ Θ=πΘ=0P n(cos Θ)p(cos Θ) d(cos Θ)	(142b)
	= ∫ Θ=0Θ=πP n(cos Θ)p(cos Θ) sin Θ dΘ	(142c)

The Legendre expansion for simple, smoothly varying and symmetric phase functions is highly accurate with just few terms. For instance, a single moment Legendre expansion exactly describes Rayleigh scattering. However, the expansions of asymmetric phase functions converge much more slowly. The strongly peaked forward scattering lobes which appear at large size parameter may require hundreds of moments for accurate representation.

Computing the χn using quadrature methods significantly increases accuracy and reduces time. Wiscombe [1977] recommends Gauss-Lobatto quadrature for highly asymmetric phase functions:

2.2.9 Henyey-Greenstein Approximation

Computation of the exact phase function of scatterers is laborious but desirable when the objective is to predict directional radiances. However, phase function approximations often suffice when hemispheric fluxes are the objective. The Henyey-Greenstein Phase Function approximates the full phase function in terms of its first Legendre moment, i.e., its asymmetry parameter g.

2.2.10 Direct and Diffuse Components

When working with half-range intensities it is convenient to decompose the downwelling intensity I- into the sum of a direct component, I- s (τ,), and a diffuse component, I- d (τ,) such that

It is plain that I- s (τ,) is zero in all directions except that of the collimated beam. Moreover, the intensity in the direct beam is, by definition, subject only to extinction by Bougher’s law (i.e., there is no emission). Thus the direct beam component of the downwelling intensity is the irradiance incident at the top of the atmosphere, attenuated by the extinction law (91),

2.2.11 Source Function

The source function for thermal emission St ν is

All terms in (152) are functions of position and scattering process. If there are n distinct scattering processes (e.g., Rayleigh scattering, aerosol scattering, cloud scattering) then we must know the properties of each individual process (e.g., σ_i, p_i) and sum their contributions to obtain to total scattering source function Ss ν = ∑ _i=1ⁱ⁼ⁿS _ν,i^s. We neglect such details for now and merely remind the reader that applications need to consider multiple extinction processes at the same time.

The total source function is obtained by adding the source functions for thermal emission (150) and for scattering (152)

2.2.12 Radiative Transfer Equation in Slab Geometry

Inserting (153) into (96) we obtain the radiative transfer equation including absorption, thermal emission and scattering

We are most interested in solving the radiative transfer equations in a slab geometry. To do this we transform from the path optical depth coordinate to the vertical optical depth coordinate τ. Applying the procedure described in (96)–(100) to (154) we obtain the radiative transfer equation for the half range intensities in a slab geometry

- μ	= I-(τ,) - (1 - ϖ)B - ∫ +I+(τ,′)p(+′,-) dΩ′
	- ∫ -I-(τ,′)p(-′,-) dΩ′	(155a)
μ	= I+(τ,) - (1 - ϖ)B - ∫ +I+(τ,′)p(+′, +) dΩ′
	- ∫ -I-(τ,′)p(-′, +) dΩ′	(155b)

Likewise, inserting (147) and (150) into (155b), we obtain the equation for the diffuse upwelling intensity

2.2.13 Azimuthal Mean Radiation Field

Often we are not concerned with the azimuthal dependence of the radiation field. In some cases this is because the azimuthal dependence is very weak. For example, heavily overcast skies, or diffuse reflectance from a uniform surface. In other cases the azimuthal dependence may not be weak, but there insufficient information to fully determine the solutions. For example, the full surface BRDF or the shapes of clouds are not available. The azimuthal dependence of a radiative quantity X(ϕ) is removed by applying the azimuthal mean operator

Applying (160) to (25) and (138), we obtain the azimuthal mean hemispheric intensity, phase function, and single-scattering source function respectively,

Applying (160) to the full radiative transfer equations for slab geometry, (156) and (158), we obtain the radiative transfer equations for the azimuthal mean intensity in slab geometry

-μ = I- d (τ,μ) - (1 - ϖ)B -F⊙e-τ∕μ0 p(-μ0,-μ)
	- ∫ +I+ d (τ,μ′)p(+μ′,-μ) dμ′- ∫ -I- d (τ,μ′)p(-μ′,-μ) dμ′	(165a)
μ = I+ d (τ,μ) - (1 - ϖ)B -F⊙e-τ∕μ0 p(-μ0, +μ)
	- ∫ +I+ d (τ,μ′)p(+μ′, +μ) dμ′- ∫ -I- d (τ,μ′)p(-μ′, +μ) dμ′	(165b)

We apply the two stream formalism introduced in §2.4 to the radiative transfer equation for the azimuthal mean radiation field (165). Recall that the two stream formalism is obtained by applying three successive approximations. The approximation is achieved in three stages. First, we operate on both sides of (165a) with the hemispheric averaging operator (161). The procedure for (165a) is identical.

Further approximations are required to simplify equations (173). These approximations allow us to decouple products of functions with μ dependence. The first approximation, (167), replaces the continuous value μ on the LHS of (165) by a suitable hemispheric mean value .

The first two terms on the RHS of (166) are simple hemispheric integrals of hemispherically-varying intensities. We will replace I±(τ,μ) by the hemispheric mean intensity I±(τ) (205). The Planck function is isotropic and so may be pulled outside the hemispheric integral which promptly vanishes since ∫ ₀¹ dμ = 1. Thus the hemispherically integrated intensities which result explicitly are replaced by the hemispheric mean intensity (205) associated with .

The final three terms on the RHS of (166) involve products of the phase function and the intensity. We approximate the scattering integrals by applying the hemispheric integral over μ to the intensities, and then extracting the intensities from original integrals over μ′.

The hemispheric mean backscattering coefficient b is the hemispheric integral of (169)

The result of these two steps is

- = I- d (τ) - (1 - ϖ)B -F⊙e-τ∕μ0 p(-μ0,-μ)
	- ∫ +I+ d (τ)p(+μ′,-μ) dμ′- ∫ -I- d (τ)p(-μ′,-μ) dμ′	(173a)
= I+ d (τ) - (1 - ϖ)B -F⊙e-τ∕μ0 p(-μ0, +μ)
	- ∫ +I+ d (τ)p(+μ′, +μ) dμ′- ∫ -I- d (τ)p(-μ′, +μ) dμ′	(173b)

2.2.14 Anisotropic Scattering

Armed with techniques introduced in solving the two stream equations for an isotropically scattering medium §2.4.1, we now solve analytically the coupled two stream equations for an anisotropic medium. First we recast (??) into a simpler notation

	= I+ ν (τ) - ϖ(1 - b)I+ ν (τ) - ϖbI- ν (τ) - (1 - ϖ)Bν - S*+ ν	(174a)
-	= I- ν (τ) - ϖ(1 - b)I- ν (τ) - ϖbI+ ν (τ) - (1 - ϖ)Bν - S*- ν	(174b)

	= -(α - β)(I+ - I-)	(175a)
	= -(α + β)(I+ + I-)	(175b)

α	≡-[1 - ϖ(1 - b)]∕	(176a)
β	≡ ϖb∕	(176b)

2.2.15 Diffusivity Approximation

The angular integration required to convert the intensity field into a hemispheric irradiance is a time consuming aspect of numerical models and should be avoided where possible. According to (56a)–(56b)

F + ν	= 2 ∫ 01I ν(+μ)μ dμ	(177a)
F- ν	= 2 ∫ 01I ν(-μ)μ dμ	(177b)

Two point full Gaussian quadrature (§10.5) tells us the optimal angles to evaluate (177) at are u = ±3^-1∕2, θ = ±54.7356∘.

F + ν	≃Iν(+3-1∕2)	(178a)
F- ν	≃Iν(-3-1∕2)	(178b)

What is often used to evaluate (177) instead of Gaussian quadrature is a diffusivity approximation. Instead of choosing optimal quadrature angles μ_k, the diffusivity approximation relies on choosing the optimal effective absorber path. Thus the diffusivity factor, D, is defined by

F + ν	≃ Iν(u = +D-1)	(180a)
F- ν	≃ Iν(u = -D-1)	(180b)

The reasons for subsuming so many factors into D are historical. Nevertheless, (180) is much more common than (178) in the literature. Heating rate calculations in isotropic, non-scattering atmospheres show that using D = 1.66 results in errors < 2% [Goody and Yung, 1989, p. 221].

In applying (180), D only affects the optical depth factor of the intensity. Thus computing the intensity in (180) at the angle θ = θD is formally equivalent to computing the intensity in the vertical direction θ = 0 but through an atmosphere in which the absorber densities have been increased by a factor of D.

2.2.16 Transmittance

The transmittance between two points measures the transparency of the atmosphere between those points. The transmittance from point P₁ to point P₂ is the likelihood a photon traveling in direction at P₁ will arrive at P₂ without having interacted with the matter in between.

(181)

We have left implicit the dependence on frequency. In a stratified atmosphere, we are interested in the transmission of light traveling from layer z′ to layer z at angle θ from the vertical, thus

(182)

Thus, in common with absorptance and reflectance R, the transmittance T [0, 1]⁹

The vertical gradient of T is frequently used to formulate solutions of the radiative transfer equation in non-scattering, thermal atmospheres.

(184)

The integral solutions for the upwelling and downwelling intensities in a non-scattering, thermal, stratified atmosphere (109)–(110) may be rewritten in terms of T . (110) from τ to z. To do this we use

Iν(z, +μ)	= T(0,z; μ)Bν(0) + ∫ 0zBν(z′) dz′	(187a)
Iν(z,-μ)	= -∫ z∞Bν(z′) dz′	(187b)

2.3 Reflection, Transmission, Absorption

Perhaps the most useful metrics of the radiative properties of entire systems are the quantities reflection, transmission, and absorption. These metrics describe normalized properties of a system and thus are somewhat more fundamental than absolute quantities (like transmitted irradiance), which may, for example, change depending on time of day. We shall define the reflection, transmission, and absorption first of the radiance field, and then integrate these definitions with the appropriate weighting to obtain the R, T , and pertinent to fluxes. We shall introduce a painful but necessary menagerie of terminology to describe the various species of R, T , and .

Molotch et al. [2004] show how remotely sensed surface reflectance differs significantly from that derived from snow-age models. Using the more realistic albedos improved models of snowmelt by providing more accurate net surface shortwave radiation estimates. This term, FsfcSW, dominates the snow-melt energy budget.

Wang et al. [2005] derived a simple functional fit to the desert surface albedo observed by MODIS.

2.3.1 BRDF

The bidirectional reflectance distribution function (BRDF) ρ is the ratio of the reflected intensity to the energy in the incident beam. As such, ρ is a function of frequency, incident angle, and scattered angle

The dimensions of ρ are sr-1 and they convert irradiance to intensity. The reflected intensity in any particular direction is the sum of contributions from all incident directions ′ that have a finite probability of reflecting into .

2.3.2 Lambertian Surfaces

Fortunately many surfaces found in nature obey simpler reflectance properties first characterized by Lambert. A Lambertian surface is one whose reflectance is independent of both incident and reflected directions. The reflectance of a Lambertian surface depends only on frequency

The irradiance reflected from a Lambertian surface F + νr is the cosine-weighted integral of the reflected intensity (191) over all reflected angles

The reflected irradiance may not exceed the incident irradiance or the requirement of energy conservation will be violated. Therefore (192) shows that

Many researchers prefer the Lambertian BRDF to have an upper limit of 1, not π^-1 (193). Thus it is common to encounter in the literature BRDF defined as r = πρ.

2.3.3 Albedo

The reflectance of a surface illuminated by a collimated source such as the Sun or a laser is of great interest in planetary studies and in remote sensing. For an incident collimated beam of intensity F⊙, the diffusely reflected intensity is

When the reflecting body is of finite size then the corresponding ratio of reflected to incident fluxes is more complicated because it contains “edge effects”, e.g., diminishing contributions from planetary limbs. First we consider the reflection, called the spherical albedo, planetary albedo or Bond albedo, of an entire planetary disk illuminated by collimated sunlight.

The directional reflectance ρ(ν,-2π,) is the

The flux reflectance of a Lambertian surface to collimated light is found by subsituting (190) in (196)

The BRDF must be integrated over all possible angles of reflectance in order to obtain the flux reflectance, ρ(ν,-′ , 2π). Representative flux reflectances of various surfaces in the Earth system are presented in Table 4.

2.3.4 Flux Transmission

The flux transmission TF ν between two layers is (twice) the cosine-weighted integral of the spectral transmission

(200)

The factor of 2 ensures that expressions for the hemispheric fluxes in non-scattering, thermal atmospheres closely resemble the equations for intensities, as will be shown next.

The vertical gradient of the flux transmission is frequently used to formulate solutions of the radiative transfer equation in non-scattering, thermal atmospheres. It is obtained by differentiating (200) then applying (602a)

The integral solutions for the upwelling and downwelling fluxes in a non-scattering, thermal, stratified atmosphere (113) and (114) may be rewritten in terms of TF ν . The procedure for doing so is analogous to the procedure applied to intensities in §2.2.16. First, we change the independent variable in (113) and (114) from τ to z.

F + ν (z)	= Bν(0)TF ν (0,z) + ∫ z′=0z′=zB ν(z′) dz′	(202a)
F- ν (z)	= -∫ z′=zz′=∞B ν(z′) dz′	(202b)

Equations (202a)–(202b) are particularly useful because they contain no explicit reference to the angular integration. The polar integration, however, is still implicit in the definition of TF ν (199). The hemispherical irradiances may be determined without any angular integration if the diffusivity approximation is applied to (202a)–(202b). This is accomplished by replacing TF ν by the vertical spectral transmission T through a factor D (180) times as much mass

(203)

Determination of T requires the use of a spectral gaseous extinction database in conjunction with either a line-by-line model, a narrow band model, or a broadband emissivity approach.

2.4 Two-Stream Approximation

The two-stream approximation to the radiative transfer equation assumes that up- and down-welling radiances travel at mean inclinations + and -, respectively. With this assumption, the hemispheric intensities Iν^± in the azimuthally averaged radiative transfer equation (165) lose their explicit dependence on μ and depend solely on optical depth. For an isotropically scattering atmosphere in a slab geometry, the two-stream approximation may be written

+	= I+ ν (τ) -I+ ν (τ) -I- ν (τ) - (1 - ϖ)Bν[T(τ)]	(204a)
--	= I- ν (τ) -I+ ν (τ) -I- ν (τ) - (1 - ϖ)Bν[T(τ)]	(204b)

We now re-define the other important properties of the radiation field to be consistent with the two-stream approximation. Perhaps most important is the hemispheric irradiance. Starting from (56) we obtain

A common mis-understanding of the two-stream approximation is that the radiance field is isotropic in each hemisphere. If two-stream radiance were isotropic in each hemisphere then the appropriate measure of hemispheric irradiance would be F± ν = πI± ν (59) as opposed to F± ν = 2π±I± ν (206). In fact, assuming the intensity field varies linearly with μ (and is thus anisotropic) is a reasonable interpretation of the two-stream approximation.

The mean intensity ν (20) for an azimuthally independent, isotropically scattering radiation field is

Continuing to write the two-stream approximations for radiative quantities of interest, we turn now to the source function. The source function Sν (153) for an azimuthally independent, isotropically scattering (p = 1) radiation field is

A comment on vertical homogeneity is appropriate here. The Planck function in (208) varies continuously with τ, whereas ϖ is assumed to be constant in a given layer over which (208) is discretized. Thus τ-dependence appears explicitly for Bν, and is absent for ϖ. More generally, the radiation field (intensity, irradiance, source function, etc.) varies continuously with τ. The discretization of τ (or, in general, all spatial coordinates) is required in order to solve for the continuously varying radiation field. The discretization determines the scale over which the properties of the matter in the medium may change.

The radiative heating rate hν in the two-stream approximation may be obtained in at least two distinct forms.

Together the relations for Fν, ν, Sν, and hν in (206)–(209) define a complete and self-consistent set of radiative properties under the two-stream approximation. We turn now to obtaining the two-stream intensities I± ν (τ) from which these quantities may be computed.

2.4.1 Two-Stream Equations

We shall simplify both the physics and nomenclature of (204) before attempting an analytical solution. First, for algebraic simplicity we set + = - = . Also, we neglect the thermal source term B ν[T(τ)]. With these assumptions, (204) becomes

	= I+ ν (τ) -I+ ν (τ) -I- ν (τ)	(210a)
-	= I- ν (τ) -I+ ν (τ) -I- ν (τ)	(210b)

For the remainder of the derivation we drop the ν subscript and the explicit dependence of I^± on τ. Adding and subtracting the equations in (210) we obtain

	= I+ + I-- ϖ(I+ + I-) = (1 - ϖ)(I+ + I-)	(211a)
	= I+ - I-	(211b)

Y ±	≡ I+ ± I-	(212a)
I±	= (Y + ± Y -)	(212b)

	= Y +	(213a)
	= Y -	(213b)

	=		= Y -		= Y -
	=		= Y +		= Y +