ABSTRACT ENVIRONMENTAL EFFECTS OF AIRCRAFT CONDENSATION TRAILS by David Randall Lyzenga Co-Chairmen: Ernest G. Fontheim, Robert H. Williams Condensation trails form as the result of the cooling of aircraft exhaust gases, which contain a high proportion of water vapor. The water vapor first liquifies and then freezes to form small spherical ice particles having radii between 1 and 10 microns. If the surrounding air is saturated with respect to ice, the contrails may spread to form an artificial cirrus-like cloud layer. The main environmental effects of such contrail layers, aside from their visual impact, are on the solar and terrestrial radiation fields in the atmosphere. These effects are evaluated for a layer of ice particles of radius r = 1, 3, and 10 microns, with an optical thickness T* = 0.2. The transmissivity of such a layer for terrestrial radiation varies from 0.97 to 0.73, depending on the particle size. Previous measurements on contrail layers indicate that the most likely value is 0.89, corresponding to a particle radius of 3 microns. The solar albedo of the layer is independent of particle size, but varies strongly with the solar zenith angle. The daily-average albedo is about five percent at low and mid

latitudes, increasing to 20 percent or more at high latitudes. Because of the low transmissivity of the atmosphere for long wave radiation, the effects at the surface of the earth are due mostly to the reflection of solar radiation by the layer. The surface effects are therefore relatively insensitive to the particle size, and are highly dependent on latitude. The reduction in the net radiative energy absorbed at the surface is on the order of.016 cal cm-2min-1 at midlatitudes for a continuous layer. Assuming radiative equilibrium of the surface, this would lead to a lowering of the of the surface temperature by about 20C. In areas of heavy air traffic the average contrail cover is probably about five percent, and the cooling effect is on the order of 0.10C. The largest potential effects occur in arctic regions, where the net radiation at the surface would be reduced by more than.02 c a.l cri-2Inirli-1 duringl tihc:.;;uirner pby;. ((,)lttI.lollo'; layer. Assuming that the only effect of this change is a reduction in the melting rate of ice, the total area of the arctic ice pack could be increased by 0.5 percent per year by a one percent contrail cover. This effect would be cumulative, and is significant in view of the possible role of the arctic ice cover as a positive feedback mechanism for climatic change.

T H E U N I VE R S IT Y O F MI C H I G A N COLLEGE OF ENGINEERING Department of Electrical & Computer Engineering Space Physics Research Laboratory ENVIRONMENTAL EFFECTS OF AIRCRAFT CONDENSATION TRAILS David Randall Lyzenga under contract with: National Aeronautics and Space Administration Grant No. NGR-23-005-015 Washington, D. C. November 1973

ACKNOWLE DGMENTS I am indebted to my doctoral committee for advice and encouragement throughout the course of this research. Without their support, and particularly that of my co-chairmen Dr. Ernest G. Fontheim and Assistant Professor Robert H. Williams, this work would never have been brought to completion. I would also like to thank Professor Julius London of The University of Colorado, Drs. Ralph J. Cicerone and Donald S. Stedman of The University of Michigan, and Mr. Richard J. Kauth of The Environmental Research Institute of Michigan for discussions on the subject matter of this thesis. Finally, I wish to express my appreciation to Mr. George R. Carignan, Director of the Space Physics Research Laboratory, for the use of the Laboratory's facilities, and to Roslyn A. Drummond for her careful typing of the manuscript. This work was carried out under a fellowship from the Institute for Environmental Quality of The University of Michigan. Additional support was provided by NASA Grant NGR-23-005-015. ii

TABLE OF CONTENTS ACKNOWLEDGMENTS........................................ ii LIST OF TABLES........................................ iv LIST OF FIGURES. v LIST OF APPENDICES..................................... vii CHAPTER 1: INTRODUCTION.............................. 1 1.1. Background 1 1.2. Previous Work on Climatic Effects of Contrails 2 1.3. Statement of the Problem 5 CHAPTER 2: PHYSICAL PROPERTIES OF CONTRAILS........ 7 2.1. Formation and Microstructure of Contrails 7 2.2. Single-Particle Scattering Properties 13 CHAPTER 3: RADIATIVE PROPERTIES OF CONTRAILS.......... 21 3.1. Radiation Fields in the Atmosphere 21 3.2. The Transfer Equation 23 3.3. Solutions for Diffuse Incident Radiation 27 3.4. Infrared Properties of Contrail Layers 39 3.5. Solutions for Direct Incident Radiation 44 3.6. Solar Albedo- of Contrail Layers 50 CHAPTER 4: CLIMATIC EFFECTS OF CONTRAILS.............. 58 4.1. Radiation Balance of the Surface and Atmosphere 58 4.2. Surface Changes Induced by Contrails 67 CHAPTER 5: CONCLUSIONS................................ 73 5.1. Extent of Contrail Formation 73 5.2. Assessment of Climatic Impact 76 5.3. Suggestions for Further Research 78 APPENDICES................................. o...... 80 REFERENCES......................................... 92 iii

LIST OF TABLES T'ABLE 1. Averaged Single-P'article Scattering Parameters for Ice Spheres......................... 19 2. Weighting Factors For Ilfrared Trarismissivity and Emissivity................................... 41 3. Average Solar Albedo of' Contrail Layer (T* = 0.2). 56 iv

LIST OF FIGURES FIGURE 1. Saturation water vapor densities over water and ice (Van Valin, 1971)........................ 8 2. Ice fog particle size distribution measured by Huffman (1968)................................ 10 3. Extinction efficiency for spheres and transparent irregular particles (adapted from Hodkinson, 1963)................... 17 4. Geometry of scattering layer with optical thickness T1, showing diffuse radiation incident on the bottom of the layer (T = T1) and emerging from the top of the layer (TX = 0).. 24 5. Transmissivity versus optical thickness Tl for the case wo = 0 (perfect absorption). E3(T1) is the third-order exponential integral.......... 36 6. Reflectivity versus (l-g)Tl for the case wo = 1 (perfect scattering). g = w1/3wo is the asymmetry factor of the phase function, and Tl is the optical thickness............................ 37 7. Reflectivity versus single-scattering albedo Wo for a semi-infinite layer with isotropic scattering....................................... 38 8. Diffuse transmissivity and emissivity of a layer of ice particles versus reduced optical thickness T* = n'rr2t, calculated from the two-stream approximation........................................ 4 3 9. Local albedo versus the cosine of the angle of incidence 0o for conservative scattering (wo = 1), asymmetry factor g =.844 and optical thickness T1 =.5, 4, and 16............................... 51 10. Local solar albedo versus the cosine of the solar zenith angle Go for T* = 0.2, computed from Equation (3.54) with Qext = 2.0 and g = 0.85..... 53 v

LIST OF FIGURES (concluded) FIGURE 11. Average solar albedo of a contrail layer with T* = 0.2, for summer (6 = 18~) and winter (6 = -180)............................ 55 12(a). Effect of a continuous contrail layer on the spring radiation balance of the atmosphere (upper curve) and surface (lower curve)...................................... 6 3 12(b). Effect of a continuous contrail layer on the summer radiation balance of the atmosphere (upper curve) and surface (lower curve)........................................ 64 12(c). Effect of a continuous contrail layer on the fall radiation balance of the atmosphere (upper curve) and surface (lower curve)....................................... 65 12(d). Effect of a continuous contrail layer on the winter radiation balance of the atmosphere (upper curve) and surface (lower curve).-6......................... 6 13. Changes in surface temperature for a continuous contrail layer over dry land..... 70 14. Relative humidity required for jet aircraft contrail formation as a function of pressure and temperature of the environment (from Appleman, 1966)........................* 74 15. Latitudes of greater-than-50% contrail probability by season at 50 mb and at 10 km (from Appleman, 1966)...................... 75 vi

LIST OF APIPENDICES APPENDIX A. Freezing Time and Growth Rate of Ice Particles 80 B. Eddington and Two-Stream Approximations....... 82 C. Derivation of Albedo Equation................. 86 D. Temperature of Contrail Layer................. 90 vii

INTRODUCTION 1. 1. BACKGROUND Aircraft condensation trails may be considered as a form of atmospheric pollution, and their effects on the environment should be evaluated. Some evidence of contrail formation can be seen almost every day, and occasionally a few contrails appear to fill the entire sky with a cirruslike cloud layer. Contrails can be seen distinctly on several ERTS-A satellite pictures taken at an altitude of 500 miles. The main effects of contrails, aside from their visual impact, are on the fluxes of solar and terrestrial radiation in the atmosphere. These are important since they determine the surface temperature of the earth and are the energy source for atmospheric motions, for evaporation of water, and for melting snow and ice in the spring. The present study attempts to assess the influence of contrails on the energy budget of the earth's surface and the atmosphere. The radiative properties of a typical contrail layer are computed from Mie scattering and radiative transfer theory, and the effects of the layer on the fluxes of solar and terrestrial radiation are evaluated. Finally, the climatic implications are discussed, although tentatively, since there are positive and negative feedback 1

2 meclianisms in the atmiiosphere which are not yet well understood. 1.". PII,;VIOUS WORiK ()N CLIIATIC EFFECTS OF CONTRAILS The possibility of climate modification by jet contr-ails was first mentioned by Appleman (1966), who speculated that either a heating or a cooling of the earth's surface niight result from a persistent contrail layer, depending on the properties of the layer and the relative amounts of' incoming and outgoing radiation. Appleman also made certain conjectures as to the dependence of the effects on the latitude and the thickness of the layer, which will be discussed in Chapter 4. Reinking (1968) estimated that the solar radiation reaching the surface of the earth might be reduced by as much as 20 percent by a cirrus deck formed solely from the spreading of numerous contrails. lie concluded that contrails could signlificantly affect the radiation budget at the earth's surface, and further suggested that contrails might conceivably be used to prevent crop damage caused by excessive sunlight and heat during the summer months. Nicodemus and McQuigg (1969) approached the problem from a purely observational viewpoint, by attempting to correlate the daily temperature range AT measured near the surface of the earth with the fractional reduction in solar radiation AS due to cirrus clouds. The relationship they obtained, for a sample of 54 summer days in central Missouri

was AT = 25.6 - 12.2 AS + e OF (1.1) where e is a random variable with a mean of zero and a standard deviation of 4.13 OF. They concluded that intentional generation of "contrail cirrus" during the day could be used to reduce daily maximum temperatures during the summer, but that this would not significantly affect the demand for irrigation water on agricultural land. It might be pointed out that this conclusion assumes that "contrail cirrus" has the same radiative properties as natural cirrus, and that the nighttime minimum temperature is not affected by the presence of the clouds. The most thorough experimental study of contrail radiative properties to date was made by Kuhn (1970). Using the NASA Convair 990 jet laboratory, he made measurements of the flux of sol]ar radiation and the intensity of infrared radiatiori inr tiLv' J_0-l2 ]i-ic(-'ron wvleIVell't~J,]t irlter'vatl. ab()OVCu llid below several contrail layers, which had an average thickness of 500 meters. From the solar radiation measurements, Kuhn deduced a solar albedo of 15% for the average contrail layer. From the infrared measurements he deduced a transmissivity of about 90% for radiation normal to the layer, and a reflectivity of less than 1%. He then constructed a model consistent with his infrared measurements. Assuming the layer to be made up of ice particles with an equivalent sphere

