Online: http://dust.ess.uci.edu/facts Updated: Thu 24th Jan, 2008, 14:14 Radiative Transfer in the Earth System by Charlie Zender University of California, Irvine Department of Earth System Science University of California Irvine, CA 92697-3100 zender@uci.edu Voice: (949) 824-2987 Fax: (949) 824-3256 Copyright c 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.ai.mit.edu/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 Table 1: Freely Available Community Texts Format URL Location Radiative Transfer in the Earth System DVI http://dust.ess.uci.edu/facts/rt/rt.dvi PDF http://dust.ess.uci.edu/facts/rt/rt.pdf Postscript http://dust.ess.uci.edu/facts/rt/rt.ps Particle Size Distributions: Theory and Application to Aerosols, Clouds, and Soils DVI http://dust.ess.uci.edu/facts/psd/psd.dvi PDF http://dust.ess.uci.edu/facts/psd/psd.pdf Postscript http://dust.ess.uci.edu/facts/psd/psd.ps Natural Aerosols in the Climate System DVI http://dust.ess.uci.edu/facts/aer/aer.dvi PDF http://dust.ess.uci.edu/facts/aer/aer.pdf Postscript http://dust.ess.uci.edu/facts/aer/aer.ps 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. ii Notes for Students of ESS 200B, Earth System Physics: 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 166 is also helpful. Notes for Students of ESS 236, Radiative Transfer and Remote Sensing: Yada yada yada. CONTENTS iii Contents Contents List of Figures List of Tables 1 Introduction 1.1 Planetary Radiative Equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Radiative Transfer Equation 2.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Intensity . . . . . . . . . . . . . . . . . . . . 2.1.2 Mean Intensity . . . . . . . . . . . . . . . . . 2.1.3 Irradiance . . . . . . . . . . . . . . . . . . . . 2.1.4 Actinic Flux . . . . . . . . . . . . . . . . . . 2.1.5 Actinic Flux Enhancement . . . . . . . . . . . 2.1.6 Energy Density . . . . . . . . . . . . . . . . . 2.1.7 Spectral vs. Broadband . . . . . . . . . . . . . 2.1.8 Thermodynamic Equilibria . . . . . . . . . . . 2.1.9 Planck Function . . . . . . . . . . . . . . . . 2.1.10 Hemispheric Quantities . . . . . . . . . . . . . 2.1.11 Stefan-Boltzmann Law . . . . . . . . . . . . . 2.1.12 Luminosity . . . . . . . . . . . . . . . . . . . 2.1.13 Extinction and Emission . . . . . . . . . . . . 2.1.14 Optical Depth . . . . . . . . . . . . . . . . . . 2.1.15 Geometric Derivation of Optical Depth . . . . 2.1.16 Stratified Atmosphere . . . . . . . . . . . . . 2.2 Integral Equations . . . . . . . . . . . . . . . . . . . . 2.2.1 Formal Solutions . . . . . . . . . . . . . . . . 2.2.2 Thermal Radiation In A Stratified Atmosphere 2.2.3 Angular Integration . . . . . . . . . . . . . . . 2.2.4 Thermal Irradiance . . . . . . . . . . . . . . . 2.2.5 Grey Atmosphere . . . . . . . . . . . . . . . . 2.2.6 Scattering . . . . . . . . . . . . . . . . . . . . 2.2.7 Phase Function . . . . . . . . . . . . . . . . . 2.2.8 Legendre Basis Functions . . . . . . . . . . . 2.2.9 Henyey-Greenstein Approximation . . . . . . 2.2.10 Direct and Diffuse Components . . . . . . . . 2.2.11 Source Function . . . . . . . . . . . . . . . . 2.2.12 Radiative Transfer Equation in Slab Geometry 2.2.13 Azimuthal Mean Radiation Field . . . . . . . . 2.2.14 Anisotropic Scattering . . . . . . . . . . . . . iii vi 1 1 1 2 4 4 4 5 7 8 13 15 15 15 16 18 20 24 24 25 26 28 30 30 31 32 33 34 37 38 39 40 40 41 42 44 46 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CONTENTS 2.2.15 Diffusivity Approximation . . 2.2.16 Transmittance . . . . . . . . . Reflection, Transmission, Absorption 2.3.1 BRDF . . . . . . . . . . . . . 2.3.2 Lambertian Surfaces . . . . . 2.3.3 Albedo . . . . . . . . . . . . 2.3.4 Flux Transmission . . . . . . Two-Stream Approximation . . . . . 2.4.1 Two-Stream Equations . . . . 2.4.2 Layer Optical Properties . . . 2.4.3 Conservative Scattering Limit Solar Heating . . . . . . . . . . . . . Chapter Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv 47 48 49 50 50 51 52 54 56 58 59 61 61 62 62 63 66 66 67 68 69 69 69 69 70 70 70 72 74 79 80 81 82 84 86 87 87 90 90 92 92 95 2.3 2.4 2.5 2.6 3 Remote Sensing 3.1 Rayleigh Limit . . . . . . . . . . . . . . . . . . 3.2 Anomalous Diffraction Theory . . . . . . . . . . 3.3 Geometric Optics Approximation . . . . . . . . . 3.4 Single Scattered Intensity . . . . . . . . . . . . . 3.5 Satellite Orbits . . . . . . . . . . . . . . . . . . 3.6 Aerosol Characterization . . . . . . . . . . . . . 3.6.1 Measuring Aerosol Optical Depth . . . . 3.6.2 Aerosol Indirect Effects on Climate . . . 3.6.3 Aerosol Effects on Snow and Ice Albedo ˚ 3.6.4 Angstrøm Exponent . . . . . . . . . . . Gaseous Absorption 4.1 Line Shape . . . . . . . . . 4.1.1 Line Shape Factor . 4.1.2 Natural Line Shape . 4.1.3 Pressure Broadening 4.1.4 Doppler Broadening 4.1.5 Voigt Line Shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Molecular Absorption 5.1 Mechanical Analogues . . . . 5.1.1 Vibrational Transitions 5.1.2 Isotopic Lines . . . . . 5.1.3 Combination Bands . . 5.2 Partition Functions . . . . . . 5.3 Dipole Radiation . . . . . . . 5.4 Two Level Atom . . . . . . . 5.5 Line Strengths . . . . . . . . . 5.5.1 HITRAN . . . . . . . . 5.6 Line-By-Line Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CONTENTS 5.6.1 6 v Literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 95 96 96 97 97 98 101 101 103 105 106 115 115 116 116 117 119 120 120 121 121 122 125 126 127 128 130 131 131 132 133 133 134 134 134 134 134 135 135 136 136 136 Band Models 6.1 Generic . . . . . . . . . . . . . . . . . . . . . 6.1.1 Beam Transmittance . . . . . . . . . . 6.1.2 Beam Absorptance . . . . . . . . . . . 6.1.3 Equivalent Width . . . . . . . . . . . . 6.1.4 Mean Absorptance . . . . . . . . . . . 6.2 Line Distributions . . . . . . . . . . . . . . . . 6.2.1 Line Strength Distributions . . . . . . . 6.2.2 Normalization . . . . . . . . . . . . . 6.2.3 Mean Line Intensity . . . . . . . . . . 6.2.4 Mean Absorptance of Line Distribution 6.2.5 Transmittance . . . . . . . . . . . . . . 6.2.6 Multiplication Property . . . . . . . . . 6.3 Transmission in Inhomogeneous Atmospheres . 6.3.1 Constant mixing ratio . . . . . . . . . . 6.3.2 H-C-G Approximation . . . . . . . . . 6.4 Temperature Dependence . . . . . . . . . . . . 6.5 Transmission in Spherical Atmospheres . . . . 6.5.1 Chapman Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Radiative Effects of Aerosols and Clouds 7.1 Single Scattering Properties . . . . . . . . . . . . 7.1.1 Maxwell Equations . . . . . . . . . . . . 7.2 Separation of Variables . . . . . . . . . . . . . . 7.2.1 Azimuthal Solutions . . . . . . . . . . . 7.2.2 Polar Solutions . . . . . . . . . . . . . . 7.2.3 Radial Solutions . . . . . . . . . . . . . 7.2.4 Plane Wave Expansion . . . . . . . . . . 7.2.5 Boundary Conditions . . . . . . . . . . . 7.2.6 Mie Theory . . . . . . . . . . . . . . . . 7.2.7 Resonances . . . . . . . . . . . . . . . . 7.2.8 Optical Efficiencies . . . . . . . . . . . . 7.2.9 Optical Cross Sections . . . . . . . . . . 7.2.10 Optical Depths . . . . . . . . . . . . . . 7.2.11 Single Scattering Albedo . . . . . . . . . 7.2.12 Asymmetry Parameter . . . . . . . . . . 7.2.13 Mass Absorption Coefficient . . . . . . . 7.3 Effective Single Scattering Properties . . . . . . 7.3.1 Effective Efficiencies . . . . . . . . . . . 7.3.2 Effective Cross Sections . . . . . . . . . 7.3.3 Effective Specific Extinction Coefficients 7.3.4 Effective Optical Depths . . . . . . . . . 7.3.5 Effective Single Scattering Albedo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LIST OF FIGURES 7.4 7.3.6 Effective Asymmetry Parameter . . . . . . . . Mean Effective Single Scattering Properties . . . . . . 7.4.1 Mean Effective Efficiencies . . . . . . . . . . 7.4.2 Mean Effective Cross Sections . . . . . . . . . 7.4.3 Mean Effective Specific Extinction Coefficients 7.4.4 Mean Effective Optical Depths . . . . . . . . . 7.4.5 Mean Effective Single Scattering Albedo . . . 7.4.6 Mean Effective Asymmetry Parameter . . . . . Bulk Layer Single Scattering Properties . . . . . . . . 7.5.1 Addition of Optical Properties . . . . . . . . . 7.5.2 Bulk Optical Depths . . . . . . . . . . . . . . 7.5.3 Bulk Single Scattering Albedo . . . . . . . . . 7.5.4 Bulk Asymmetry Parameter . . . . . . . . . . 7.5.5 Diagnostics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi 136 136 137 138 138 138 138 139 139 139 139 139 140 140 140 140 150 150 152 152 153 154 154 154 155 155 158 159 160 166 7.5 8 9 Global Radiative Forcing Implementation in NCAR models 10 Appendix 10.1 Vector Identities . . . . . . . . . . . 10.2 Legendre Polynomials . . . . . . . 10.3 Spherical Harmonics . . . . . . . . 10.4 Bessel Functions . . . . . . . . . . 10.4.1 Spherical Bessel Functions . 10.4.2 Recurrence Relations . . . . 10.4.3 Power Series Representation 10.4.4 Asymptotic Values . . . . . 10.5 Gaussian Quadrature . . . . . . . . 10.6 Gauss-Lobatto Quadrature . . . . . 10.7 Exponential Integrals . . . . . . . . Bibliography Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Figures 1 2 3 4 5 6 7 Cross Section and Quantum Yield of Nitrogen Dioxide Vertical Distribution of Photodissociation Rates . . . . Climatological Mean Absorbed Solar Radiation . . . . Climatological Mean Emitted Longwave Radiation . . ENSO Temperature and OLR . . . . . . . . . . . . . . Seasonal Shortwave Cloud Forcing . . . . . . . . . . . Zonal Mean Shortwave Cloud Forcing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 12 141 142 143 144 145 8 9 10 11 Seasonal Longwave Cloud Forcing . . . Zonal Mean Longwave Cloud Forcing . Climatological Mean Net Cloud Forcing ENSO Cloud Forcing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 147 148 149 List of Tables 1 2 3 4 5 6 7 8 9 10 . . . . . . . . . . . . . . . . . . . Wave Parameter Conversion Table . . . . Actinic Flux Enhancement . . . . . . . . Surface Albedo . . . . . . . . . . . . . . Temperature Dependence of αp . . . . . . Pressure-Broadened Half Widths . . . . . Mechanical Analogues of Important Gases HITRAN database . . . . . . . . . . . . . Full-range Gaussian quadrature . . . . . . Half-range Gaussian quadrature . . . . . FACT s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i 4 13 52 77 77 84 93 157 157 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 F ASR and F OLR are the absorbed solar radiation and emitted longwave radiation, respectively, then ∂E = F ASR − F OLR (1) ∂t On timescales longer than about a year the Earth as a whole is thought to be in planetary radiative equilibrium. That, is, the global annual mean planetary temperature is nearly constant because the 1 INTRODUCTION 2 absorbed solar energy is exactly compensated by thermal radiation lost to space over the course of a year. Thus F ASR = F OLR (2) The total amount of solar energy available for the Earth to absorb is the incoming solar flux (or irradiance) at the top of Earth’s atmosphere, F (aka the solar constant), times the intercepting area 2 2 of Earth’s disk which is πr⊕ . Since Earth rotates, the total mean incident flux πr⊕ F is actually distributed over the entire surface area of the Earth. The surface area of a sphere is four times its cross-sectional area so the mean incident flux per unit surface area is F /4. The fraction of incident solar flux which is reflected back to space, and thus unable to heat the planet, is called the aplanetary albedo or spherical albedo, R. Satellite observations show that R ≈ 0.3. Thus only (1 − R) of the mean incident solar flux contributes to warming the planet and we have F ASR = (1 − R)F /4 (3) Earth does not cool to space as a perfect blackbody (41a) of a single temperature and emissivity. Nevertheless the spectrum of thermal radiation F OLR 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 TE ≡ (F OLR /σ)1/4 (4) For a perfect blackbody, T = TE . For a planet, the difference between TE and the mean surface temperature Ts is due to the radiative effects of the overlying atmosphere. The insulating behavior of the atmosphere is more commonly known as the greenhouse effect. Substituting (3) and (4) into (2) 4 (1 − R)F /4 = σTE (5) 1/4 TE = (1 − R)F 4σ (6) For Earth, R ≈ 0.3 and F ≈ 1367 W m−2 . Using these values in (6) yields TE = 255 K. Observations show the mean surface temperature Ts = 288 K. 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 E = hν (7) 1 INTRODUCTION 3 where h is Planck’s constant. Regrettably, almost no radiative transfer literature expresses quantities in 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 ω = 2πν (8) Since radians are considered dimensionless, the units of ω are s−1 . However, angular frequency is also rarely used in radiative transfer. Thus some authors use the symbol ω to denote the element of solid angle, as in dω. The reader should be careful not to misconstrue the two meanings. In this text we use ω only infrequently. 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 λν = c where c is the speed of light. Since c depends on the medium, λ also depends on the medium. The wavenumber ν m−1 , is exactly the inverse of wavelength ˜ ν ≡ ˜ ν 1 = λ c (10) (9) Thus ν measures the number of oscillations per unit distance, i.e., the number of wavecrests per ˜ meter. Using (9) in (10) we find ν = ν/c so wavenumber ν is indeed proportional to frequency ˜ ˜ (and thus to energy). Historically spectroscopists have favored ν rather than λ or ν. Because of this ˜ history, it is much more common in the literature to find ν expressed in CGS units of cm−1 than ˜ in SI units of m−1 . The CGS wavenumber is used analogously to frequency and to wavelength, i.e., to identify spectral regions. The energy of radiative transitions are commonly expressed in CGS wavenumber units. The relation between ν expressed in CGS wavenumber units (cm−1 ) and ˜ energy in SI units (J) is obtained by using (10) in (7) E = 100hcν (11) 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. k ≡ 2π˜ ν (12) The wavenumber k is set in Roman typeface as an additional distinction between it and other symbols 1 . Table 2 summarizes the relationships between the fundamental parameters which describe wave-like phenomena. 1 The script k is already used for Boltzmann’s constant, absorption coefficients, and vibrational modes 2 RADIATIVE TRANSFER EQUATION Table 2: Wave Parameter Conversion Tableab Variable ν λ ν ˜ ω k Units s−1 m cm−1 s−1 m−1 c ν 2πν ν − 2πν ν 100c c λ ν ˜ ω k τ c λ 100c˜ ν ω 2π kc 2π 1 τ − 1 100˜ ν 2πc ω 2π k cτ 1 100λ − ω 200πc k 200π 1 100cτ 2πc λ ν ˜ 200πc − ck 2π τ 2π λ 200π˜ ν ω c − 2π cτ 4 τ s 1 ν λ c 1 100c˜ ν 2π ω 2π ck − a b The speed of light is c m s−1 . Table entries express the column in terms of the row. 2 2.1 2.1.1 Radiative Transfer Equation Definitions 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 ˆ dE = Iν (r, Ω, t, ν) dS dt dΩ dν (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 n, then we ˆ 2 RADIATIVE TRANSFER EQUATION 5 ˆ must alter (13) to account for projection of dS onto dA. If the angle between n and Ω is θ then ˆ cos θ = n · Ω ˆ ˆ and the projection of dS onto dA yields dA = cos θ dS so that ˆ dE = Iν (r, Ω, t, ν) cos θ dA dt dΩ dν (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 n. 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 n 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 (15) (14) ¯ The mean intensity is an integrated measure of the radiation field at any point r. Mean intensity Iν is defined as the mean value of the intensity field integrated over all angles. ¯ Iν = I Ω ν dΩ dΩ Ω (17) 2 RADIATIVE TRANSFER EQUATION 6 The solid angle subtended by Ω is the ratio of the area A enclosed by Ω on a spherical surface to the square of the radius of the sphere. Since the area of a sphere is 4πr2 , there must be 4π steradians in a sphere. It is straightforward to demonstrate that the differential element of area in spherical polar coordinates is r2 sin θ dθ dφ. Thus the element of solid angle is Ω = A/r2 dΩ = r−2 dA = r−2 r2 sin θ dθ dφ = sin θ dθ dφ (18) The field of view of an instrument, e.g., a telescope, is most naturally measured by a solid angle. ¯ The definition of Iν (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 θ=π φ=2π dΩ = Ω θ=0 φ=0 θ=π sin θ dθ dφ sin θ dθ θ=0 θ=π = [φ]2π 0 = 2π θ=0 sin θ dθ = 2π [− cos θ]π 0 = 2π[−(−1) − (−1)] = 4π (19) As expected, there are 4π steradians in a sphere, and 2π steradians in a hemisphere. 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 ¯ Iν = 1 4π Iν dΩ Ω (20) ¯ Iν has units of radiance, W m−2 sr−1 Hz−1 . If the radiation field is azimuthally independent (i.e., Iν does not depend on φ), then 1 ¯ Iν = 2 π Iν sin θ dθ 0 (21) Let us simplify (21) by introducing the change of variables u = cos θ du = − sin θ dθ (22) (23) 2 RADIATIVE TRANSFER EQUATION 7 This maps θ ∈ [0, π] into u ∈ [1, −1] so that (21) becomes 1 ¯ Iν du Iν = − 2 1 1 1 ¯ Iν = Iν du 2 −1 −1 (24) The hemispheric intensities or half-range intensities are simply the up- and downwelling components of which the full intensity is composed ˆ Iν (τ, Ω) = Iν (τ, θ, φ) = ˆ Iν (τ, Ω) = Iν (τ, u, φ) = + Iν (τ, θ, φ) : 0 < θ < π/2 − Iν (τ, θ, φ) : π/2 < θ < π + Iν (τ, u, φ) − Iν (τ, u, φ) (25a) (25b) : 0≤u<1 : −1 < u < 0 + Iν (τ, µ, φ) = Iν (τ, +µ, φ) = Iν (τ, 0 < θ < π/2, φ) = Iν (τ, 0 < u < 1, φ) − Iν (τ, µ, φ) = Iν (τ, +µ, φ) = Iν (τ, π/2 < θ < π, φ) = Iν (τ, −1 < u < 0, φ) (26a) (26b) 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 ˆ · Ω , thus ˆ efficiency cos Θ = Ω Fν = = θ=0 φ=0 Iν cos θ dΩ Ω θ=π φ=2π Iν cos θ sin θ dθ dφ (27) In a plane-parallel medium, this defines the net specific irradiance passing through a given vertical level. Note the similarity between (20) and (27). The former contains the zeroth moment of the intensity with respect to the cosine of the polar angle, the latter contains the first moment. Also note 0 that (27) integrates the cosine-weighted radiance over all angles. If Iν is isotropic, i.e., Iν = Iν , then Fν = 0 due to the symmetry of cos θ. Let us simplify (27) by introducing the change of variables u = cos θ, du = − sin θ dθ. This maps θ ∈ [0, π] into u ∈ [1, −1]: u=−1 φ=2π Fν = u=1 u=1 φ=0 φ=2π Iν u (−du) dφ Iν u du dφ u=−1 φ=0 = (28) 2 RADIATIVE TRANSFER EQUATION 8 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 Fν dν = −Fλ dλ dλ Fν = −Fλ dν d c = −Fλ dν ν c = −Fλ − 2 ν c λ2 Fν = 2 Fλ = Fλ ν c c ν2 Fλ = Fν = Fν λ2 c Thus Fν and Fλ are always of the same sign. 2.1.4 Actinic Flux (29) (30) A quantity of great importance in photochemistry is the total convergence of radiation at a point. This quantity, called the actinic flux, F J , 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 J Fν = 4π Iν dΩ (31) ¯ = 4πIν J ¯ Thus the actinic flux is simply 4π times the mean intensity Iν (20). Fν has units of W m−2 Hz−1 which are identical to the units of irradiance Fν (27). Although the nomenclature “actinic flux” is J somewhat appropriate, it is also somewhat ambiguous. The “flux” measured by Fν at a point P is the energy convergence (per unit time, frequency, and area) through the surface of the sphere containing P . This differs from the “flux” measured by Fν , which is the net energy transport (per unit time, frequency, and area) through a defined horizontal surface. Thus it is safest to use J the terms “actinic radiation field” for Fν and “irradiance” for Fν . Unfortunately the literature is permeated with the ambiguous terms “actinic flux” and “flux”, respectively. J The usefulness of actinic flux Fν becomes apparent only in conjunction with additional, speciesdependent data describing the probability of photon absorption, or photo-absorption. Photoabsorption 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 J field Fν , and the efficiency with with each molecule absorbs photons. This efficiency is called the 2 RADIATIVE TRANSFER EQUATION 9 absorption cross-section, molecular cross section, or simply cross-section. The absorption crosssection 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 J a photon with frequency in [ν, ν + dν] is proportional to Fν (ν)α(ν)2 . 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 α J Fν (ν) = NA Fν (ν)α(ν) (32) where NA m−3 is the number concentration of A. 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 AB + hν −→ A + B ν>ν0 AB + hν −→ A + B 0 AB + hν −→ A + B λ<λ (33) Both forms indicate that the efficiency with which reaction (33) proceeds is a function of photon energy. The second form makes explicit that the photodissociation reaction does not proceed unless ν < ν0 , where ν0 is the photolysis cutoff frequency. In any case, photon energy is conventionally written hν, rather than the less convenient hc/λ. 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 dimensionless3 . 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 . Jν = J Fν (ν)α(ν)φ(ν) hν ¯ν (ν)α(ν)φ(ν) 4πI = hν (34) J The proportionality constant is (hν)−1 , i.e., when the actinic radiation field Fν (ν) is converted from energy −1 −2 −1 −1 J (J m s Hz ) to photons (# m s Hz ) then Fν (ν)α(ν) is the probability of absorption per second per unit frequency 3 Azimuthal angle and quantum yield are both dimensionless quantities denoted by φ. The meaning of φ should be clear from the context. 2 −2 −1 2 RADIATIVE TRANSFER EQUATION 10 J The photon energy in the denominator converts the energy per unit area in Fν to units of photons per unit area. The factor of α turns this photon flux into a photo-absorption rate per unit area. The final factor, φ, converts the photo-absorption rate into a photodissociation rate coefficient. Note that each factor in the numerator of (34) requires detailed spectral knowledge, either of the radiation field or of the photochemical behavior of the molecule in question. This complexity is a hallmark of atmospheric photochemistry. The total photolysis rate coefficient J is obtained by integrating (34) over all frequencies which may contribute to photodissociation J = ν>ν0 Jν (ν) dν J Fν (ν)α(ν)φ(ν) dν hν ν>ν0 ¯ 4π Iν (ν)α(ν)φ(ν) dν = h ν>ν0 ν = (35) As mentioned above, evaluation of (35) requires essentially a complete knowledge of the radiative and photochemical properties of the environment and species of interest. Moreover, J is notoriously difficult to compute where there is any uncertainty in the input quantities. 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 J of the difficulty is due to the steep but opposite gradients of Fν 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 NO2 + hν λ<420 nm −→ NO + O (36) If [NO2 ] denotes the number concentration of NO2 in a closed system where photolysis is the only sink of NO2 , then d[NO2 ] = −JNO2 [NO2 ] + SNO2 dt (37) where SNO2 represents all sources of NO2 . The terms in (37) all have dimensions of m−3 s−1 . The first term on the RHS is the photolysis rate of NO2 in the system. 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 2 shows the vertical distribution of JNO2 [s−1 ] for the conditions shown in Figure (1). J In this section we have assumed the quantities Fν , α, 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 2 RADIATIVE TRANSFER EQUATION 11 Figure 1: (a) Spectral distribution of actinic flux F J [# 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). 2 RADIATIVE TRANSFER EQUATION 12 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. 2 RADIATIVE TRANSFER EQUATION Table 3: Actinic Flux Enhancement by Scatteringa − + ˆ Fν ρL (ν) Iν (Ω) Fν Fν Fν Fν Fν 0 0 1 1 ˆ ˆ Fν δ(Ω − Ω ) − Iν = Fν /π + Iν = 0 − ˆ ˆ Iν = Fν δ(Ω − Ω ) + Iν = Fν /π 13 Description Collimated, non-reflecting Isotropic, non-reflecting Collimated, reflecting 0 0 Fν Fν ¯ Iν Fν 4π Fν 2π 3Fν 4π Fν π J Fν Fν 2Fν 3Fν 4Fν Isotropic, reflecting a Fν /π − Terminology: Fν is downwelling irradiance of source, ρL (ν) is Lambertian reflectance of surface (or cloud), + J ¯ ˆ is resulting intensity field, Fν is upwelling irradiance, Iν is mean intensity, Fν is actinic flux. Iν (Ω) 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 J The actinic flux Fν (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. 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 J J order of increasing actinic flux Fν , shown in the final column. Fν is seen to increase with both the number and the brightness of reflecting surfaces. The reflectivity of natural surfaces is no more J than 90% nor less than 5%4 , so that the Fν 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. J The actinic flux Fν = 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 Glaciers are the most reflective surfaces in nature, with ρL (ν) < 0.9. Maximum cloud reflectance is < 0.7. Dark ∼ ∼ forests and ocean have ρL (ν) > 0.05. See discussion in §2.3.2 and Table 4. ∼ 4 2 RADIATIVE TRANSFER EQUATION 14 incident irradiance J − S ≡ Fν /Fν (38) Table 3 shows S ranges from one for the direct solar beam to a maximum of four in a completely isotropic radiation field. Thus a collimated beam is the least efficient configuration of radiant energy for driving photochemistry. Scattering processes diffuse the radiation field, and, as a result, always enhance the photolytic efficiency of a given irradiance. Absorption in a medium (i.e., in J the atmosphere or by the surface) always reduces Fν and may lead to S < 1. 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 J clouds depends on whether the actinic flux enhancement factor Fν = 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 E ≡ J Fν − 4Fν (39) Clearly 0 ≤ E ≤ 1. Table 3 shows E = 0.25 for a collimated beam, and E = 1 for isotropic radiation. Just as scattering of the solar beam is required to increase E above 0.25, absorption must be present to reduce E beneath 0.25. 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 J that photochemistry itself is driven by molecular absorption which reduces Fν 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, 2 RADIATIVE TRANSFER EQUATION 15 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 J ¯ flux Fν (31) and thus to the mean intensity Iν . Uν = 4π dUν 4π ¯ Iν c (40) = The units of Uν are J m−3 Hz−1 . 2.1.7 Spectral vs. Broadband Until now we have considered only spectrally dependent quantities such as the spectral radiance Iν , J spectral irradiance Fν , spectral actinic flux Fν , 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: ∞ I = 0 ∞ Iν (ν) dν Fν (ν) dν 0 ∞ J Fν (ν) dν 0 ∞ α Fν (ν) dν 0 ∞ F = FJ = Fα = J = 0 ∞ Jν (ν) dν Uν (ν) dν 0 ∞ U = h = 0 hν (ν) dν 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 2 RADIATIVE TRANSFER EQUATION 16 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 c2 (ehν/kT − 1) 2hc2 Bλ (T, λ) = 5 hc/λkT λ (e − 1) Bν (T, ν) = 2hν 3 (41a) (41b) Blackbodies emit isotropically—this considerably simplifies thermal radiative transfer. The correct predictions (41a) resolved one of the great mysteries in experimental physics in the late 19th century. In fact, this discovery marked the beginning of the science of quantum mechanics. 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 Bν dν = −Bλ dλ (42) Once again, the negative sign arises as a result of the opposite senses of increasing frequency 2 RADIATIVE TRANSFER EQUATION 17 versus increasing wavelength. We may derive (41b) from (41a) by using (9) in (42) Bν = −Bλ = Bλ c ν2 λ2 c2 c Bλ ν2 c Bν λ2 (43) (44) dλ dν = Bλ c Bν = Bλ λ2 Bλ = c ν2 Bν = = c These relations are analogous to (30). 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 hνM In Wien’s limit, (41) becomes 2hν 3 −hν/kT e c2 2hc2 Bλ (T, λ) = 5 e−hc/λkT λ Bν (T, ν) = (46a) (46b) kT (45) In the very low energy limit, known as the Rayleigh-Jeans limit, the photon energy is much less than the ambient thermal energy hνM kT (47) Thus in the Rayleigh-Jeans limit the arguments to the exponential in (41) are less than 1 so the exponentials may be expanded in Taylor series. Starting from (41a) hν 2hν 3 1+ −1 Bν (T, ν) ≈ 2 c kT 2hν 3 kT = c2 hν 2ν 2 kT = c2 Similar manipulation of (41b) may be performed and we obtain 2ν 2 kT c2 2ckT Bλ (T, λ) ≈ λ4 Bν (T, ν) ≈ (49a) (49b) −1 (48) 2 RADIATIVE TRANSFER EQUATION 18 The frequency of extreme emission is obtained by taking the partial derivative of (41a) with respect to frequency with the temperature held constant ∂Bν ∂ν = T h 2h 1 × hν/kT × 3ν 2 (ehν/kT − 1) − ν 3 ehν/kT 2 2 c (e − 1) kT To solve for the frequency of maximum emission, νM , we set the RHS equal to zero so that one or more of the LHS factors must equal zero 2 2hνM 3(ehνM /kT − 1) − hνM ehνM /kT kT 2 (ehνM /kT − 1)2 c hνM hνM /kT 3(ehνM /kT − 1) − e kT hνM hνM /kT e kT hνM kT = 0 = 0 = 3(ehνM /kT − 1) = 3(1 − e−hνM /kT ) (50) An analytic solution to (50) is impossible since νM cannot be factored out of this transcendental equation. If instead we solve x = 3(1 − e−x ) numerically we find that x ≈ 2.8215 so that hνM ≈ 2.82 kT νM ≈ 2.82kT /h ≈ 5.88 × 1010 T Hz (51) where the units of the numerical factor are Hz K−1 . Thus the frequency of peak blackbody emission is directly proportional to temperature. This is known as Wien’s Displacement Law. A separate relation may be derived for the wavelength of maximum emission λM by an analogous procedure starting from (41b). The result is λM ≈ 2897.8/T µm (52) where the units of the numerical factor are µm K. Note that (51) and (52) do not yield the same answer because they measure different quantities. The wavelength of maximum emission per unit wavelength, for example, is displaced by a factor of approximately 1.76 from the wavelength of maximum emission per unit frequency. 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 2 RADIATIVE TRANSFER EQUATION 19 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 u=1 Fν = 2π = 2π Iν u du u=−1 u=0 u=1 Iν u du + u=−1 u=0 Iν u du (53) We now introduce the change of variables µ = |u| = | cos θ|. Referring to (23) we find µ = | cos θ| = µ = |u| = cos θ − cos θ u −u : 0 < θ < π/2 : π/2 < θ < π (54a) (54b) : 0≤u<1 : −1 < u < 0 Most formal work on radiative transfer is written in terms of µ rather than u or θ. Substituting (54) into (53) µ=0 µ=1 Fν = 2π µ=1 µ=0 Iν (−µ) (−dµ) + µ=0 µ=1 Iν µ dµ = 2π µ=1 Iν µ dµ + µ=0 µ=1 Iν µ dµ µ=1 = 2π − = = µ=0 − + −Fν + Fν + − Fν − Fν Iν µ dµ + µ=0 Iν µ dµ (55) where we have defined the hemispheric fluxes or half-range fluxes 1 + Fν = 2π 0 1 − Fν = 2π 0 Iν (+µ)µ dµ Iν (−µ)µ dµ (56a) (56b) The hemispheric fluxes are positive definite, and their difference is the net flux. − + Fν = Fν − Fν (57) + − Fν = Fν − Fν (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. 2 RADIATIVE TRANSFER EQUATION 20 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 1 B Fν = 2π 0 1 Bν µ dµ µ dµ µ 2 0 2 µ=1 µ=0 = 2πBν = 2πBν = = πBν 2πBν ( 1 2 − 0) (59) Thus the upwelling blackbody irradiance tranports π times the constant intensity of the radiation. Given that the upper hemisphere contains 2π steradians, one might naively expect the upwelling irradiance to be 2πBν . In fact the divergence of blackbody radiance above the emitting surface is 2πBν . But the vertical flux of energy is obtained by cosine-weighting the radiance over the 1 hemisphere and this weight introduces the factor of 2 difference between the naive and the correct solutions. 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 ≡ 0 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 FB in terms fundamental physical constants and the temperature of the body. ∞ + FB = π 0 ∞ Bν dν 2hν 3 = π dν c2 (ehν/kT − 1) 2πh ∞ ν3 = dν c2 0 ehν/kT − 1 0 (60) 2 RADIATIVE TRANSFER EQUATION 21 To simplify (60) we make the change of variables hν kT kT x ν = h kT dν = dx h h dν dx = kT x = This change of variables maps ν ∈ [0, ∞) to x ∈ [0, ∞). Substituting this into (60) we obtain + FB = 2πh ∞ c2 0 2πk 4 T 4 = c2 h3 kT h ∞ 0 3 x3 kT dx ex − 1 h (61) x3 dx ex − 1 The definite integral in (61) is π4 /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 ζ(x) ≡ 1 Γ(x) ∞ 0 ux−1 du eu − 1 (62) Comparing (61) with the Riemann zeta function definition (62), we see that x = 4, i.e, + FB = 2πk 4 T 4 Γ(4)ζ(4) c2 h3 (63) 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: ux−1 e−u ux−1 = u × −u eu − 1 e −1 e ux−1 e−u = 1 − e−u = ux−1 e−u × 1 1 − e−u (64) The integration limits in (62) are [0, ∞) so it is always true that e−u < 1 in (64), and thus in the integrand of (62). Recall that the sum of an infinite power series with initial term a0 and ratio r is ∞ a0 r k = k=0 k→∞ lim a0 + a0 r + a0 r2 + · · · + a0 rk−1 + a0 rk (65) = a0 /(1 − r) for |r| < 1 2 RADIATIVE TRANSFER EQUATION 22 Hence the last term in (64) is the sum of a power series (65) with initial term a0 = 1 and ratio r = e−u . 1 = 1 − e−u Substituting (66) into (64) we obtain ux−1 = ux−1 e−u e−ku u−1 e k=0 ∞ ∞ ∞ e−ku k=0 (66) = k=0 ∞ ux−1 e−(k+1)u ux−1 e−ku k=1 = (67) where the last step shifts the initial index to from zero to one. Using the inifinite series representation (67) for the integrand of (62) yields 1 ζ(x) = Γ(x) ∞ 0 ∞ ux−1 e−ku du k=1 Integration and addition are commutative operations. Interchanging their order yields 1 ζ(x) = Γ(x) ∞ 0 ∞ ux−1 e−ku du k=1 (68) We change variables from u ∈ [0, +∞] to y ∈ [0, +∞] with y dy u du so that (68) becomes 1 ζ(x) = Γ(x) 1 = Γ(x) = = k=1 ∞ 0 ∞ = = = = ku k du y/k k −1 dy (69) k=1 ∞ y k k x−1 e−y k −1 dy ∞ k k=1 ∞ −(x−1) −1 0 y x−1 e−y dy 1 Γ(x) ∞ k −x [Γ(x)] k=1 k −x (70) 2 RADIATIVE TRANSFER EQUATION 23 where we replaced the integral in brackets with the Gamma function it defines. Hence the Riemann zeta function of a positive integer n is the sum of the reciprocals of the postive integers to the power n. Contour integration in the complex plane gives analytic closed-form solutions to Σ∞ k −n (70) 1 for positive, even integers n. Our immediate concern is n = 4. Consider the complex function [Carrier et al., 1983, p. 97] f (z) = This function 1. 2. 3. 4. Is analytic throughout the complex plane Has first order poles at all integer values on the real axis (except the origin) Has a fifth order pole at the origin Satisfies lim|R|→∞ f (z) = 0 where z = Reiθ π cot πz z4 (71) Therefore (71) obeys the residue theorem for suitably chosen contours. In other words, fxm. Having shown ζ(4) = π4 /90 Using Γ(4) = 3! = 3 × 2 × 1 = 6 and ζ(4) = π4 /90 from (72), Equation (63) becomes π4 2πk 4 ×6× × T4 2 h3 c 90 2π5 k 4 4 = T 15c2 h3 This is known as the Stefan-Boltzmann Law of radiation, and is usually written as + FB = + FB = σT 4 where (73) 5 4 2π k (74) σ ≡ 15c2 h3 σ is known as the Stefan-Boltzmann constant and depends only on fundamental physical constants. The value of σ is 5.67032 × 10−8 W m−2 K−4 . The thermal emission of matter depends very strongly (quartically) on T (73). This rather surprising result has profound implications for Earth’s climate. + We derived FB (73) directly so that the Stefan-Boltzmann constant would fall naturally from the derivation. For completeness we now present the broadband blackbody intensity B ∞ (72) B = 0 Bν (ν) dν 2π4 k 4 4 T 15c2 h3 + = FB /π = σT 4 /π = (75) + The factor of π difference between B and FB is at first confusing. One must remember that B is the + broadband (spectrally-integrated) intensity and that FB is the broadband hemispheric irradiance (spectrally-and-angularly-integrated). 2 RADIATIVE TRANSFER EQUATION Luminosity 24 2.1.12 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. L = 4πD2 F (76) The meaning of the solar constant F is made clear by (76). F is the solar irradiance that would be measured normal to the Earth-Sun axis at the top of Earth’s atmosphere. The mean EarthSun distance is ∼ 1.5 × 1011 m and F ≈ 1367 W m−2 so L = 3.9 × 1026 W. Note that L has dimensions of power, i.e., J s−1 . 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) 4πr2 σT 4 = 4πD2 F D2 F r2 = σT 4 r = D T2 F σ (77) For the Sun-Earth system, T ≈ 5800 K so r = 6.9 × 108 m. Solar radiation received by Earth appears to originate from a portion of the solar atmosphere known as the photosphere. 2.1.13 Extinction and Emission Radiation and matter have only two forms of interactions, extinction and emission. Lambert5 , 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 dIν = −k(ν)Iν ds Extinction only (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 Beer is usually credited with first formulating this law, but actually the similarly-name Bougher deserves the credit. 5 2 RADIATIVE TRANSFER EQUATION 25 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 interrelated by κe [m2 molecule−1 ] ρ [kg m−3 ] ψe [m2 kg−1 ] × M [kg mol−1 ] κe [m2 molecule−1 ] = N [molecule mol−1 ] ke [m−1 ] = κe [m2 molecule−1 ] × n [molecule m−3 ] ψe [m2 kg−1 ] = (79) 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 dIν = k(ν)Sν Emission only (80) ds 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, dIν ds 1 dIν ke ds = −ke Iν + ke Sν = −Iν + Sν (81) Equation (81) is the equation of radiative transfer in its simplest differential form. 2.1.14 Optical Depth We define the optical path τ between points P1 and P2 as ˜ P2 τ (P1 , P2 ) = ˜ P1 ke ds (82) d˜ = ke ds τ 2 RADIATIVE TRANSFER EQUATION 26 The optical path measures the amount of extinction a beam of light experiences traveling between two points. When τ > 1, the path is said to be optically thick. ˜ 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 z2 > z1 , and then allow z2 → ∞ z =z2 τ (z1 , z2 ) = z =z1 z =∞ ke dz ke dz z =z z =∞ τ (z, ∞) = = z =z (83) (84) ke dz Equation (84) is the integral definition of optical depth. The differential definition of optical depth is obtained by differentiating (84) with respect to the lower limit of integration and using the fundamental theorem of differential calculus dτ = [ke (∞) − ke (z)] dz = −ke (z) dz (85) where the second step uses the convention that ke (∞) = 0. By convention, τ is positive definite, but (85) shows that dτ may be positive or negative. If this seems confusing, consider the analogy with atmospheric pressure: pressure increases monotonically from zero at the top of the atmosphere, and we often express physical concepts such as the temperature lapse rate in terms of negative pressure gradients. 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 scattering6 . 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. The effects of diffraction will not be explicitly included either, but are implicit in the assumption of the horizontal homogeneity of the radiative flux. 6 2 RADIATIVE TRANSFER EQUATION 27 Since each particle absorbs incident energy in proportion to its geometric cross-section, the total absorption of radiation per particle is proportional to πr2 Qe , where Qe = 1 for perfectly absorbing particles. If Qe < 1, then each particle removes fewer photons than suggested by its geometric size7 . Conversely, if Qe > 1 each particle removes more photons than suggested by its geometric size. Each particle encountered removes πr2 Qe F (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 ∆F = −πr2 Qe F N ∆z ∆F = −πr2 Qe N ∆z F If we take the limit as ∆z → 0 then (86) becomes 1 dF = −πr2 Qe N dz F d(ln F ) = −πr2 Qe N dz Let us define the volume extinction coefficient k as k ≡ πr2 Qe N Finally we define the optical depth τ in terms of the extinction τ ≡ k∆z τ = πr2 Qe N ∆z (89) (90) (88) (86) (87) The dimensions of k are m−1 and therefore τ is dimensionless. All quantities which compose τ are positive by convention, therefore τ itself is positive definite. 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 τ F (τ ) = F e−τ (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. r 7 lim Qe = 2 λ (92) Although intuition suggests a spherical particle should not remove more energy from a collimated light beam than it can geometrically intercept, this is not the case. As discussed later, diffraction around the particle is important. In fact, for r λ, Qe → 2. However, most of the extinction due to diffraction occurs as scattering, not absorption. 2 RADIATIVE TRANSFER EQUATION 28 Thus particles larger than 5–10 µm extinguish twice as much visible light as their geometric crosssection 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 4 (93) M = πr3 ρN 3 If we substitute 3M πr2 N = 4rρ into (90) we obtain τ= 3M Qe ∆z 4rρ (94) Typical cloud particles have r ∼ 10 µm so that for visible solar radiation with λ ∼ 0.5 µm we may employ (92) to obtain 3M ∆z τ= (95) 2rρ 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) dIν d˜ τ = −Iν + Sν (96) A stratified atmosphere is one in which all atmospheric properties, e.g., temperature, density, vary only in the vertical. As shown in the non-existant Figure, the photon path increment ds at polar angle θ in a stratified atmosphere is related to the vertical path increment dz by dz = cos θ ds or ds = u−1 dz. In other words, the optical path traversed by photons is proportional to the vertical path divided by the cosine of the trajectory. d˜ = u−1 dτ τ dτ = u d˜ τ Substituting (97a) into (96) we obtain u dIν dτ = −Iν + Sν (98) (97a) (97b) This is the differential form of the radiative transfer equation in a plane parallel atmosphere and is valid for all angles. The solution of (98) is made difficult because Sν depends on Iν . A more tractable set of equations may be obtained by considering the form of the boundary conditions. For many (most) problems of atmospheric interest, we know Iν over an entire hemisphere at each boundary of a “slab”. Considering the entire atmosphere as a slab, for example, we know that, 2 RADIATIVE TRANSFER EQUATION 29 at the top of the atmosphere, sunlight is the only incident intensity from the hemisphere containing the sun. Or, at the surface, we have constraints on the upwelling intensity due to thermal emission or the surface reflectivity. When combined, these two hemispheric boundary conditions span a complete range of polar angle, and are thus sufficient to solve (98). However, in practice it is difficult to apply half a boundary condition. Moreover, we are often interested in knowing the hemispheric flows of radiation because many instruments (e.g., pyranometers) are designed to measure hemispheric irradiance and many models (e.g., climate models) require hemispheric irradiance to compute surface exchange properties. For these reasons we will decouple (98) into its constituent upwelling and downwelling radiation components. Using the definitions of the half-range intensities (25) in (98) we obtain − dIν − − = Iν − Sν 0 < θ < π/2, u = cos θ > 0 (99a) dτ dI + + + −u ν = Iν − Sν π/2 < θ < π, u = cos θ < 0 (99b) dτ where we have simply multiplied (98) by −1 in order to place the negative sign on the LHS for − + reasons that will be explained shortly. The definitions of Sν and Sν are exactly analogous to (25). 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 −u − dIν − − = Iν − Sν (100a) dτ dI + + + µ ν = Iν − Sν (100b) dτ These are the equations of radiative transfer in slab geometry for downwelling (0 < θ < π/2) and upwelling (π/2 < θ < π) intensities, respectively. The only mathematical difference between (100a) and (100b) is the negative sign. A helpful mnemonic is that the negative sign on the LHS + − is associated with Iν while the implicit unary positive sign is associated with Iν . Of course the Iν and Sν terms on the RHS are prefixed with opposited signs since they represent opposing, but positive definite, physical processes (absorption and emission). 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 −µ dIν = −Iν + Bν (101) dτ 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. u 2 RADIATIVE TRANSFER EQUATION 30 2.2 2.2.1 Integral Equations 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ν . dIν dτ u dIν u dIν + Iν dτ Iν dIν + dτ u u = −Iν + Sν = −Iν dτ + Sν dτ = Sν dτ = u−1 Sν dτ (102) Multiplying (102) by the integrating factor eτ /u eτ /u dIν + eτ /u Iν dτ = u−1 eτ /u Sν dτ u d(eτ /u Iν ) = u−1 eτ /u Sν dτ d(eτ /u Iν ) = u−1 eτ /u Sν dτ (103) The LHS side of (103) is a complete differential. The boundary condition which applies to this first degree differential equation depends on the direction the radiation is traveling. Thus we denote the solutions to (103) as Iν (τ, +µ) and Iν (τ, −µ) for upwelling and downwelling radiances, respectively. We shall assume that the upwelling intensity at the surface, Iν (τ ∗ , +µ), and the downwelling intensity at the top of the atmosphere, Iν (0, −µ), are known quantities. Since µ is positive definite (54), +µ and −µ uniquely specify the angles for which these boundary conditions apply. The solution for upwelling radiance is obtained by replacing u in (103) by −µ because u < 0 for upwelling intensities. + d(e−τ /µ Iν ) + = −µ−1 e−τ /µ Sν dτ (104) We could have arrived at (104) by starting from (100b), and proceeding as above except using e−τ /µ as the integrating factor. We now integrate from the surface to level τ (i.e., along a path of decreasing τ ) and apply the boundary condition at the surface e−τ /µ τ =τ τ =τ Iν (τ , +µ) −τ ∗ /µ τ =τ ∗ = −µ−1 τ τ∗ τ =τ ∗ e−τ e−τ ∗ /µ /µ + Sν dτ e−τ /µ Iν (τ, +µ) − e Iν (τ ∗ , +µ) = µ−1 Iν (τ, +µ) = e −τ ∗ /µ + Sν dτ τ∗ −1 τ e −τ /µ Iν (τ , +µ) + µ ∗ e−τ τ∗ −1 τ /µ + Sν dτ )/µ + Sν dτ Iν (τ, +µ) = e (τ −τ ∗ )/µ Iν (τ , +µ) + µ e(τ −τ 2 RADIATIVE TRANSFER EQUATION 31 Note that τ ∗ > τ and τ > τ so that both of the transmission factors reduce a beam’s intensity between its source (at τ ∗ or τ ) and where it is measured (at τ ). The physical meaning of the transmission factors is more clear if we write all transmission factors as negative exponentials. Iν (τ, +µ) = e −(τ ∗ −τ )/µ τ∗ Iν (τ , +µ) + µ ∗ −1 τ e−(τ −τ )/µ Sν (τ , +µ) dτ (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. − d(eτ /µ Iν ) − = µ−1 eτ /µ Sν dτ eτ /µ τ =τ τ =τ Iν (τ , −µ) = µ−1 e eτ τ =0 τ τ /µ 0 /µ τ =0 − Sν dτ eτ /µ Iν (τ, −µ) − Iν (0, −µ) = µ−1 e τ /µ − Sν dτ τ −1 0 τ Iν (τ, −µ) = Iν (0, −µ) + µ Iν (τ, −µ) = e −τ /µ eτ −1 /µ − Sν dτ Iν (0, −µ) + µ e(−τ +τ 0 τ )/µ − Sν dτ Iν (τ, −µ) = e−τ /µ Iν (0, −µ) + µ−1 0 e−(τ −τ )/µ Sν (τ , −µ) dτ (106) 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 2 RADIATIVE TRANSFER EQUATION 32 blackbody radiation into the upper hemisphere Iν (τ ∗ , +µ) = Bν [T (τ ∗ )] (107) ∗ For brevity we shall define Bν = Bν [T (τ ∗ )]. At the top of the atmosphere, we assume there is no downwelling thermal radiation. Iν (0, −µ) = 0 (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 Iν (τ, +µ) = e −(τ ∗ −τ )/µ τ τ∗ Bν (τ ) + µ )/µ ∗ −1 τ e−(τ −τ )/µ Bν (τ ) dτ (109) (110) Iν (τ, −µ) = µ−1 0 e−(τ −τ Bν (τ ) dτ The first term on the RHS of (109) is the thermal radiation emitted by the surface, attenuated by absorption in the atmosphere until it contributes to the upwelling intensity at level τ . The second term on the RHS contains the upwelling intensity arriving at τ contributed from the attenuated atmospheric thermal emission from each parcel between the surface and τ . The µ−1 factor in front of the integral accounts for the slant path of the thermally emitting atmosphere. The RHS of (110) is similar but contains no boundary contribution since the vacuum above the atmosphere is assumed to emit no thermal radiation. The upwelling and downwelling intensities in a stratified, thermal atmosphere are fully described by (109) and (110). 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). µ=1 − Fν (τ ) = 2π µ=0 Iν (τ, −µ)µ dµ (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 µ=1 − Fν (τ ) = 2π µ=0 τ =τ µ=1 τ =τ µ−1 τ =0 e−(τ −τ e−(τ −τ )/µ )/µ Bν (τ ) dτ µ dµ = 2π τ =0 τ =τ µ=0 Bν (τ ) dµ dτ )/µ µ=1 = 2π τ =0 Bν (τ ) µ=0 e−(τ −τ dµ dτ (112) Notice that two factors of µ cancelled each other out: The reduction in irradiance due to nonnormal incidence (µ) exactly compensates the increased irradiance due to emission by the entire slant column which is µ−1 times greater than emission from a vertical column. 2 RADIATIVE TRANSFER EQUATION 33 In terms of exponential integrals defined in Appendix 10.7, the inner integral in parentheses in (112) is E2 (τ − τ ) (603c). τ =τ − Fν (τ ) = 2π τ =0 Bν (τ )E2 (τ − τ ) dτ (113) Similar terms arise when we consider the horizontal upwelling flux obtained by substituting (109) into (56a) and we obtain µ=1 + Fν (τ ) = 2π µ=0 µ=1 Iν (τ, +µ)µ dµ e µ=0 µ=1 ∗ −(τ ∗ −τ )/µ τ∗ = 2π = 2πBν (τ ) Bν (τ ) + µ ∗ −1 τ e−(τ τ∗ −τ )/µ Bν (τ ) dτ µ=1 µ dµ −τ )/µ e µ=0 −(τ ∗ −τ )/µ µ dµ + 2π τ τ =τ Bν (τ ) µ=0 e−(τ dµ dτ (114) = 2πBν (τ ∗ )E3 (τ ∗ − τ ) + 2π τ =τ ∗ Bν (τ )E2 (τ − τ ) dτ Subtracting (113) from (114) we obtain the net flux at any layer in a non-scattering, thermal, stratified atmosphere + − Fν (τ ) = Fν (τ ) − Fν (τ ) τ∗ τ = 2π Bν (τ )E3 (τ − τ ) + τ ∗ ∗ Bν (τ )E2 (τ − τ ) dτ − 0 Bν (τ )E2 (τ − τ ) dτ (115) Equations (113) and (114) may not seem useful at this point but their utility becomes apparent in §2.3.4 where we define the flux transmissivity. 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 1 + Fν = 2π 0 1 Bν (T )µ dµ µ dµ 0 1 = 2πBν (T ) + Fν µ2 = 2πBν (T ) 2 0 1 = 2πBν (T ) −0 2 = πBν (T ) (116) 1 + F = Bν (T ) π ν 2 RADIATIVE TRANSFER EQUATION 34 Notice the isotropy of the Planck function allows the factor of 2 from the azimuthal integration to cancel the mean value of the cosine weighting function over a hemisphere. Moving the remaining factor of π from the azimuthal integration to the LHS conveniently sets the RHS equal to the Planck function. We integrate (116) over frequency to obtain the total upwelling thermal irradiance 1 π ∞ + Fν dν = 0 0 ∞ Bν (T ) dν ∞ 1 + Bν (T ) dν F = π 0 1 + σT 4 F = π π + F = σT 4 (117) Thus the factor of π from the azimuthal integration nicely cancels the factor of π from the StefanBoltzmann Law. Equation (117) applies to any surface whose emissivity is 1. Consider, e.g., a thick cloud with cloud base and cloud top temperatures T (zB ) = TB and T (zT ) = TT , respectively. Then the upwelling thermal flux at cloud top and the downwelling flux at cloud bottom will be 4 4 F − (zB ) = σTB and F + (zT ) = σTT , respectively. 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 2E3 (τ ∗ − τ ) = exp − τ∗ − τ µ ¯ (118) Inspection (or differentiation) shows that the atmosphere within one optical depth makes the dominant contribution to (118). For ∆τ = τ ∗ − τ = 1, the diffusivity factor µ−1 ≈ ¯ 5 3 (119) The hypothetical collimated beam of radiation is inclined to the zenith by arccos(3/5) ≈ 53.13◦ . This is equivalent to a collimated beam of radiation travelling vertically through an optical depth equal to five thirds the vertical optical depth traversed by the diffuse radiation. 2 RADIATIVE TRANSFER EQUATION 35 Using this assumption (118), the upwelling hemispheric irradiance (56a) for blackbody radiation is 1 + Fν = 2π 0 1 Iν (+µ)µ dµ Iν (+¯)µ dµ µ 0 1 ≈ 2π = 2π¯Iν (+¯) µ µ 0 2 µ dµ 1 µ 2 0 1 −0 = 2πIν (+¯) µ 2 = πIν (+¯) µ = 2πIν (+¯) µ (120) An analogous relationship holds for the downwelling irradiance. 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) ¯ −¯ µ − dFν − B = Fν − πFν dτ dF + + B µ ν = Fν − πFν ¯ dτ (121a) (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˜ = µ−1 dτ = τ ¯ 5 dτ . Second, we drop the frequency dependence, which is equivalent to integrating over a broad 3 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. − dF − = F − − πB d˜ τ dF + = F + − πB d˜ τ (122a) (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 rate8 or cooling rate. The Conventionally h is called the heating rate, which is a misnomer since the dimensions of dT /dt are K s−1 . Thus we use the less common but more accurate terminology “warming rate”. We reserve “heating rate” for a measure power dissipation, i.e., energy per unit time, in J s−1 , J m−3 s−1 , or J kg−1 s−1 . 8 2 RADIATIVE TRANSFER EQUATION 36 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 ] h≡ dT dt = 1 dF ρcp dz kg J 3 × m kg K −1 (123) × 1 J 2 × m s m K = s The forcing term on the RHS, dF/dz is the radiative flux divergence, the vertical gradient of net radiative flux. Absorption and emission are the only mechanisms which contribute to the flux divergence. In terms of hemispheric fluxes, dF dz = d (F − − F + ) dz (124) 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 F − − F + = F0 The net radiative flux F (z) = F0 is constant in radiative equilibrium. Adding and subtracting (122), we obtain d (F + − F − ) = F + + F − − 2πB d˜ τ d (F + + F − ) = F + − F − = F0 d˜ τ Defining ψ = F + − F − and φ = F = F + − F − , dφ = ψ − 2πB d˜ τ dψ = φ = F0 d˜ τ Since φ = F + − F − is constant, (127b) shows that dφ/d˜ = 0 and thus τ ψ = 2πB Substituting this into (127b), d (2πB) = F0 d˜ τ dB F0 = d˜ τ 2π F0 τ ˜ B(˜) = τ + C1 2π (128) (127a) (127b) (126a) (126b) (125) (129) 2 RADIATIVE TRANSFER EQUATION 37 We evaluate the constant of integration by using the boundary condition at the top of the atmosphere. By definition F − = 0 and τ = 0 at TOA. This implies that φ = ψ(0) = F + (0) = F0 . We ˜ must therefore have B(0) = φ/2π by (128). Using this result in (133), B(0) = C1 F0 (0) φ = + C1 2π 2π φ F0 = = 2π 2π (130) Substituting (130) back into (133) shows F0 τ ˜ F0 + 2π 2π F0 (˜ + 1) τ (131) = 2π The thermal absorption and emission in a grey atmosphere increase linearly with optical depth from TOA to the surface. The upwelling irradiance at the surface is F + (˜∗ ) = πB(Ts ) where Ts [K] is the surface skin τ temperature. The atmospheric temperature just above the surface is given by (131) as B(˜∗ ) = τ F0 (˜∗ + 1). Thus there is a temperature discontinuity between the near-surface air and the ground. τ 2π Climate models typically express h (124) in terms of the flux gradient with respect to pressure by invoking the hydrostatic equilibrium condition (456) B(˜) = τ h≡ dT dt 1 dF ρcp [−(ρg)−1 dp] g dF = − cp dp − + − + g (Fk − Fk ) − (Fk+1 − Fk+1 ) = − cp pk − pk+1 − − + + g (Fk − Fk+1 ) + (Fk+1 − Fk ) = cp pk+1 − pk = (132) (133) where the subscript denotes the kth vertical interface level in the atmosphere. 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. ˆ ˆ cos Θ = Ω · Ω (134) 2 RADIATIVE TRANSFER EQUATION 38 The case where incident and emergent directions are equal, i.e., Ω = Ω corresponds to Θ = 0. When the scattered direction continues moving in the forward hemisphere (relative to the plane ˆ defined by Ω ), it is called forward scattering, and corresponds to Θ < π/2. When scattered radiation has been reflected back into the hemisphere from whence it arrived, it is called back scattering, and corresponds to Θ > π/2. The case where incident and emergent directions are ˆ ˆ opposite, i.e., Ω = −Ω, corresponds to Θ = π. The Cartesian components of Ω and Ω are straigtforward to obtain in spherical polar coordinates. ˆ ˆ Ω = sin θ cos φ ˆ + sin θ sin φ ˆ + cos θ k i j ˆ ˆ Ω = sin θ cos φ ˆ + sin θ sin φ ˆ + cos θ k i j (135) (136) ˆ ˆ The scattering angle Θ is simply related to the inner product of Ω and Ω by the cosine law cos Θ = = = = 2.2.7 ˆ ˆ Ω ·Ω sin θ cos φ sin θ cos φ + sin θ sin φ sin θ sin φ + cos θ + cos θ sin θ sin θ (cos φ cos φ + sin φ sin φ) + cos θ cos θ sin θ sin θ cos(φ − φ) + cos θ cos θ (137) 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 1 4π 1 4π φ=2π φ=0 θ=π θ=0 p(cos Θ) dΩ = 1 Ω (138a) (138b) p(θ , φ ; θ, φ) sin θ dθ dφ = 1 The dimensions of the phase function are somewhat ambiguous. If the (4π)−1 factor in (138) is assumed to be steradians, then p is a true probability and is dimensionless. However, if the (4π)−1 factor is considered to be numeric and dimensionless (i.e., a probability), then p has units of (dimensionless) probability per (dimensional) steradian, sr−1 . The latter convention best expresses the physical meaning of the phase function and is adopted in this text. It is therefore important to remember that factors of (4π)−1 which multiply the scattering integral in the radiative transfer equation are considered to be dimensionless in the formulations which follow, e.g., (152). Furthermore, the units of p are probability per steradian, sr−1 . ˆ ˆ 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 Θ. 2 RADIATIVE TRANSFER EQUATION 39 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 1 2 θ=π p(θ ; θ) sin θ dθ = 1 θ=0 (139) In accord with the discussion above, the factor of 1/2 is dimensionless, as is the RHS. 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 n=N p(cos Θ) = n=0 (2n + 1)χn Pn (cos Θ) (140) An expansion of order N contains N + 1 terms. The zeroth order polynomial and coefficient are identically 1. 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: 1 δn,k 2n + 1 u=−1 1 u=1 1 Pn (cos Θ)Pk (cos Θ) d(cos Θ) = δn,k 2 u=−1 2n + 1 Pn (u)Pk (u) du = 1 2 u=1 (141a) (141b) Some other properties of Legendre polynomials are discussed in §10.5. The expansion coefficients χn are defined by the projection of the corresponding Legendre polynomial onto the phase function. χn = 1 2 1 = 2 1 = 2 u=1 Pn (u)p(u) du u=−1 Θ=0 (142a) (142b) (142c) Pn (cos Θ)p(cos Θ) d(cos Θ) Θ=π Θ=π Pn (cos Θ)p(cos Θ) sin Θ dΘ Θ=0 map the phase function into a series of Legendre polynomials. The χn are also called the phase function moments. Radiative transfer programs like DISORT [Stamnes et al., 1988] usually require that directional information be specified as a Legendre polynomial expansion. 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 2 RADIATIVE TRANSFER EQUATION 40 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: χn = 1 2 l=L Hl Pn (cos Θl )p(cos Θl ) sin Θl l=1 (143) where L is the number of Gauss-Lobatto abscissae. Lobatto quadrature is discussed more thorought in § 10.6. 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. p(cos Θ) = 1 − g2 (1 + g 2 − 2g cos Θ)3/2 (144) The n’th moment in the Legendre expansion of (144) is, conveniently, g n . χn = g n Applying (145) in (146) gives n=N (145) p(cos Θ) = n=0 (2n + 1)g n Pn (cos Θ) (146) This convenient property (145) makes numerical quadrature of the Legendre expansion coefficients unnecessary once the first moment, g = χ1 , is known. 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, Is (τ, Ω), and a diffuse component, Id (τ, Ω) such that − − − ˆ ˆ ˆ Iν (τ, Ω) = Is (τ, Ω) + Id (τ, Ω) (147) where we have suppressed the ν subscript on the RHS to simplify notation. The direct component refers to any photons contributing from a collimated source which have not (yet) been scattered. Typically the collimated source is solar radiation, and so we subscript the direct component with s for “solar”. There may be a corresponding direct component of intensity in the upwelling direction 2 RADIATIVE TRANSFER EQUATION 41 + Is if the reflectance at the lower surface is specular, e.g., ocean glint. In such situations it is straightforward to define + + + ˆ ˆ ˆ Iν (τ, Ω) = Is (τ, Ω) + Id (τ, Ω) (148) Of course there exist planetary atmospheres somewhere which are illuminated by multiple stars. We shall neglect upwelling solar beams for the time being. For consistency, though, we shall shall + − use Id rather than I + in equations in which Id also appears. − ˆ It is plain that Is (τ, Ω) 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), − ˆ ˆ ˆ Is (τ, Ω) = F e−τ /µ0 δ(Ω − Ω ) = F e−τ /µ0 δ(µ − µ0 )δ(φ − φ0 ) (149) 2.2.11 Source Function t The source function for thermal emission Sν is t Sν = j(ν) k(ν) α(ν) = Bν (T ) k(ν) k(ν) − α(ν) Bν (T ) = k(ν) = (1 − )Bν (T ) (150) ˆ ˆ With knowledge of the phase function p(Ω , Ω) as well as the scattering and absorption coefficients of the particular extinction process, we can determine the contribution of scattering to the total ˆ source function Sν . Consider the change in radiance dIν (Ω) occuring over a small change in path ˆ ˆ ˆ d˜ due to a scattering process with phase function p(Ω , Ω). The scattering contribution to dIν (Ω) τ ˆ ˆ from every incident direction Ω is proportional the radiance of the incident beam, Iν (Ω ), the ˆ ) at location τ is due to scattering (not to absorption), and to the probability that extinction of Iν (Ω ˜ ˆ ˆ normalized phase function p(Ω , Ω)/4π (138). ˆ ˆ ˆ dIν (Ω) σ(ν) p(Ω , Ω) ˆ = Iν (Ω ) d˜ τ σ(ν) + α(ν) 4π ˆ ˆ ˆ = Iν (Ω )p(Ω , Ω) 4π (151) s The source function for scattering Sν is obtained by integrating (151) over all possible incident ˆ ˆ directions Ω that contribute the exiting radiance in direction Ω s Sν = 4π Iν (˜, Ω )p(˜, Ω , Ω) dΩ τ ˆ τ ˆ ˆ Ω (152) 2 RADIATIVE TRANSFER EQUATION 42 The LHS and RHS of (152) must have equal dimensions. This is easily verified: is dimen−1 s sionless, (4π) is dimensionless (cf. §2.2.7), Iν and Sν both have dimensions of intensity, the dimensions of p are sr−1 and the cancel the dimensions of dΩ which are sr. 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 , pi ) and sum their contributions to obtain s s to total scattering source function Sν = i=n Sν,i . We neglect such details for now and merely i=1 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) t s Sν = Sν + Sν Sν = (1 − )Bν (T ) + 4π Iν (˜, Ω )p(˜, Ω , Ω) dΩ τ ˆ τ ˆ ˆ Ω (153) Although (153) is complete for sources within a medium, additional terms may need to be added to account for boundary sources, such as surface reflection. Boundary sources are discussed in §2.3. 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 dIν d˜ τ = −Iν + Sν = −Iν + (1 − )Bν (T ) + 4π Iν (˜, Ω )p(˜, Ω , Ω) dΩ τ ˆ τ ˆ ˆ Ω (154) This is a general form for the equation of radiative transfer in one dimension, and the jumping off point for our discussion of various solution techniques. The physics of thermal emission, absorption, and scattering are all embodied in (154). The most appropriate solution technique for (154) depends on the boundary conditions of the particular problem, the fields required in the solution (e.g., if irradiance is required but radiance is not), and the required accuracy of the solution. 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 −µ ˆ dI − (τ, Ω) ˆ = I − (τ, Ω) − (1 − dτ − µ 4π − )B − 4π ˆ ˆ ˆ I + (τ, Ω )p(+Ω , −Ω) dΩ + ˆ ˆ ˆ I − (τ, Ω )p(−Ω , −Ω) dΩ )B − ˆ ˆ ˆ I + (τ, Ω )p(+Ω , +Ω) dΩ + (155a) ˆ dI + (τ, Ω) ˆ = I + (τ, Ω) − (1 − dτ − 4π − 4π ˆ ˆ ˆ I − (τ, Ω )p(−Ω , +Ω) dΩ (155b) 2 RADIATIVE TRANSFER EQUATION 43 If we substitute (147) into (155a) we obtain −µ − ˆ ˆ dId (τ, Ω) dI − (τ, Ω) −µ s = dτ dτ − ˆ ˆ I − (τ, Ω) + Is (τ, Ω) − (1 − d )B − 4π − ˆ ˆ ˆ Is (τ, Ω )p(−Ω , −Ω) dΩ − − ˆ ˆ ˆ Id (τ, Ω )p(−Ω , −Ω) dΩ − 4π + ˆ ˆ ˆ Id (τ, Ω )p(+Ω , −Ω) dΩ − + 4π − − − The direct beam (150) satisfies the extinction law (91) so that −µ dIs /dτ = Is . Thus the second terms on the LHS and the RHS of (156) cancel eachother and we are left with − ˆ dId (τ, Ω) − ˆ ˆ = Id (τ, Ω) − (1 − )B − S ∗ (τ, −Ω) dτ ˆ ˆ ˆ ˆ ˆ ˆ I + (τ, Ω )p(+Ω , −Ω) dΩ − I − (τ, Ω )p(−Ω , −Ω) dΩ − 4π + d 4π − d −µ (156) where S ∗ is called the single-scattering source function. S ∗ is defined by the integral containing the direct beam on the RHS of (156) ˆ S ∗ (τ, −Ω) = = = 4π 4π 4π − ˆ ˆ ˆ Is (τ, Ω )p(−Ω , −Ω) dΩ − ˆ ˆ ˆ ˆ F e−τ /µ0 δ(Ω , Ω )p(−Ω , −Ω) dΩ − ˆ ˆ F e−τ /µ0 p(−Ω , −Ω) (157) Note that (156) is an equation for the diffuse downwelling radiance, not for the total downwelling radiance. The relation between the diffuse and total radiance is given by (147). The direct component is always known (150) once the optical depth has been determined. Likewise, inserting (147) and (150) into (155b), we obtain the equation for the diffuse upwelling intensity + ˆ dId (τ, Ω) + ˆ ˆ = Id (τ, Ω) − (1 − )B − S ∗ (τ, +Ω) µ dτ ˆ ˆ ˆ ˆ ˆ ˆ − I + (τ, Ω )p(+Ω , +Ω) dΩ − I − (τ, Ω )p(−Ω , +Ω) dΩ 4π + d 4π − d (158) where ˆ S ∗ (τ, +Ω) = = = 4π 4π 4π − ˆ ˆ ˆ Is (τ, Ω )p(−Ω , +Ω) dΩ − ˆ ˆ ˆ ˆ F e−τ /µ0 δ(Ω , Ω )p(−Ω , +Ω) dΩ − ˆ ˆ F e−τ /µ0 p(−Ω , +Ω) (159) 2 RADIATIVE TRANSFER EQUATION Azimuthal Mean Radiation Field 44 2.2.13 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 ¯ Xφ = 1 2π 2π X(φ) dφ 0 (160) For future reference we present also the polar angle integration operator which will be used to formulate the two stream equations. ¯ Xθ = 0 1 π X(θ) sin θ dθ X(µ) dµ 0 ¯ Xµ = (161) Applying (160) to (25) and (138), we obtain the azimuthal mean hemispheric intensity, phase function, and single-scattering source function respectively, 1 I (τ, µ) = 2π 1 p(±µ ; ±µ) = 2π ± 2π I ± (τ, µ, φ) dφ 0 2π (162) (163) (164) p(±µ , φ ; ±µ, φ) dφ 0 S ∗ (τ, ±µ) = 4π F e−τ /µ0 p(−µ0 , ±µ) No additional symbols are used to indicate azimuthal mean quantities. The presence of only the polar angle (θ or µ) on the LHS of (162) indicates that azimuthal mean quantities are involved. 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 −µ − − dId (τ, µ) − = Id (τ, µ) − (1 − dτ )B − 2 4π − F e−τ /µ0 p(−µ0 , −µ) (165a) 2 + dI + (τ, µ) + µ d = Id (τ, µ) − (1 − dτ − 2 + + Id (τ, µ )p(+µ , −µ) dµ − − Id (τ, µ )p(−µ , −µ) dµ )B − 2 4π − F e−τ /µ0 p(−µ0 , +µ) − Id (τ, µ )p(−µ , +µ) dµ + Id (τ, µ )p(+µ , +µ) dµ − (165b) The azimuthal mean equations (165) are very similar to the full equations (156) and (158). While the remainder of the terms in (165) contain only one azimuthally dependent intensity, the scattering integrals contain two, p and I ± . This causes the factor of /4π in front of the scattering integrals has to become a factor of /2. 2 RADIATIVE TRANSFER EQUATION 45 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. 1 − 0 dI − (τ, µ) dµ = µ d dτ 1 1 − Id (τ, µ) dµ 0 1 1 − (1 − ) 0 1 B dµ − 4π F e −τ /µ0 0 p(−µ0 , −µ) dµ − 2 + Id (τ, µ )p(+µ , −µ) dµ dµ − 0 + 2 0 − − Id (τ, µ )p(−µ , −µ) dµ dµ (166) 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 µ. ¯ 1 µ 0 − dId (τ, µ) d dµ ≈ µ ¯ dτ dτ 1 − Id (τ, µ) dµ = µ ¯ 0 − dId (τ ) dτ (167) Extracting µ from the integral on the LHS of (167) is an approximation whose validity worsens as ¯ ± the correlation between µ and Id (τ, µ) increases. The first two terms on the RHS of (166) are simple hemispheric integrals of hemisphericallyvarying 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 1 vanishes since 0 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 µ . 1 + p(+µ , −µ)Id (τ, µ ) dµ 0 + 1 ≈ 0 + Id (τ, µ ) dµ + p(+µ , −µ) dµ (168) + = Id (τ ) + p(+µ , −µ) dµ These terms define the backscattered fraction of the radiation field. Equation (168), for example, represents the fraction of upwelling energy backscattered into the downwelling direction. This is one of four similar terms that appear in the coupled equations (166). We define the azimuthal mean backscattering function b(µ) as b(µ) ≡ 1 2 1 p(−µ , µ) dµ = 0 1 2 1 p(µ , −µ) dµ 0 (169) The complement of the backscattering function is the forward scattering function f . f (µ) = 1 − b(µ) ≡ 1 2 1 p(−µ , −µ) dµ = 0 1 2 1 p(µ , µ) dµ 0 (170) 2 RADIATIVE TRANSFER EQUATION 46 Pairs of the terms in (169)–(170) are identical due to reciprocity relations. Reciprocity relations state that photon paths are reversible. The hemispheric mean backscattering coefficient b is the hemispheric integral of (169) 1 b ≡ 0 b(µ) dµ 1 2 1 0 0 1 (171) (172) = p(−µ , µ) dµ dµ The hemispheric mean forward scattering coefficient is defined analogously. The result of these two steps is −¯ µ − − dId (τ ) − = Id (τ ) − (1 − dτ )B − 4π F e−τ /µ0 p(−µ0 , −µ) (173a) I + (τ )p(+µ , −µ) dµ − I − (τ )p(−µ , −µ) dµ 2 + d 2 − d + dId (τ ) + = Id (τ ) − (1 − )B − F e−τ /µ0 p(−µ0 , +µ) µ ¯ dτ 4π − 2 + Id (τ )p(+µ , +µ) dµ − + 2 − − Id (τ )p(−µ , +µ) dµ (173b) These are the hemispheric mean, azimuthal mean, two stream equations for the radiation field in slab geometry. It should be clear that these equations (173) are not derived simply by performing a hemispheric integral (161) on (165). A strict hemispheric averaging operation would applies to ± the product of the phase function and and the intensity, not to each separately, so that extracting Id from the scattering integrals is an approximation, as is the definition of µ. ¯ 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 + dIν (τ ) + + − ∗+ = Iν (τ ) − (1 − b)Iν (τ ) − bIν (τ ) − (1 − )Bν − Sν (174a) µ ¯ dτ dI − (τ ) + ∗− − − −¯ ν µ (174b) = Iν (τ ) − (1 − b)Iν (τ ) − bIν (τ ) − (1 − )Bν − Sν dτ For the remainder of the derivation we drop the ν subscript and the explicit dependence of I ± on τ . Adding and subtracting we obtain d(I + + I − ) = −(α − β)(I + − I − ) dτ d(I + − I − ) = −(α + β)(I + + I − ) dτ where we have defined α ≡ −[1 − β ≡ b/¯ µ (1 − b)]/¯ µ (175a) (175b) (176a) (176b) 2 RADIATIVE TRANSFER EQUATION Diffusivity Approximation 47 2.2.15 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) 1 + F =2 π ν 1 − F =2 π ν 1 Iν (+µ)µ dµ 0 1 (177a) (177b) Iν (−µ)µ dµ 0 The factor of π has been moved to the LHS so that, when the source function is thermal, the RHS is the integral Planck function (116). These angular integrals may be reduced to a calculation of the intensity at certain quadrature points (zenith angles) in each hemisphere with surprising accuracy. The location of the optimal quadrature points may be arrived at through both theoretical and empirical methods. Two point full Gaussian quadrature (§10.5) tells us the optimal angles to evaluate (177) at are u = ±3−1/2 , θ = ±54.7356◦ . 1 + F π ν 1 − F π ν 2 √ Iν (+3−1/2 ) 3 2 √ Iν (−3−1/2 ) 3 (178a) (178b) The relation in Equation (178) is exact when Iν (µ) is a polynomial of degree < 2. For two point Gaussian quadrature the quadrature weights happen to be unity, but not so for three point (or higher order) Gaussian quadrature. Here is the four point quadrature version of (178) u1 = −u4 u2 = −u3 A1 = A4 A2 = A3 1 + F 2π ν 1 − F 2π ν = = = = 0.339981 0.861136 0.652145 0.347855 A1 u1 Iν (u1 ) + A2 u2 Iν (u2 ) A3 u3 Iν (u3 ) + A4 u4 Iν (u4 ) (179) 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 1 + F Iν (u = +D−1 ) (180a) π ν 1 − Fν Iν (u = −D−1 ) (180b) π Comparing (180) to (178) and (177) shows that D replaces, simultaneously, the quadrature angle, its weight, the factor of µ in (177), and a factor of 2 from the azimuthal integration. The angle 2 RADIATIVE TRANSFER EQUATION 48 θD = arccos D−1 may be interpreted as the mean slant path of radiation in an isotropic, nonscattering atmosphere. 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 T from point P1 to point P2 is the likelihood a photon traveling in ˆ at P1 will arrive at P2 without having interacted with the matter in between. direction Ω ˜ T (P1 , P2 ) = exp[−˜(P1 , P2 )] τ (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 T (z , z; µ) = exp[−τ (z , z)/µ] (182) Thus, in common with absorptance A and reflectance R, the transmittance T ∈ [0, 1]9 The vertical gradient of T is frequently used to formulate solutions of the radiative transfer equation in non-scattering, thermal atmospheres. ∂ ∂ T (z , z; µ) = exp[−τ (z , z; µ)/µ] ∂z ∂z ∂T (z , z; µ) 1 ∂τ (z , z; µ) = − exp[−τ (z , z; µ)/µ] ∂z µ ∂z 1 ∂τ = − T (z , z; µ) µ ∂z Inverting the above, we obtain µ−1 T (z , z; µ) = − ∂T (z , z; µ) ∂z ∂z ∂τ (184) (183) 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 text distinguishes between transmittance and transmission. The former measures the potential for the latter to occur. Thus a transmittance of unity means the transmission of energy is very high. A similar rule holds for reflectance and reflection and for absorptance and absorption. 9 2 RADIATIVE TRANSFER EQUATION 49 this we use z(τ = τ ∗ ) z(τ = 0) T (0, z; µ) T (z , z; µ) T (z, z ; µ) ∂T (z , z; µ) ∂z = = = = = 0 ∞ ∗ e−(τ −τ )/µ e−(τ −τ )/µ e−(τ −τ )/µ ∂T (z, z ; µ) = − ∂z (185) The last relation simply states that the transmission decreases as the distance between z , z increases and visa versa. Substituting the above into (109) and (110) yields 0 Iν (z, +µ) = T (0, z; µ)Bν (0) + µ z −1 z T (z , z; µ)Bν (z ) (−dz ) Iν (z, −µ) = µ−1 ∞ T (z, z ; µ)Bν (z ) (−dz ) Rearranging terms we obtain z Iν (z, +µ) = T (0, z; µ)Bν (0) + µ−1 ∞ 0 T (z , z; µ)Bν (z ) dz (186a) (186b) Iν (z, −µ) = µ−1 z T (z, z ; µ)Bν (z ) dz Substituting (184) in the above leads to z Iν (z, +µ) = T (0, z; µ)Bν (0) + 0 ∞ ∂T (z , z; µ) Bν (z ) dz ∂z (187a) (187b) Iν (z, −µ) = − z ∂T (z, z ; µ) Bν (z ) dz ∂z If T is known, then Equations (187a)–(187b) are suitable for use in quadrature formulae to obtain irradiances. 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 A pertinent to fluxes. We shall introduce a painful but necessary menagerie of terminology to describe the various species of R, T , and A. 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 2 RADIATIVE TRANSFER EQUATION 50 sfc by providing more accurate net surface shortwave radiation estimates. This term, FSW , 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 ˆ ˆ ρ(ν, −Ω , Ω) ≡ + ˆ dIνr (Ω) − ˆ Iν (Ω ) cos θ dΩ (188) 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 Ω. + ˆ Iνr (Ω) = ˆ −Ω ˆ −Ω + ˆ dIνr (Ω) − ˆ ˆ ˆ Iν (Ω )ρ(ν, −Ω , Ω) cos θ dΩ = (189) ˆ ˆ The dependence of ρ(ν, −Ω , Ω) on two directions and on frequency makes it a difficult property to measure. 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 ˆ ˆ ρ(ν, −Ω , Ω) = ρL (ν) (190) The intensity reflected from a Lambertian surface is given by inserting (190) into (188) and then using (56b) + ˆ Iνr (Ω) = ˆ −Ω − ˆ Iν (Ω )ρL (ν) cos θ dΩ − ˆ Iν (Ω ) cos θ dΩ = ρL (ν) ˆ −Ω − = ρL (ν)Fν (191) As expected, the intensity reflected from a Lambertian surface depends only on the incident irradiance, and not at all on the details of the angular distribution of the incident intensity field. 2 RADIATIVE TRANSFER EQUATION 51 + The irradiance reflected from a Lambertian surface Fνr is the cosine-weighted integral of the reflected intensity (191) over all reflected angles + Fνr = + ˆ Iν (Ω) cos θ dΩ − ρL (ν)Fν cos θ dΩ ˆ Ω ˆ Ω = − ˆ Neither Fν and ρL (ν) depend on the emergent angle Ω so the integral reduces to the familiar integral of cos θ over the hemisphere (59) which is π. + − Fνr = πρL (ν)Fν (192) The reflected irradiance may not exceed the incident irradiance or the requirement of energy conservation will be violated. Therefore (192) shows that ρL (ν) ≤ π−1 (193) with equality holding only for a perfectly reflective (non-absorbing) Lambertian surface. 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 + ˆ Iνr (Ω) = ˆ ˆ ˆ ˆ F δ(Ω − Ω )ρ(ν, −Ω , Ω) cos θ dΩ ˆ −Ω ˆ ˆ = F ρ(ν, −Ω , Ω) cos θ0 (194) Integrating (194) over all reflected angles we obtain the diffusely reflected hemispheric irradiance + Fνr = ˆ +Ω ˆ ˆ F ρ(ν, −Ω , Ω) cos θ0 cos θ dΩ ˆ +Ω ˆ ˆ ρ(ν, −Ω , Ω) cos θ dΩ (195) = F cos θ0 The flux reflectance or plane albedo is the ratio of the reflected (diffuse) irradiance (195) to the incident (collimated) irradiance (150) + Fνr ˆ ρ(ν, −Ω , 2π) = − Fν = ˆ +Ω ˆ ˆ ρ(ν, −Ω , Ω) cos θ dΩ (196) The 2π indicates that the plane albedo pertains to the entire reflected hemisphere. 2 RADIATIVE TRANSFER EQUATION Table 4: Surface Albedoa Surface Type Reflectance Glacier 0.9 Snow Sea-ice Clouds Desert Savannah Forest Ocean Planetary a 52 0.8 0.6 0.2–0.7 0.2–0.3 0.2–0.25 0.05–0.1 0.05–0.1 0.3 Sources: Bird Entrails (2000) 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 + ˆ Fνr (Ω) = ˆ −Ω − ˆ ˆ Fν (−Ω )ρ(ν, −2π, Ω) cos θ dΩ (197) The flux reflectance of a Lambertian surface to collimated light is found by subsituting (190) in (196) ˆ ρ(ν, −Ω , 2π) = ˆ +Ω ρL