radius of 50 microns, and using Mie theory to calculate their scattering and absorbing properties, he found that a particle -3 concentration of.027 cm would give reasonable agreement with his measureirnerits In the 10-12 micron spectral interval. Onrc difficulty of' this model is that it predicts a much sinaller value for thie solar albedo than Kuhn claims to have mencasured. The agreement is improved by assuming a slmaller particle size, as will be shown in Chapter 3. Using this model to calculate the properties of the layer over the entire terrestrial spectrum, Kuhn then computed the change in the net radiation at the surface of the earth due to the presence of the contrail layer. The result of this computation was a 7% decrease in the total radiant power absorbed at the surface at a latitude of 13~N, and an 8% decrease at 390N. Assuming the surface to be in radiative equilibrium (which it is not, since there are convective heat exchanges between the surface and the atmosphere) Kuhn concluded that over a long period of tinie this would lead to a decrease in the surface temperature of about 6~C multiplied by the fraction of the time that such a contrail layer exists. As these calculations werle made for clear sky conditions, Kuhn arbitrarily reduced the results by one-half to account for tile masking effect of lower clouds. Thus his final result is, for a contrail persistence of 5%, a reduction in the surface temperature of 0.15~0C.

1.3. STATEMENT OF THE PROBLEM The present treatment of the contrail problem dififers from Kuhn's study in two respects. First, the solar and terrestrial radiative properties of the layer are both calculated from a physical model, which is based upon groundlevel observations of contrail microstructure. Since these observations indicate that the ice particle size varies, depending on atmospheric conditions, parallel calculations are carried out for three values of the particle radius: 1, 3, and 10 microns. Second, these radiative properties are used to calculate the change in the radiation balance, or the net radiation at the top of the atmosphere, as well as at the surface of the earth. The calculations are carried out at ten degree latitude intervals in the northern hemisphere, and both the seasonal variation and the latitude dependence of the effects are examined. The organization of this report is as follows: Chapter 2 contains a discussion of the formation and physical properties of contrails, as far as these are known. The single-scattering properties of the constituent particles are then derived. In Chapter 3 the radiative properties of a typical contrail cirrus layer are calculated, by means of the radiative transfer equation. In Chapter 4 these results are applied to the radiation fluxes in the earth's atmosphere,

6 to determine tile effcctlsn of a contrail layer on the radiation budget of the atmosphere and the surface. Conclusions anrid recommnendations for further research are presented in Chapter 5. This research was partly inspired by the reports of tile Study of Critical Environmental Problems (1970) and the St udy of' Man's Impact on Climate (1971). Both of these reports mentioned the possibility of climate modification by jet contrails, and made recommendations for research in this area. Somne specific recommendations were that the radiative properties of representative contrails and contrail-produced cirrus clouds be determined, that the changes in cirrus cloud population be monitored, and that the significance, if any, of ice crystals falling from contrail clouds as a source of freezing nuclei for lower clouds be investigated. The present research is essentially a response to the first recommendation. The role of contrails in seeding lower clouds is tile subject of papers by Murcray (1970) and Dingle (1971), and will tiot be dealt with in this report.

CHAPTER 2 PHYSICAL PROPERTIES OF CONTRAILS 2.1. FORMATION AND MICROSTRUCTURE OF COiTRAILS The combustion of fuel in jet aircraft engines produces water vapor, which is expelled at high temperatures from the exhaust nozzle. As the exhaust expands and mixes with the ambient air, it cools until saturation is reached and the water begins to condense in the form of very small droplets. Since the ambient temperature is far below freezing (from -40~C to -60~C) at cruise altitudes, the droplets rapidly freeze, forming roughly spherical ice particles which continue to grow by condensation as long as the environment is supersaturated with respect to ice. Two factors cause the resulting size distribution to be fairly narrow. First, the freezing time for the droplets is inversely proportional to their volume (see Appendix A) so the larger droplets freeze before the smaller ones do. Since the saturation vapor pressure of ice is lower than that of liquid water (cf. Figure 1), the ice particles then grow rapidly in the ice-supersaturated environment while the remaining water droplets evaporate as the water vapor pressure falls below the liquid-saturation point. This mechanism effectively eliminates particle sizes below a certain limit. This was confirmed in a laboratory study of contrail formation 7

'(IL6T'UTTIA UeA) aOT puP aGnxM IGAO s9T;TsuGp JodeA LIa4eM uoTeaxnqes 1 aT nS T 30'aJnloaduwjal 001- 06- ot- 09- 00- O.- 0,- oa -01 CD (U): 0 is\ I I 6-01 -" \\S J-a:M! J aoo g-01 on (~sa) JaloMv JAO o UOJOA S -lt\ 9 t A _ _ i —_-_ * _ - -- | l9 _08

9 by Pili6 and Jiusto (1958), in which aircraft fuel was burnled in an ice-supersaturated environment at -520C and 400 mb pressure. The condensate was collected and photographed at 30 second intervals: during the first two intervals the droplets with diameters less then 1.5 microns disappeared altogether, while the larger particles increased in size by as much as 0.5 micron. The second factor which tends to narrow the size distribution is the dependence of the growth rate on the particle size. The rate of increase in the ice particle radius due to condensation is roughly inversely proportional to the radius (see Appendix A). Thus the smaller particles grow faster and tend to "catch-up" to the larger particles in time. Huffman (1968) measured the particle size distributions occurring in ground-level ice fogs produced by jet aircraft during the Alaskan winter. Iliis measurements showed a size distribution peaking at about 4 microns diameter (cf. Figure 2) for fogs produced by aircraft under a variety of atmospheric conditions. Ice fogs produced by lower-temperature sources (heating plants, open water) had a larger mode diameter, up to 16 microns, depending more sensitively on the air temperature. A theoretical model for the formation and growth of ice fog particles has been developed by Huffman and Ohtake (1971). Simultaneous equations for the saturation ratio, nucleation rate, particle growth rate, and freezing rate were numerically solved as a function of time to yield the final size distribution. According to this model, the size distribution is

40 Method: Impaction Location: Eielson Date: 14 Feb. 1967 T ime: 0030 30 Temp. ("C): -40 N = 1740 0 20 X I' 20 o, I I, _ 1 1,11,,.. L 0 2 4 6 8 10 12 14 16 DIAMETER,,U Figure 2. Ice fog particle size distribution measured by Huffman (1968).

11 determined by the temperature difference between the source and the environment (the greater the temperature difference, the smaller the particle sizne). Agreement between this theory and the obsc.rvattions is good for lower-temperature sources, but the model predicts too-small particle sizes for fogs produced by automobile exhaust. Jet aircraft exhaust is not considered in this paper, but the model would very likely break down for that case also, since the source temperature is higher for jet exhaust than for automobile exhaust. The authors speculated that the reason for the breakdown of the theory in the case of hightemperature sources is the incomplete freezing of the very small droplets produced by these sources. If all the particles freeze simultaneously, the frozen particles can grow only by condensation of the water vapor in the air. however, if a significant fraction of the droplets remain unfrozen, these will evaporate as soon as the water vapor pressure falls below the liquid-saturation point, thus supplying more water vapor for the growth of the frozen particles. The result is that the average diameter of the ice particles is larger than it would have been if the freezing had been complete. The problem of theoretically predicting the exact size distribution of ice particles in jet contrails has not yet been completely solved. However, the relative success of Huffman and Ohtake's model seems to justify the conclusion that contrail s formed at cruise altitude have roughly the same microstructure as the ice fogs observed by liuffman at

12 ground level. Pilie and Jiusto's work on laboratory contrails formed at reduced pressure seems to support this conclusion, as their photographs show particles of nearly the same size as those observed by lluffman. We lnay there Fore conclude that contrails are composed initially of roughly spherical ice particles of approximately 2 microns radius. Under the proper conditions (i.e., supersaturation with respect to ice), these particles may grow furthler as tlhe contrail spreads to form what has been called "contrail cirrus". Tlhe ice particles in contrail cirrus probably never grow as large as those found in natural cirrus clouds, since cirrus clouds are formed on longer time scales and with more available water vapor than contrails. Hansen and Pollack (1970) made calculations of near-infrared light scattering from cirrus clouds for various values of the particle radius, and found that the best agreement with observations was obtained for an equivalent particle radius of 10 micronrs, wiichi they claimed to be "consistent with values typically qbtained from direct sampling measurements of ice clouds." As a result of these calculations, it was decided to make parallel calculations of the contrail radiative properties at three values of the particle radius: r = 1, 3, and 10 microns. The first two values are probable limits for the particle sizes encountered in newly formed contrails, while the last value is taken to be the maximum size attained in fully developed contrail cirrus layers. The fact that

13 some particles may grow larger than 10 microlns is not crucial since the radiative properties for solar aild- ter restrial radiation do not charge greatly as the part ic le Iadius is increased above 10 microns. 2.2. SINGLE-PARTICLE SCATTERING PROPERTIES When electromagnetic radiation is incident on a particle of matter, it is partly absorbed and partly scattered, in a manner which depends on the wavelength of the radiation and the size, shape, and composition of the particle. The problem of calculating the scattering properties for spherical particles of arbitrary size was first solved by G. Mie in 1908. Essentially, the problem is solved by expanding the incident plane wave and the scattered wave (inside and outside of the sphere) in a series of spherical Bessel functions, and applying the proper boundary conditions at the surface of the sphere. The single-scattering properties of interest include the scattering cross section as and the absorption cross section a a. Multiplying these cross sections by the concentration or number density of particles gives the probability per unit path length that an incident photon will undergo scattering or absorption, respectively. The probability that a photon will be either scattered or absorbed is given by the total or extinction cross section at = as + a a multiplied by the particle con cent ration. The results of the Mie theory are usually written in terms of the corresponding efficiency factors Qsca, Qabs' and

1l Qext' These are equal to the cross sections for scattering, absorption, and extinction, respectively, divided by the geometrical cross section wrr2 (for particles of radius r). These quantities are expressed as infinite series involving the parameter x 2rrr/A, wlhere X is the wavelength of the radiation, as Follows (van de Hulst, 1957, and Irvine, 1965a): Qext =2 j (2nr + 1) Re(an + bri) et n= 1 Qsca 2= 2 (2n + l){IanI2 + Ibn 2} nl Qabs Qext Qsca (2.1) where''(mx) ~n(X) - min(mx) *n(X) n n mn n n- n m'n(mx) in(x) - 4n(mX) n'(x) b r n (2 2) n mi'(mx) n (x) - kn(mx) Cn(x) and il is the conmpllex index of refraction of thc A pai ti.i.. The functions in and rn (and their derivatives'n and C ) are def i led irn t~er mls oFL the spherical Bessel functions iJn and h n(2) as foJlows: d(z) - zj2(z) 1 (z) - zh (z). (2.3) n ~n

15 The angular distribution of the scattered radiation is described by the phase function P(O,P), where 0 aiid P are the polar and azimuthal angles of the scattered rtadiation, nleasured relative to the incident direction. The angle 0 between the incident and scattered radiation is called the scattering angle. For unpolarized radiation incident on spherical particles, the scattered radiation is independent of the azimuthal angle P, and the dependence on the scattering angle may be approximated by the following function (Henyey and Greenstein, 19 41): HG() = Wo(l - g2) (2.4) [1 - g2 + 2g cose]3/2 where wo is the single-scattering albedo, and g is the asymmetry factor of the phase function. These quantities are defined as follows: o = Qsca/Qext =;-ff (e)dQ g = (cos-> = f+-w f ()cosOdQ (2.5) The single-scattering albedo can be calculated from Equations (2.1), and the asymmetry factor is given by the following expression (Irvine, 1965a): g=E Re ni+2 b b* ) + jn+l.+ ~x~~2Q&~ n+l ~(nn+ n n+l n(n+l) n.n sca n=l (2.6)

16 1ti1 tile limit x<<l, the Mie solution reduces to the ]Rayleigh scattering case, in which the scattering efficiency i, proportional to x4 (or to -4 for fixed radius) and the absorptionl efficiency is proportional to x. This is the fanlilar explanatiorl of tie blue color of the sky, since the silortetL wave rlengt}is or tihe sunlitgllt are scattered much more strongllgy tilan the longer wavelergtiis s in the atmosphere. In the limit x>>l, the MlVie solution reduces to the geometric optics case, in which the scattering properties can be calculated by ray tracing, taking into account the effects of refraction, reflection, and diffraction at the surface of the sphere. The extinction efficiency in this case has the asymptotic value 2, of which one unit is due to direct interception of rays by the particle and one unit is due to the diffraction of rays passing near the edge of the particle. The asymmetry factor also reaches an asymptotic limit, but its value depends 4 on the index of refraction of the particle. For m = 3, which approximates the index of refraction of ice for visible wavelengths, the asymptotic value of g is about 0.85 (Irvine, 1965a). The behavior of the scattering properties in the transition region between these two limits is illustrated by the curve of Qext versus p = 2x(m-1) in Figure 3. For nonabsorbing, perfectly spherical particles the curve for Qext exhibits an oscillatory behavior due to resonance of the internally reflecatecd rays. If there is absorption or if the particles are rnot perfectly sphelrical, this oscillatory behavior is smoothed and the extirlction efficiency rises steadily from zero to its

17 MIE THEORY FOR SPHERES 0L3 IL W-J 2 0 Xr - MEASUREMENT ON QUARTZ PARTICLES!w I 0 5 10 15 20 25 p = 2x(m-I) Figure 3. Extinction efficiency for spheres and transparent, irregular particles (adapted from Hodkinson, 1963).

18 asymptotic value of 2. The asymmetry factor reaches its asymptotic value at about x = 4, according to calculations by Irvine (1965a). Thlle transition is smooth if the index of refraction is near unity, or if there i-; absorption within the particle, but thetre is some osciillation for nonabsorbing particles with n > The solar spectrum ranges, for practical purposes, from about 0.3 to 3.0 microns wavelength and peaks at about 0.55 microns. Over this range of wavelengths the index of 4 refraction of ice is approximately m = 3, and absorption may be neglected. For ice particles with radii between 1 and 10 microns, the scattering properties may be considered to have reached their asymptotic values over most of this wavelength interval, especially since the particles are not perfect spheres. Therefore, the following values were selected to describe the scattering of solar radiation by contrail ice particles: Qex 2.0 ext o = 1.0 g = 0.85 (2.7) The terrestrial radiation spectrum extends from roughly 4 to 100 microns and peaks at about. 10 microns wavelength. The scattering properties vary widely over this interval, so the spectrum was subdivided into the five bands indicated in Table 1, and the radiative properties were calculated

L9'TABLE 1 AVERAGED SINGLE-P ARITl'ICLE SCATTERING PARAMETER}S 1'OR ICE SPHERES Band 1 2 3 4 5 Wavelength 4-8 8-12 12-20 20-40 40-100 (microns) r=lp Qext.190.103.078.027.060 to.398.102.117.035.001 g.232.074.035.011.002 r= 3 QxtQ 1.40.636.784.142.196 e xt ow.709.474.639.362.024 g.806.632.360.096.021 r=10 Qext 2.73 2.59 3.37 1.86 1.06 o.701.688.729.756.314 g.909.916.787.687.274

20 separately for each band. The values for Qext' Wo, and g listed in Table 1 represent the results of Mie theory calculations by Irvine and Pollack (1968) for ice spheres, averaged over each of the bands indicated. Note that the extinction coefficient varies over more th;rn an order of mnagnitude with particle size, in contrast withl the extinction coefficient for solar radiation which has the same value in each case. For r = 1l, the scattering properties are in the Rayleigh scattering limit over most of the terrestrial spectrum, while for r = 10 the extinction coefficient is near its maximum resonance value. The relative magnitude of the effects on the solar and terrestrial radiation thus depends strongly on the particle size. Calculations of these effects are carried out separately for each particle size in the following chapter.

CHAPTER 3 RADIATIVE PROPERTIES OF CONTRAILS 3.1. RADIATION FIELDS IN THE ATMOSPHERE At every point in the earth's atmosphere there is radiation of both solar and terrestrial origin. These radiation fields differ in wavelength, as was mentioned in Chapter 2, but they also differ in their directional dependence. Terrestrial radiation is diffuse, consisting of photons traveling at all directions, whereas the solar radiation incident on the top of the atmosphere is in the form of a parallel beam. In the atmosphere part of the incident solar radiation is scattered in all directions, so at the surface of the earth the sunlight has both a direct and a diffuse component. At the altitude at which contrails form, the incident sunlight may be considered to have only a (Jirect cornpotlernt. When this direct radiation strikes the top of a contrail layer it is scattered in all directions. The fraction of it which is scattered back to space is called the solar albedo of the layer. Since the absorption of solar radiation by the ice particles can be neglected, the albedo of the layer is equivalent to the fractional reduction in the solar radiation received at the surface of the earth. The solar albedo of a contrail layer will be calculated in this chapter, using the single-scattering properties 21

22 derived in Chapter 2. To simplify the calculations, the layer will be assumed to be uniform and of infinite horizontal extent. Most of the terrestrial radiation originates at the surface of the earth and in the troposphere. A contrail layer located near the tropopause absorbs part of the outgoing terrestrial radiation, and also emits radiation in both directions. The net effect of the layer is to decrease the amount of long wave radiation lost to space and to increase slightly the amount of long wave radiation absorbed at the earth's surface. To evaluate the magnitude of these effects the infrared transmissivity and emissivity of the layer must be calculated. Generally the backscattering of long wave radiation by ice particles is small enough so that the reflectivity of the layer may be neglected. The infrared transmissivity and emissivity of the solar albedo are calculated by solving the radiative transfer equation with the proper boundary conditions. The radiative transfer equation is discussed in the following section. In Section 3.3 two approximate solutions are derived for the case of diffuse radiation incident on one side of a scattering and absorbing layer. One of these (the two-stream approximation) is then used to calculate the infrared transmissivity and emissivity of a contrail layer. In Section 3.5 an approximate solution is derived for the case of direct radiation incident on a scattering layer. This solution is then used to calculate the solar albedo of a contrail layer.

23 3.2. THE TRANSFER EQUATION The fundamental variable describing a radiation field is the intensity, I, which is a function of position and direction. (It is also a function of the vwavecillgtl, but in the following the wavelength dependence will be igrlored, and it will be assumed that all quantities refer to the same wavelength interval.) The problem to be considered involves a plane parallel layer containing a uniform distribution of scattering and absorbing particles. Thus the appropriate position variable is the distance z from one side of the layer, and the direction variables are the polar angle 0, referred to the z axis, and the azimuthal angle $. The intensity is defined such that I(90,)dQ is the amount of radiant power confined within a solid angle dQ about the direction (0,P) and passing through a unit area normal to that direction. The amount of power passing through a given unit area from all directions is called the net flux; If the unit area is parallel to the sides of the layer, the hemispheric flux is given by 27r ~1 F+ = d d/l (3.1) where p = cosO and the two signs refer to the two sides of the area. The radiative transfer equation describes the change in the intensity due to scattering and absorption within the

24 z T=O T =-TI I(T, pL) Figure 4. Geometry of scattering layer with optical thickness T1, showing diffuse radiation incident on the bottom of the layer (T= T1) and emerging from the top of the layer (T = 0).

25 layer. The layer is assumed to be made up of scattering and absorbing particles, with a volume number density n. The single-scattering properties of each particle are described by the phase function ~(9,~), and the cross sections os, a'3 and at discussed in Chapter 2. For spherical particles and unpolarized radiation there is no azimuthal dependence, so the phase function may be written 1(cos6). ~( cosG)dQ is proportional to the probability that a photon interacting with a given particle will be scattered into the solid angle dQ about the scattering angle 0. It is usually normalized as follows: f( cosO)dQ = 4r w 3.2) where wo = as/at is the single-scattering albedo, and the range of integration is over all angles. The probability per unit solid angle for a photon to be scattered from the initial direction (p,4) into the final direction (',' ) is given by the scattering function p(,4;lp',') = f {p' + (1-112) (1-' 2)~ cos(4-' )} (3.3) The argument of P in this expression is just the cosine of the scattering angle, defined in Chapter 2. Notice the symmetry in the primed and unprimed variables, which is required by electromagnetic theory. The radiative transfer equation is written by Chandrasekhar (1960) in the following form:

26 2f 1 Iujy I(Iji) = i(T,H,() - d J P(IJ,;l' I(Tp',p )dp 0 -1 (3.4) The first term on the right-lland side of this equation accounts for reduction of the itltensity in thle (p,4) direction by scattering and absorption, while tlhe second term accounts for scattering of radiation into this direction from all other directions. Notice that the geometrical variable z has been replaced by the optical depth T = noz. The sign convention used here is such that a photon traveling in the positive T direction has P < 0. An equation for the intensity averaged over azimuthal angles may be obtained by integrating (3.4) over 0 < $ < 2rr. The result is: -dT I(T,P) = I(T,p) - -- P(i,') I(T't)dp' (3.5) -1 whe re 2rr I (~,u) = +fI(T,p,P)dq (3.6) 0 and 2r P(,pJ' ) =',1, fP(p,;?,)'' )dq (3.7) () NoLt that thlis irntegral is indep(endent of 9' since the scalttering; functiorl is periodic in ~-~'

27 If the phase function is expanded in Legendre polynomials, 00 ( P) = EkP (P) (3.8) Z=0 then by the addition theorem of spherical harmonics (Chandrasekhar, 1960, p. 150), the scattering function averaged over azimuthal angles may be written as follows: 00 P(p,p') = kQP (9) P('). (39) Q=0 In the following sections, the intensity is used only to calculate the hemispheric flux via Equation (3.1). Therefore, the azimuthal dependence of the intensity is not needed and the transfer equation in the form (3.5) will be used thlrougho ut. 3.3. SOLUTIONS FOR DIFFUSE INCIDENT RADIATION Solving the radiative transfer equation is in general a rather difficult procedure, and no exact analytical solution exists for a finite layer with arbitrary scattering properties. One method for obtaining an approximate solution is to convert the original equation in the two continuous variables 1I and T into a set of N equations in the single variable T. This may be done in either of two ways, either by expand

28 intg thle intensity in a series of Legendre polynomials N- 1 [(L,) = I ) P () (3.10),= 0 and solving for the coefficients I (T), or by solving the radiative transfer equation at N discrete values of p: I (T) = I(T,pQ) = 1,2,...N (3.11) In the first case, the orthogonality of the Legendre polynomials is used to remove the integral in the radiative transfer equation. In the second case, the,9~ values are chosen to be the Gaussian quadrature points (Chandrasekhar, 1960, p. 71) and the integral in the transfer equation is approximated by Gaussian quadrature. In order to calculate the diffuse transmissivity and reflectivity of a layer of optical thickness T1, the radiative transfer equation is solved for 0<T<T1 with the boundary conditions of unit flux incident on the bottom of the layer (T = T1) and zero flux incident on the top (T = 0). The transmissivity is then equal to the flux leaving the top of the layer, and the reflectivity is equal to the flux leaving tile bottom. Simple arialytical expressions for the transmissivity arid reflectivity miay bo oota;ined by using eithe~ of the abovementioned procedures, witi1 N = 2. Tihe first method reduces

29 to the Eddington approximation (Shettle and Weiniiian, 1970) and the second becomes the two-stream approximation (Sagan and Pollack, 1967). Consider first the Eddington approximation I(T,p) = IO(T) + pIl(T) (3.12) Substituting this into the transfer Equation (3.5) and taking the first two moments of the transfer equation with respect to p yields the following pair of simultaneous differential equations: dIl(T ) 1dT — 3(1-w ) Io (T) dI (T) W1 d = (1- ) I1(T) (3.13) where w0 and wl are the first two coefficients in the Legendre expansion of the phase function. Making the change of variable T = [3(1-wo) (1- - )] T (3.14) Equations (3.13) become

30 dli(T' ) i = al (T') d-t' whe re 2 3(1-w ) a2 = (3.16) Wi TIle solution of this set of equations may be written o(T') = e A2 e I () =ae - aA2 e (3.17) The boundary condition F_(O) = 0 yields the following relation: A1 + A = a(A-A2) (3.18) The reflectivity t 3 and transmissivity T are then given by the following expression (see Appendix B):

31 TI -T1 (1-b2) (e 1 - e 1(.19) (l+b)2e (1b )2e 4b T t= T. 4b -T (3.20) (1+ b)2e 1 (1-b)2e 1 where = [3(1-wo) (1- )]T (3. 21) and b2 4 2 4 1- 0~ b2 = a - (3.22) 1- - From energy conservation considerations, the absorptivity is given by A = 1 - R - T. (3. 23) The Eddington approximation gives a negative reflectivity when b > 1, or equivalently when wo < (i+f1). That is, the approximation breaks down when absorption is large (Lo is small) or the scattering is strongly peaked in the forward direction (Uo is large). These conditions tend to

32 produce a high degree of anisntropy in the intensity which cannot be accurately represeint ed by the truncated Legendre exparlnsion (3.12). If tile intensity represented by this e xpanlsion is positive for all 1, as it must be, the maximum ratio of the flux in the f'orward hemisphere to that in the reverse hemisphere is 5. (The intensity is non-negative for all pl only if 1I11 < 10, which implies that F+/F < 5). Whenl either the absorption is large or the scattering is;;trorigly peaked in thie forward direction, very little of the incident flux is scattered into the reverse direction, and this maxiiriuIl ratio is exceeded. This leads to negative intensities, thereby causing the approximation to break down. This difficulty is avoided in the two-stream approxiiiationl, where the intensity is obtained at two discrete values of p and there is no restriction on the relative values of the intensity at these points. The points are chosenl as 1;Ale fir'st-order Gaussian quadrature points (~+ K1 = + 1/3) so thlat the integral in the transfer equation carl be approximated by the corresponding quadrature formula. Gauss's formula is thie optimum quadrature method in the sen"se that it is exact for integrands which are polynomials of the Iighest possible degree, for a given number of quadrature points (Chandrasekhar, 1960). The first-order Gaussian quadrature formula, whlich is exact for polynomials of the third degree mnay be written as follows:

33 1 fF( P)d1 ~ F(-i1) + F() (.24) The quadrature poirtlts are also the zeros of the second-order Legendre polynomial: P2 (+~1) = 0. Neglecting terms of order >_ 3 in the Legendre expression of the phase function (3.8), the scattering function evaluated at these points is then W P(l'+U =l) P(-1,'T1)1) =o - (3.25) The transfer equation evaluated at these two points may therefore be written as follows: dI (T) 1 W1 _1 ( d() = I(T) W _ - -(W +-)) I(T) 1 dTo 3 2 3+ 3 dI (T) t t =1(T Wo1 d =) I_ (T) T) - o+ I (T) (3. 26) where I +(T) - I(T,~1) (3.27) These equations are exactly equivalent to the Eddington Equations (3.13) under the transformation I+(T) = Io(T) + piI1(T) (3.28)

3)4'l'he so()lutions may tiler el'ore be taken over directly from Equations (3.17) and rewritten in the form I+(') = A(1+b) e' A2(l-b) eeT T t -T' (3.29) I_ (T) = A1(1-b) e' + A2 (l+b) e-T' (29 where 1-w = 2 a2 0 (3.30) 3 tt) In the two-stream approximation, the flux in either heiriisphere is assumed to be proportional to the intensity at,I = + ~l' so the boundary condition is taken as I (O) = 0 and thle definitions of reflectivity and transmissivity redu eC to ) i (T1) I+(O) + I+( I+(T1) (331 Applying this boundary condition yields the relation A1 b+l 1 (3.32) A2 b - 1 Substituting the solutions (3.29) into the definitions of R and Tri then results in the same expressions (3.19) and (3.20)

35 for the reflectivity and transmissivity as obtained for the Eddington approximation. The results of the two approximations are not the same, however, because the parameter b is defined by Equation (3.30) in the case of the two-stream approximation, and by (3.22) for the Eddington approximation. The Eddington and two-stream approximations are compared with exact results for three special cases in Figures 5 through 7. Figure 5 shows the transmissivity plotted versus T1 for the case = 0 (complete absorption). In this case the transfer equation can be solved exactly to yield the transmissivity T = 2 pexp(-T3l/)dP (3.33) The Eddington and two-stream approximations give very nearly the same results for the transmissivity in this case, and both are quite close to the exact result. Figure 6 shows the reflectivity plotted versus T1 for the case = 1 (conservative scattering, no absorption). The approximations in this case are compared with the numerical computations of Hansen (1969). Note that the reflectivity in this case is a function only of (1 - g)Tl, as pointed out by Piotrowski (1956). The reflectivities calculated by the two approximate methods are again quite close, with the twostream approximation being somewhat more accurate for small

36.8.6 T WoO:.2 -EXACT SOLUTION: T 2E3('z,) —. —EDDINGTON AND TWO-STREAM APPROXIMATIONS " oL l.01 0.1 1.0 Figure 5. Transmssvity versus optical thickness Figure 5. Transmissivity versus optical thickness T for the case w = 0 (perfect absorption). E3(T1) is the third-order exponential integral.

37 1.0 I |-HANSEN'S NUMERICAL COMPUTATION, - |. -TWO- STREAM APPROXIMATION.8 — EDDINGTON APPROXIMATION -.6 / R.4:/ /.2,1 0 O.I I I 0.1 1.0 10. (I-g) T, — Figure 6. Reflectivity versus (1-g)T1 for the case wo = 1 (perfect scattering). g = w1/3wo is the asymmetry factor of the phase function, and T1 is the optical thickness.

38 1.0.8 - EXACT SOLUTION (CHANDRASEKHAR). —- TWO-STREAM APPROXIMATION - -- EDDINGTON APPROXIMATION T.- /;,.6 R scpattering..4.2 0.2.4 Wo.6.8 1.0 Figure 7. Reflectivity versus single-scattering albedo w for a semi-infinite layer with isotropic scattering.

39 T1 and the Eddington approximation better for large r1. The third case is the reflectivity of a semi-infinite layer (T1 = a) with isotropic scattering (w = 90, > 1). This problem can also be solved exactly, with tile following result: R = 1 - 2(1-w) ml(o ) (3.34) where 1l(wo) is the first moment of the H-function tabulated by Chandrasekhar (1947, 1960). This is plotted in Figure 7 along with the values calculated from the Eddington and twostream approximations. The two-stream approximation is considerably better in this case; the reflectivity given by the Eddington approximation is too small for all values of wO, and is negative when w < 1 3.4. INFRARED PROPERTIES OF CONTRAIL LAYERS The results of Section 3.3 indicate that the twostream approximation is somewhat more accurate than the Eddington approximation when absorption is appreciable. The two-stream approximation was therefore used to calculate the properties of a contrail layer for terrestrial radiation. The results of this approximation have been expressed in terms of the first two Legendre expansion coefficients wO and l of the phase function and the optical thickness T1 of the layer. These are related to the single-scattering parameters defined in Chapter 2 as follows:

1 10 3 o T1 = QextT* (3.35) whle:re'r* = n7rr2t for a layer of geometrical thickness t 1havingt a concentration n of particles with average radius r.'1the transinissivity Ti and absorptivity A. were calculated from ]equations (3.20), (3.21), (3.23), and (3.30) using thie single scattering parameters listed in Table 1. Calculatiosns were carried out separately for each particle size (r = L, 3, and 10 microns), and for each incident wavelenlgth inlterval (i = 1,...,5). The average transmissivity of the layer over the entire terrestrial spectrum is given by = 1 iH Ti (3. 36) i=l where Ti is the transinissivity for band i, Hi is the incident flux of terrestrial radiation in band i, and Ht is the total flux in all five bands. The values of Hi/Ht for a typical terrestrial spectrum arce indicated in Table 2. According to Kirchhoff's law, the flux emitted by the layer in band i is equal to the absorptivity A. multiplied by 1

41 TABLE 2 WEIGHITING FACTORS FOR INFRARED TRANSMISSIVITY AND EMISSIVITY Band i 1 2 3 4 5 Wavelength 4-8 8-12 12-20 20-40 40-100 (microns) Hi/Ht.043. 359.286.236.068 Bi(Tc)/oTc.035. 157.357.335.105

4 2 tthe rFlux c-mnitted by a bl]aclk body at the same temperature.'rlhe averajge emissivity of the layer is therefore given by 5 f =,1 — ~A1 (T) A(3.37) Col 4 c wh',r'e a i; thtee i; l a,(Ftrl-l(l ~ tznl;arln (con:ta.nta, T it, the temperalturle of the corntrail layer, a"n(i.i(T] ) is tile black body 1 c l ux in band i. Trhe values of B.i(T )/aT c are indicated in')ta'le 2 for:a cont;rail temperature of T = 2200K.'lThe average tr;ranismissivity T and emissivity E are plotted versus T* in Figure 8, for each of the three particle sizes considered. By comparing these with the transmissivity measured by Kuhn (1970), an estimate of the value of T* for:un average or typical contrail layer can be obtained. Kuhln measured an average transmissivity of about 0.90 for the contrail layers he observed over California and in the West Indies. Assuming that the thickness of these layers was 500 meters, and the particle radius was 50 microns, he cal-culated a particle concentration of 0.027 cm-3. This implies tlhat T* = 0.106 and the condensed water content is 7 x 10-4gmn/cm2 for these layers. From the discussion in Chapter 2 it appears that a better estimate of the particle radius in contrails is about 3 microns. For tlhis particle size, the observed transmissivity is reacied. at about i[* = ().2, or twice the value calculated by Kuhn. Because of the smaller particle size, however, the water content of suchi a layer would be only 8 x 10-5gm/cm2,

43 1.0 -.6 I.R. TRANSMISSIVITY r -OL I.R. EMISSIVITY l r 3 0.1.2.3.4.5 Figure 8. Diffuse transmissivity and emissivity of a layer of ice particles versus reduced optical thickness T* = nTr2t, calculated from the two-stream approximation.

44 almost an order of magnitude smaller than that implied by Kuhn's model. Although 3 microrls seems to be the most likely particle size, calculations are carried out for r = 1 micron and r = 10 microns as well. The results of these calculations give an indication of the uncertainty in the effects of a conltrail layer due to variations in the particle size. The samei value of T* is used in each calculation, so that the effects on the solar radiation are the same in each case. The infrared transmissivities corresponding to T* = 0.2 are 0.968 for r = 1 micron, 0.887 for r = 3 microns, and 0.726 for r = 10 microns. The reflectivities (not shown in Figure 8) are 0.001, 0.023, and 0.042, respectively. ~. *. OLuIJr[lIONS FOrl DTIREliCT INCIDENT RADIATION A beacnl off dir ect radiation, in which the photons are,ll t rlveilinl: In the direction (-p,4 ), may be described by I1 ho' i]nten:I. ty 00 0 (,= I6(P+ [1) 6( ffO) (3.38).If such a beam falls on a horizontal layer of scattering particles, the flux incident on the top of the layer (T = 0) is F(O) = p (3.39) 0

45 The flux of direct radiation at an optical depth T below the surface of the layer is then F(T) = vTrE e-T/ ~ (3.40) The flux removed from the direct berr is, partly absorbed, if (o < 1, and converted into diffuse -radiation. The intensity of diffuse radiation at any point is obtained by solving the radiative transfer equation. If the azimuthally-averaged diffuse intensity is indicated by I(r,i), the transfer equation may be written (c.f. Chandrasekhar, 1960, p. 22) as 1 d = I(, ) P(pp') I(T, p)dp' WdyI(T,)) -1 T - E[ E e ~p(~ -uJ) (3.111) where p(P,-Po) is the azimuthally-averaged scattering function, defined by Equation (3.7). The last term represents the scattering of radiation from the direct beam into the direction P, while the integral represents scattering from the direction I' into the direction 4. If there is no diffuse radiation incident on the layer,

46 the boundary conditions are I((O,) = 0, for p < 0 I(1l,1) = 0, for P > 0 (3.42) whire r1 is the total optical thickness of tile layer. IPThe flux of diffuse radiation escaping from the top of the layer is 1 = 27T PI(0,p)du (3.143) o The albedo of the layer is then defined as A(T1, ) - IF',/F(0) = 2 ffI(O, 0 )d (3. 44) 1o Ute solution of the transfer Equation (3.41) is the Neumann series n is byIrvine= E1 I Ths my be p -(3.45) n==1 which is discussed by Irvine (1965b, 1968). This may be physically interpreted as an expansion in successive orders of scattering. That is, I(n) (T,) is the intensity of radiation which h[las been scattered n times. Th-e derivation of the Neurnann series may be followed most easily by r(eferrirng2 to the transfer equation in the form of E:quatiorn (3.5).'I'he integral term in this equation accounts

47 for the scattering of radiation from the direction p' into the direction 11. The first-order intensity is therefore obtained by substituting the intensity of direct radiation -1/pa Id(Ti) = 2 E 6(p+io) e (3 46) into the integral term of the transfer equation. This gives the equation d (1) (1) udT I (T,v) = I ( E e ) (3.47) which may be solved by using the integrating factor e T/P The solution, with the boundary conditions (3.42), is (T ) =e ( 1-exp + )( T (3.48) where T, for p < 0 (349) ~ T1 ) for u > 0 The nth order intensity I(n)(T,u) is obtained by substituting I(n-l)(T,p) into the integral term of the transfer equation:

1 8 d1()_10, ( 5n-1 ) (3.50) for thin layers, most of the radiation escapes after having been scattered once. Therefore, the single-scattering Intensity I (,1)(I) dominates and the Neumann series converges rap i dly. b'or thick layers, multiple scattering effects dominate and the Neumann series converges very slowly. In either case, the terms of order n > 2 must be numerically computed, since there is no analytical solution for Equation (3.50). As a result, this method is time-consuming and impractical for most applications. An approximate solution may be obtained by using the Eddington approximation for the multiple scattered radiation, while retaining thle exact solution for the first-order (sinlgly-scattered) intensity. That is, the intensity is written in the forin Ik( r,) = I(1) (TP) + Im(T,P) (3.51) where Im(T,d) is assumed to be of the form Im(T,U)= I((T) )+ Il(T) (3.52)

49 and the transfer equation is solved for I (T) and II(T), using the boundary conditions (3.42). This approximation is sul;gested by the fact thlat the first-order illtellti;ity I (1)T,p) is highly anisotropic, while the higher-order terms have a smoother angular variation. The results are substituted into the definition of the albedo: 1 A(Ti,Po) I Ef iI (o,p)dp + E I o() + o I1(0) (3.53) O0 Details of the derivation are given in Appendix C. The results may be written as follows: 1 1'(Tl'o=) =1 + 00 I f + 6f2 p - 2 0 (3.54) whe re 3 3 (1W13 x=-g) - 3-e x =- e -T /0o x= 1 - eFl/l f= [(1- p) 00 ~ 2+-2~

30 Tl'is function was numerically computed using Simpson's rule witht N = 40 subintervals. The phase function was approximated by tile fenyey-Greenstein function (2.4), which has the followiln-, l,et?;enddre exparnsiorn coefficients (van de Hulst, 1971): o = o0 (2 + 1)gZ (3.56)'t'hie first 40 terms of' the Legendre expansion were used in ttie numerical evaluation of the albedo (wO = 1 for the case of conserlvative scattering under consideration). In Figure 9 the results of this approximation are compared with exact numerical computations of the albedo by Hansen (1969) for g =.844 and T1 = 0.5, 4, and 16. The approximation gives quite accurate results, considering the mirtinlal computer time required for its evaluation. 3.6. SOLAR AL1EDO OF CONTRAIL LAYERS In Chapter 2, the single scattering properties of contralil ice particles for sunlight were taken to be Q =2.0 e xt co = 1.0 g = 0.85 (3.57) for all three particle sizes considered. For a typical layer, with T* = 0.2, the optical thickness for solar radiation is then

51 HANSEN J969) 1.0 VAP,%SEN IMAXION 4S.644 ~~~CI ~ ~ sf'75 t+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~at 0.5~~~~~~~~~~~~~~~~~~~~~~ lVp;

r = ext T=.= (3.58) The solar' albedo of sucth a layer was computed by the approximate method described in Section 3.5. The results are plotted versus the cosine of the solar zenith angle in Figure 10. Kuhn (1970) made airborne radiometric observations of;several conltrail layers, and reported an average solar albedo of' 15 percent. T'l'his agrees with the calculations presented here if it is assumed that the solar zenith angle was about 700 at the time of his measurements. The solar zenith angle 0O = cos-'(,o) is given in terms of' the latitude $, the solar declination 6, and the hour angle 11 by the following formula:,> = sirnd sin6 + cost cos6 cosh. (3.59) At local nIoon, h = 0 and therefore 11 = cos(q-6). (3.60) At sunset, 0o = O and h = II _ cos-'(-tan~ tan6). (3.61) Trle solar declination varies from -23.5~ at winter solstice to +23.5~ at summer solstice. If d is the day number, beginning with dC = O at the vernal equinox, thie solar declina

53 O8 QEXT s 2.0 Gs 0.85 ~* g 0.2.6*.6.8 0 -M i 0.-.2 p 2:COS e0 Fligure 10. Local solar albedo versus the cosine of the solar zenith angle 6o for -t 0.2, comPuted from Equation (3'54) with Qext 2.0 and g 085

5" tiorn is givenl by 27rd sin6 = sin(23.5)si n( —). (3.62) The amount of solar power falling on the top of the atmosphere at any given time is PoIo, where Io is the solar -2 — 1 constant (I = 2.0 cal cm min ). Neglecting absorption of sunlight in the atmosphere above the tropopause, this is also equal to the power falling on the top of a contrail layer. The amount of solar power reflected from the layer is then PoIoA(T10,o)The daily average albedo of the layer (assuming that it persists all day long) is therefore given by H fO0A(T1},po) dh A_= (3. 6 3) H P odh This average albedo was computed at 10~ latitude intervals in the northern hemisphere for the spring (6 = 90), summer (6 = 18~), fall (6 = -9~), and winter (6 = -18~). The results for the summer and winter are plotted in Figure 11, and all four seasonal values are tabulated in Table 3.

55.5 WINTER o.3 a MER 2 0,_ j, I, 8,,SU MM 0 01' 20 30 40 50 60 70 80 9(r LATITUDE Figure 11. Average solar albedo of contrail layer with T* = 0.2, for sumrmer (6 = 18~) and winter (6 = -18o).

56 TABLE 3 AVERAGE SOLAR ALBEDO OF CONTRAIL LAYER (T* = 0.2) Latitude (ON) Spring Summer Fall Winter 5.045.048.048.053 156 046.046.054.064 25.049.047.065.082 35.057.051.084.112 115.070.060.117.167 55.092.075.176.272 65.130.100.286.456 75.200.157.479 -- 85.298.184

57 The large latitude and seasonal variation in the average albedo of the layer is due to the strong dependence of the local albedo on the solar zenith angle (cf. Figure 10). This dependence results from the increase in the path length of direct radiation through the layer and the decrease in the fraction of radiation scattered downward with increasing solar zenith angle. The average albedo at low latitudes is smaller than the figure of 0.15 used by Kuhn (1970), but increases to values much larger than Kuhn's albedo at high latitudes. This has important implications for the effects of contrails in arctic regions, as will be shown in Chapter 4.

CHAPTER 4 CLIMATIC EFFECTS OF CONTRAILS 4l.1. RADIATION BALANCE OF THE SURFACE AND ATMOSPHERE The solar and terrestrial radiation fields in the atmosphlere are coupled with each other and with the physical state of the atmosphere. They also determine the rates of processes on the surface of the earth, such as the evaporation of water and the melting of snow and ice. The interactionl of all these variables is extremely complex and is not well understood. In this section, the initial changes in the radiation fields due to the formation of a contrail layer are calculated, assuming that no other changes take place in the atmosphere or at the surface of the earth. In the followill?- section the resulting changes in the surface temperature, ev;aporation rate, and melting rate will be estimated using cor tain other simplifying assumptions. In order to evaluate the radiative changes caused by a contrail layer, some estimate of the preexisting conditions must first be made. The radiation fields in the atmosphere vary with time and place, having both a periodic component due to the motion of the earth relative to the sun and a random component due to changes in the physical state of the atrnosphere. The radiation fluxes at the bottom and the top of the atmosphere averaged over a three-month period were estimated by London (1957) in a study of the atmospheric heat 58

59 balance. These estimates are used here as initial conditions to calculate the changes induced by contrails. Among the variables estimated by London are: the incoming radiation absorbed in the atmosphere, incoming radiation absorbed by earth's surface, net long wave radiation from the earth's surface (or nocturnal radiation), and long wave radiation lost to space. These are tabulated at ten degree latitude intervals for each season, in units of cal cm-2 min- 1 The net downward flux of solar radiation is denoted here by S, and the net upward flux of terrestrial (long wave) radiation by L. The flux through the top of the atmosphere is indicated by the subscript t, and through the bottom of the atmosphere (the surface of the earth) by the subscript s. Therefore, the variables tabulated by London are symbolized as follows: Incoming radiation absorbed by earth's surface = S Incoming radiation absorbed in the atmosphere (Sa) = S -S t s Long wave radiation lost to space = Lt Net long wave radiation from the earth's surface = Ls. The net radiation of the atmosphere (Ra) is defined as the net amount of radiative energy gained per unit time by a vertical column of unit cross section extending from the bottom to the top of the atmosphere. The net radiation of the atmosphere may therefore be written as

60 a Sa + Ls - Lt (4.1) The net radiation of the earth's surface (Rs) is given by R s - L. (4.2) In most cases, the net radiation of the surface is positive, while that of the atmosphere is negative. That is, the surface absorbs more radiative energy than it emits, and the atmosphere emits more than it absorbs. In order to maintain thermal equilibrium, energy is transferred from the surface to the atmosphere by convection and by evaporation of water which condenses in the atmosphere. A contrail layer at the tropopause may be considered for the present purposes to be effectively located at the top of thie atmosphere, since most of the radiative energy is absorbed and emitted in the troposphere. Therefore the effect of such a layer on the solar radiation is to decrease St by the factor A, where A is the average solar albedo of the layer. If all wavelengths of the solar radiation are equally reflected by the layer, and the absorption properties of the atmosphere do not change, the net flux of solar radiation at the surface (Ss) is also reduced by the factor s A. The effect of the layer on the long wave radiation is more complex, since the layer both absorbs and emits long wave radiation. To simplify the calculation, it is assumed

61 that reflection of long wave radiation may be neglected (cf. Chapter 3) and that the layer is in radiative equilibrium (cf. Appendix D). The layer then absorbs an amount of long wave radiation (1-T)Lt, where T is the average transmissivity of the layer, and emits half of this amount in each direction. The net decrease in Lt is therefore ~(1-T)Lt. The decrease in Ls due to radiation from the layer is given by ~(l-T)LtTa, where Ta is the transmissivity of the atmosphere for long wave radiation. This transmissivity is only about five percent under average conditions, according to London (1957). Therefore most of the radiation emitted downward by the layer is absorbed in the atmosphere. Combining these results, the increase in the net radiation of the atmosphere due to the formulation of a contrail layer is AR = ~(1-Ta)(l-T)Lt - A Sa (4.3) while the increase in the net radiation of the earth's surface is ARs = ~Ta(1-T)Lt - A Ss (4.4) These quantities were calculated using London's tabulated values for Lt,, S Ss, and the transmissivity T and daily average solar albedo A obtained in Chapter 3. Calculations were carried out at ten degree latitude intervals

62 for each season. An atmospheric transmissivity (Ta) of five percent was used throughout. The results of these calculations are plotted in Figures 12(a) through 12(d). The upper curve represents the change in the net radiation, or radiation balance, of the atmosphere for a layer of particles with r = 3 microns. The lower curve represents the change in the surface radiation balance for the same layer. The upper error bars in eaclh case represent the effects of a layer with r = 10 microns, while the lower error bars are for r = 1 micron particles. Thus the error bars indicate the uncertainty in the effects due to variations in particle size. The change in the atmospheric radiation balance is dominated by infrared effects for large particle sizes, and by solar effects for small particle sizes. In most cases ARa is positive, indicating an atmospheric heating effect, but the magnitude of this effect cannot be estimated without furthier knowledge of the particle size. The change in the surface radiation balance is dominated by solar effects in all cases, and is always negative. This would imply a reduction in temperature and/or the rates of evaporation and melting beneath the contrail layer. These changes are estimated in the following section, neglecting changes in the convective heat transfer between the surface and atmosphere. The seasonal and latitude variations of the effects calculated here do not agree with Appleman's (1966) conjecture

63.05 -.04 SPRING.03 02 NE.01 -.01 -.02 -.03 -.04 0 10 20 30 40 50 60 70 80 90 LATITUDE (~N) Figure 12(a). Effect of a continuous contrail layer on the spring radiation balance of the atmosphere (upper curve) and surface (lower curve).

64.05.*-04 -. SUM MER.04.03 02 N.01'E C) 1 — 01 -.02 -.03 0 10 20 30 40 50 60 70 80 90 LATITUDE (ON) Figure 12(b). Effect of a continuous contrail layer on the summer radiation balance of the atmosphere (upper curve) and surface (lower curve).

.05 _.04 FALL.03 02 N-.01'E -.01 -02 -.03 -.04 I o 10 20 30 40 o 60 7o so80 90 LATITUDE ("N) Figure 12(c). Effect of a continuous contrail layer on the fall radiation balance of the atmosphere (upper curve) and surface (lower curve).

66.05.04 WINTER 03 02 0C -.02 -.03 -.04 0 10 20 30 40 50 60 70 80 90 LATITUDE (ON) Figure i2(d). Effect of a continuous contrail layer on the winter radiation balance of the atmosphere (upper curve) and surface (lower curve).

67 that a contrail layer would cause a warming of those latitudes where the outgoing radiation exceeds the incoming, and a cooling where the reverse situation exists. The variation observed in Figures 12(a) through 12(d) is more nearly in the opposite direction. This is due to the increase in the solar albedo of the layer with increasing solar zenith angle. In winter, and at high latitudes, the incoming radiation is smaller but the solar albedo is larger, so the net effect of the layer on the incoming radiation can actually be larger in the winter than the summer, or at high latitudes than at low latitudes. 4.2. SURFACE CHANGES INDUCED BY CONTRAILS The calculations presented in the previous section indicate that the presence of a contrail layer in most cases increases the radiation balance of the atmosphere and decreases the radiation balance of the earth's surface. As a result of these perturbations, the temperature, rate of evaporation, and other climatic variables will change in order to restore the thermal equilibrium of the surface and atmosphere. The general problem of predicting these changes has not been solved, due to the complexity of the interactions between climatic variables. However, in some cases the change in one climatic variable may be estimated by assuming that the others remain constant. For example, over water or wet soil it might be assumed that the change in the surface radiation balance produces a change only in the evaporation rate, with

68 out affecting the suirface temperature, The change in the evaporation rate hE would then be given by AE = AR /L (4.5) whe-re I. = 590 cal gm-lis the latent heat of vaporization of water. For ARs -.016 cal cm-1min-', a typical value at midlatitudes, the evaporation rate would be reduced by about 0.4 mm per day (about a 20 percent reduction in the average evaporation rate). Over dry land, it might be assumed (as is done by Kuhn, 1970) that the presence of a contrail layer causes only a heating or cooling of the surface, without affecting the other climatic variables. In this case the surface temperature would change just sufficiently to return the surface radiation balance to its original value. The outgoing long wave radiation at the surface is given by L, () = aT4 (4.6) S s where a = 8.12 x 10-1l cal cm-lmin-l~K-4 and T is the mean surface temperature. If the mean surface temperature changes by ATs, the emitted radiation changes by ALs () = 4oTsATs (4.7) The change in the mean surface temperature required to return the radiation balance to its original value is therefore

69 ATs AR /4aT. (4. 8) 5 s S This formula gives the maximum temperature change which could be expected, disregarding the possibility of positive feedback machanisms which could intensify the effects of a contrail layer. The change is on the order of 2 ~C at rnidlatitudes, increasing to a maximum of 3.5 ~C in the arctic spring and summer (note, however, that this calculation assumes a dry surface, which rarely occurs in the arctic). It should be repeated that these results are for a continuous and complete contrail cover. In actuality, of course, the cover is never continuous or complete. In order to estimate the average temperature change over a period of time T (at least one day), the temperature change indicated in Figure 13 should be multiplied by the factor f flf2, where fl is the fraction of the time period under consideration during which contrails were observed and f2 is the fraction of the sky covered by the contrails while they lasted. If N distinct contrails were observed during this period, the factor f may be expressed as follows: N f =1ti C (4.9) i=l where ti is the lifetime of contrail i and Ci is the fraction of the sky covered by this contrail. A single contrail is sometimes observed to cover a tenth of the sky, and conditions suitable for the formation of such contrails may persist all day. Values of f from one-tenth to one-half can

70 -4 -3 SP ING - SUMMER I~ -2 -I - \FALL- WINTER 0 10 20 30 40 50 60 70 80 90 LATITUDE (ON) Figure 13. Changes in surface temperature for a continuous contrail layer over dry land.

71 therefore occasionally occur in some locations over a period of one or two days. Another climatic effect, which might be of importance especially in arctic regions, is the influence of contrail formation on the rate of melting ice and snow. In the spring, when ice and snow are being melted, the surface temperature would be expected to remain fairly constant until melting is complete. Therefore any changes in the surface radiation balance would have little effect on the surface temperature, but would directly influence the rate of melting. If the radiation balance is changed by an amount ARs during a certain period, the average melting rate M will change during this period by the amount AM = AR /F (4.10) where F = 80 cal gm-' is the heat of fusion of water. During the arctic spring and summer, Figures 12(a) and 12(b) indicate that the surface radiation balance would be reduced by about 0.02 cal cm-2min-1 by a continuous contrail cover of moderate thickness. [Jacobs (1971) measured a decrease in the surface radiation balance of 0.1 cal cm-2minfor a period of seven minutes caused by a single dense contrail over Baffin Island.] For a coverage factor f = 0.01 during the six-month period from March to September, the total reduction in energy density available to melt ice during this period would be about 50 cal cm-2. According to Flohn (1969), the area of ice in the

72 Arctic Ocean decreases from 18 million km in March to 10 2 million kmrl In September. This ice is mostly in the form of floes about 250 cm thick, so the amount of energy required to melt this ice is about 104cal cm2. From the preceding calculations it appears that a one-percent contrail over the arctic region could increase the arctic ice cover by about one-half of one-percent per year. This could have an appreciable climatic impact in view of the possible role of the arctic ice cover as a positive feedback mechanism for climatic change (Budyko, 1972). Briefly, this is due to the fact that snow and ice have a much higher solar albedo than water. Any growth in the ice cover therefore leads to a reduction in the amount of solar radiationl absorbed, which causes a further growth of the ice cover and a further cooling of the earth as a whole. This theory has been proposed by Budyko as an explanation for the appearance of ice ages. If the theory is correct, it is of vital importance to guard the arctic ice cover from external influenlces, because of the possibly catastrophic results which might ensue from a change in its size.

CHAPTER 5 CONCLUSIONS 5.1. EXTENT OF CONTRAIL FORMATION The atmospheric conditions required for contrail formation were derived by Appleman (1953), and are summarized in Figure 14. According to this figure, contrails will form at cruise altitude even in a completely dry atmosphere if the temperature is below -53 ~C. In this case, however, the contrails will evaporate as soon as the exhaust mixes with the ambient atmosphere. The contrails will persist for a significant length of time only if the atmosphere is saturated over ice (roughly 60 percent relative humidity). The latitude regions where there is a probability of contrail formation greater than 50 percent are shown in Figure 15. Conditions favorable to contrail formation occur most frequently at high latitudes because of the low atmospheric temperatures there. This is significant because of the relatively large magnitude of the contrail effects at these latitudes, and the climatic importance of the arctic ice cover. Figure 15 does not indicate the frequency with which persistent contrails form. However, a recent study by Beckwith (1972) provides some information on contrail persistance. According to this study, which is based upon 395 observations over the continental United States, persistent 73

74 Pressure Pressure Altitude (mb) (1000's ft) 40 SST FLIGHT LEVEL 50 60 60 80 90 55 100 Relative Humidity (%) 6i 90 100 i 50 45 4 0 200 \ SUBSONIC JET Yes' \Possible No FLIGHT LEVEL 250 35 3\00 \ r-U.S. STANDARD 3001" - ATMOSPHERE, 30 1962 425 50o0 20 800 15 700 -O 800 900 1000 0 -70 -60 -50 -40 -30 -20 Temperature (0C) Figure 14. Relative humidity required for jet aircraft contrail formation as a function of pressure and temperature of the environment (from Appleman, 1966).

75 LATITUDE JAN APR JUL OCT 7004 60 5"0 40 30... 2 0 NORTHERN HEMISPHERE I0 _________ EQUATOR ~ ~ ~ ~..~.~.......0 SOUTHERN HEMISPHERE -30 30~~~~~~~~~~~~~~~~~~ - 40.... 50 60 T00 S~ JAN APR JUL Ocr KEY 50mb 10 km. Figure 15. Latitudes of greater-than-50% contrail probability by season at 50 mb and at 10 km (from Appleman, 1966).

76 corltrails (lasting over five minutes) were found to occur 25 percent of the time, while nonpersistent contrails (lasting under five minutes) occurred in 38 percent of the observations. Assuming the average duration of the persistent contrails to be one-half hour and their average width to be one kin, the total area of contrails formed by a single jet flying continuously at a cruise speed of 900 km/hr would be approximately 450 km2 if persistent contrails were always formed. Since they are actually formed only 25 percent of the time, the average contrail cover produced by one jet is about 100 km2. United States airlines currently fly an average distanlce of about 240 x 10 Km per month in their domestic operatIons. This is about 350 times the distance that could be traveled by a single jet flying continuously at cruise speed. Therefore, the average contrail cover generated over the United States is about 35,000 km2, or about 0.5 percent of the land area. Locally, near airports and heavily traveled air corridors, the average contrail cover is probably 10 times larger than this, or about five percent. This figure is also supported by recent observations (SMIC, 1971). 5.2. ASSESSMENT OF CLIMATIC IMPACT According to the calculations presented in Chapter 4 a continuous contrail cover would cause a decrease in the surface temperature of about 2 OC over dry land at mid

77 latitudes, or a decrease in the evaporation rate of 0.14 nmm per day over wet soil or water. In regions of heavy air traffic the average contrail cover is probably about five percent over the long term. In these regions, therefore, the average surface temperature of dry areas is probably about 0.1 ~C lower than it would be if there were no jet flights there. Alternately stated, it is likely that on one day out of ten there will be generated a contrail layer of sufficient density to lower the surface temperature by 1 OC, or reduce the water evaporation by 0.2 mm. These conclusions are very tentative, since not enough is known about the dynamics of the atmosphere to confidently predict tile effects of arbitrary perturbations. For example, it is possible that a reduction in the evaporation rate due to an increased cirrus cloudiness could result in a decrease in the amount of lower cloudiness, which would negate the effects of the initial perturbation. On the other hand, there may also be positive feedback mechanisms which tend to magnify the initial perturbation. If contrails were generated over the arctic at the same rate that they are currently formed over the United States, there could be a significant and cumulative growth of the arctic ice cover over a period of a few years. This might have serious implications for the global climate as well, if Budyko's theory is correct. Contrail formation is at least as noticeable a phenomenon as the injection of dust into the stratosphere by volcanoes, and the connection

78 between volcanic activity and climatic change has already been demonstrated. 5.3. SUGGESTIONS FOR FURTHER RESEARCH More research is needed in establishing the radiative properties of contrails. This could be accomplished by means of radiometric studies similar to. those of Kuhn (1970), or by means of further investigations into the size, shape, and concentration of ice particles in contrails. This report has largely ignored atmospheric effects, since very little can be said about them until more is known about radiative properties of contrails. However, even if these properties are known accurately, the atmospheric effects will not be predictable without further research in atmospheric dynamics and surface-atmosphere interactions. One difficulty in studies of climatic impact is that the results are often difficult to confirm by observation because there are too many uncontrollable variables. It is possible, however, that an observational program could be devised to at least determine upper limits to some of the effects mentioned in this report. In connection with this, an index combining the area of observed contrails with their thickness and density would have to be developed. The feasibility of establishing a program for monitoring the extent of contrail coverage using such an index should also be looked into.

79 Additional data on the water-vapor content and temperature at cruise altitudes would be useful for predicting and possibly reducing contrail formation by means of routing changes. Such data is especially needed in arctic regions, in order to assess the impact of increasing jet traffic there. Devices for the suppression of contrails have been developed for military use; if sufficient evidence of the adverse effects of contrails is found, these could be installed on commercial planes, although they would be likely to have environmental effects of their own.

APPENDIX A FREEZING TIME AND GROWTH RATE OF ICE PARTICLES The dependence of the rate of freezing of water droplets on their temperature and volume was empirically studied by Bigg (1953). His results were expressed in the following form: Rn(l-P) = -Ar3t(exp AT - 1) (A.1) where P is the probability that a droplet will have frozen after t seconds, r is the radius of the droplet, and AT is the temperature (in degrees Centigrade) below the freezing point. A is a constant depending on the concentration of particles which. act as freezing nuclei in the water. Equation (A.1) may be rewritten in the form P = 1 - e kt (A.2) where k = Ar3(exp AT - 1). (A.3) The time constant, or freezing time, T=l/k is inversely proportional to the number of freezing nuclei, and therefore to the volume of the droplet. Thus the larger particles in any distribution tend to freeze first, causing the smaller particles to disappear by the mechanism discussed in Section 2.1. 80

81 Once the droplets have frozen, they grow by diffusion of water vapor from the environment to the surface of the droplet. The classical expression for the rate of deposition of water onto the particle is dm 4= (S-l) (A. 4) dt u+v where m is the mass of the particle, r is its radius, S is the saturation ratio of water vapor with respect to ice in the environment, and u and v are constants depending on the temperature, the thermal conductivity of air, and the diffusivity of water vapor in air. Since the mass of the particle is proportional to the cube of the radius, Equation (A.4) implies that dr = k (S-1) (A.5) dt r where k is a constant for given environmental conditions. Thus, the smaller particles in any distribution grow faster than the larger ones, and the width of the distribution decreases while the average radius increases.

APPENDIX B EDDINGTON AND TWO-STREAM APPROXIMATIONS The hemispheric flux in the Eddington approximation is defined as - 1 F+(T) = 2 [2Io (T) ~ I (T)1 (B.1) 31 where I~ and I1 are given by Equations (3.17) of the text. Substituting these equations into the above definition yields 2 2t Tt F+(T) = rA1(1+3a)eT + rA2(3a)e 2) The boundary condition F_(0) = 0 gives the following relation: Al(1-3a) + A2(1+ 3a) = o (B 3) The solutions of this equation may be written as A = c(2 a+l) A = C(2a-) (B.4) where c is an arbitrary constant. The flux may therefore be written as F+(T) = rc(l+b)2 e'l - rrc(l-b)2 e 82

83 and F (T) = Tc(l-b2)(e - e-') (B.6) 2= - where b = 3a. The incident flux is F+(T1), so the reflectivity is F(T1) (1-b )(e - e ) R = _F+(i (l+b)aeT' (B. 7) + 1 1+b ) 2e - (1-b) e-T and the transmissivity is F+(0) T E ) 4b _ (B.8) F - T- I + 1 (l+b)2e - (l-b) e In the two-stream approximation, the hemispheric flux is defined as F_+(T) - 2rplI+(T) (B.9) where I+(T) are given by the following expressions: I+ = Al(l+b)e + A2 (l-b)e I = Al(1-b)eT' + A2(l+b)erT (B.10) (c.f. Equations (3.29) of the text). The boundary condition F (O) = O then yields the relation Al(1-b) + A2(l+b) = O (B.11)

84 which has the solutions (for arbitrary c): A c(b+l), A c(b-l) (B.12) Substituting these into Equations (B.10) gives I+ = c(l+b)2 e - c(l-b)2e-T' I = c(l-b2)(e'- e (B.13) The reflectivity is therefore again T V T1 -T1 F (T)' FR(__1) (1 -b2) (e - e ) Ta -T, (B.14) P+(T1) (l+b)2e 1 (1-b)2 e and the transmissivity is F+(O) _4b T =, b (B.15) F+(1l) (l+b)2eTl - (1-b)2 e 1 These are the same expressions as obtained for the Eddington approximation. Note, however, that the value of b is different for the two cases, and that b can be larger than 1 in the Eddington approximation but not in the two-stream approximation. In the case =o 1, both b and T1 are equal to zero and the expressions for R and T are indeterminate. For very small b and T1, these expressions reduce to

85 T R = 1 (B.16) 2b 1 2b + T 1 +,B1 1 2b + 1 1 + T1 where 3 = (1 - 3) in the Eddington approximation, and = (1 - ) in the two-stream approximation. 2 3 The reflectivity and transmissivity for the case w = 1 can be evaluated from these equations.

APPENDIX C DERIVATION OF ALBEDO EQUATION In Section 3.5 an approximate solution of the transfer equation was obtained by splitting the radiation field into two components. The first component consists of the radiation which has been scattered only once inside the layer. This component is given by the expression E 1 - P -1 + 1) (TO p(, 0-+ 0 1-exp (~ C (C.1) 0 for i < 0 where T = whre T1 for k > 0 The second component consists of the radiation which has been scattered two or more times. It may be written as I m(T) = E I n)(T,>) (C.2) n=2 From Equation (3.50), the transfer equation for the multiplyscattered radiation is

87 1 whe re IM( P) PP dP f -1 Using the Legendre expansion of the phase function (3.9), this may be written as =o whe re -1 The Eddington approximation is applied to the multiplyscattered radiation field by assuming Im(T,,) to be of the form Im(T,P) = Io(T) + VI1(T) (C.7) and taking the first two moments of Equation (C.3). The result is the following set of equations:

88 2 dI 3 dT 2 ( 1- 0) o to o() 2 dIo 2 W1 1 2d(1- )I, Wa(T) (C.8) 3 dT =31 3-1 - 3 11(T) In the case w = 1 these equations reduce to the following: dI d 1 2 o - -2ao(T) dI dT ( l-g)I1- 2gl(T) (C.9) where g = c1/3. By simple integration, the solutions are I1(T) = (0) - o ()dT (C.10) T T I(T) = I (0) + (l-g) 1(0)T- 32(1-g) dT' rj (T)d T - g ac 1(T)dT (C.ll) The boundary conditions (zero incident diffuse flux) for the multiply-scattered radiation are

89 2 I (0) _ 2- Il(O) = 0 Io( T1) + 2 I (T1) =0 (C. 12) Substituting Equations (C.10) and (C.11) into the second boundary conditions yields Io(O) + (l-g)IL(O)T1 - 3(1-g)f dT f % (T)dT 0 0 3 g l(T)dT + 2 I (0) - a (-T)dT = 0 (C.13) I1(0) may be eliminated from this equation using the first boundary condition, and the integrations over T can be carried out using Equations (C.1) and (C.6). The result is an expression for Io(0) involving an integral over p of a known function. The albedo Equation (3.54) is finally obtained by substituting this result into the definition of the albedo (3.53) which, from the first boundary condition, may be written as 1) 2 r)21 (O0) A(T1 0o) = E 03)dp +-oE (C.14) 0 PO -1le

APPENDIX D TEMPERATURE OF CONTRAIL LAYER In calculating the effects of a contrail layer on the terrestrial radiation, it was assumed that the layer emits radiation equally in both directions and that the total energy emitted is equal to the energy absorbed. This is equivalent to assuming that the layer is at the uniform temperature Te, given by 2ET e4 = aLt (D.1) where c is the emissivity of the layer, a is the absorptivity, a is the Stefan-Boltzmann constant, and Lt is the flux of long-wave radiation incident on the bottom of the layer. Kirchhoff's law states that the emissivity is equal to the absorptivity, so the radiative equilibrium temperature Te is independent of the cloud thickness. For the global-average value of Lt =.324 cal cm-2min(London, 1957) the radiative equilibrium temperature Te is about 2110K. The standard atmospheric temperature at the tropopause is about T = 2170K. The actual temperature of a the layer would be somewhere between these two limits, closer to Te or to Ta depending on whether radiative energy transfer is more or less important than sensible heat transfer. Since the difference between Te and Ta is small, and the effects at the surface of the earth are not sensitive to 90

91 this difference, it was decided to use the simpler assumption of radiative equilibrium. The assumption that the temperature inside the layer is constant is justified by the fact that the layer is optically thin, so radiation is absorbed and emitted nearly uniformly throughout the layer. In thicker clouds the radiation is absorbed and emitted near the surface of the cloud, so there is a strong heating effect at the bottom and a strong cooling effect at the top of the cloud. In this case the temperature of the cloud would not be uniform, and there would be significantly less radiation emitted from the top than from the bottom of the cloud. However, for transmissivities on the order of 90 percent, it seems justifiable to assume a uniform heating and cooling rate throughout the layer, and therefore a uniform temperature.

REFERENCES Appleman, H. S., "The formation of exhaust condensation trails by jet aircraft," Bull. Am. Meteorol. Soc., 34, 14-20, 1953. Appleman, H. S.., "Effect of supersonic aircraft on cirrus formation and climate," AMS/AIAA Paper No. 66-369, 1966. Beckwith, W. B., "Future patterns of aircraft operations and fuel burnouts with remarks on contrail formation over the United States, "International Conference on Aerospace and Aeronautical Meteorology, May 22-26, Washington, D.C., 1972. Bigg, E. K., "The supercooling of water," Proc. Phys. Soc. B, 66, 688-694, 1953. Budyko, M. I., "The future climate," Trans. Am. Geophys. Union, 53, 868-874, 1972. Chandrasekhar, S., "On the radiative equilibrium of a stellar atmosphere, XX," Astrophys. J., 106, 145-151, 1947. Chandrasekhar, S., Radiative Transfer, Dover Publications, New York, 1960. Dingle, A. N., "Man-made climatic changes: Seeding by contrails," Science, 173, 461, 1971. Flohn, H., Climate and Weather, McGraw-Hill, New York, 1969. Hansen, J. E., "Exact and approximate solutions for multiple scattering by cloudy hazy planetary atmospheres," J. Atmos. Sci., 26, 478-487, 1969. Hansen, J. E. and J. B. Pollack, "Near-infrared light scattering by terrestrial clouds," J. Atmos. Sci., 27, 265281, 1970. Henyey, L. C. and J. L. Greenstein, "Diffuse radiation in the galaxy," Astrophys. J., 93, 70-83, 1941. Hodkinson, J. R., "Light scattering and extinction by irregular particles larger than the wavelength," Electromagnetic Scattering, M. Kerker, ed., 1963. 92

93 Irvine, W. M., "Light scattering by spherical par'ticles: Radiation pressure, asymmetry factor, and extinction cross section," J. Opt. Soc. Am., 55, 16-21, 1965a. Irvine, W. M., "Multiple scattering by large particles," Astrophys. J., 142, 1563-1575, 1965b. Irvine, W. M., "Multiple scattering by large particles II. Optically thick layers," Astrophys. J., 152, 823-834, 1968. Irvine, W. M. and J. B. Pollack, "Infrared optical properties of water and ice spheres," ICARUS, 8, 324-360, 1968. Jacobs, J. D., "Aircraft contrail effects on the surface radiation budget in an arctic region," Bull. Am. Meteorol. Soc., 52_, 1101-1102, 1971. Kuhn, P. M., "Airborne observations of contrail effects on the thermal radiation budget," J. Atmos. Sci., 27, 937942, 1970. Huffman, P. J., "Size distribution of ice fog particles," M.S. Thesis, University of Alaska, College, Alaska, 1968. Huffman, P. J. and T. Ohtake, "Formation and growth of ice fog particles at Fairbanks, Alaska," J. Geophys. Res., 76, 657-665, 1971. van de Hulst, H. C., Light Scattering B Small Particles, Wiley & Sons, Inc., New York, 1957. van de Hulst, H. C., "Multiple scattering in planetary atmospheres," J. Quant. Spect. Rad. Trans., 11, 785-795, 1971. London, J. "A study of the atmospheric heat balance," Final Report, Contract AF19(122)-165, Department of Meteorology and Oceanography, New York University, 99 pp., 1957. Murcray, W. B., "On the possibility of weather modification by aircraft contrails," Mon. Wea. Rev., 98, 745-748, 1970. Nicodemus, M. L. and J. D. McQuigg, "A simulation model for studying possible modification of surface temperature," J. Appl. Meteorol., 8, 199-204, 1969. Pilie, R. J. and J. and J. E. Jiusto, "A laboratory study of contrails," J. Meteorol., 15, 149-154, 1958. Piotrowski, S., "Asympotic case of the diffusion of light through an optically thick scattering layer," Acta Astron., 6, 61-73, 1956.

94 Reinking, R., "Insolation reduction by contrails," Wea., 23, 171-173, 1968. Sagan, C. and J. B. Pollack, "Anisotropic nonconservative scattering and the clouds of Venus," J. Geophys. Rev., 72, 469-477, 1967. Shettle, E. P. and J. A. Weinman, "The transfer of solar irradiance through inhomogeneous turbid atmospheres evaluated by Eddington's approximation," J. Atmos. Sci., 27, 1048-1055, 1970. Study of Critical Environmental Problems, Man's Impact on the Global Environment, MIT Press, Cambridge, Massachusetts, 1970. Study of Man's Impact on Climate, Inadvertent Climate Modification, MIT Press, Cambridge, Massachusetts, 1971. Van Valin, C. C., "Effects on the stratosphere of SST operation," American Geophysical Union 52nd Annual Meeting, April 12-16, Washington, D.C., 1971.

3 9015 125