2.1 Radiation



2.1.1 Ultraviolet (UV) and visible radiative region (117 - 730 nm)

The 2-D model has a total of 170 wavelength bands: a resolution of 500 cm-1 wavenumber between 116.3 and 307.7 nm, a resolution of 5 nm in wavelength between 307.7 and 655 nm, and a resolution of 10 nm in wavelength between 655 and 735 nm (see Table 1 ). At each wavelength, the actinic flux is obtained by solving the equation of radiative transfer with the delta-Eddington method (Shettle and Wienman 1970, Joseph et al.1976), which is a two-stream method (Meador and Weavor 1980; Toon et al., 1989).

The general equation for absorption and scattering of solar radiation in a plane-parallel atmosphere for each wavelength interval can be written as:

eq11

where I is the mean radiance, tau is the optical depth, omega is the single scattering albedo (ratio of scattering extinction and total extinction), p is the scattering phase function, µ is the cosine of the zenith angle, µ0 is the cosine of the zenith angle for the direct solar beam, and I0 is the mean radiance at the upper boundary of the model. The scattering phase function p is a non-dimensional quantity that describes the probability of angular distribution of the scattered energy. The first moment of the phase function is the asymmetry factor g, which gives the overall directionality of the phase function. The first term in the right-hand side of Equation (11) represents attenuation by absorption, the second term represents the diffusive radiation increase by multiple scattering, and the last term is the contribution from single-scattering of the direct solar beam.

Before obtaining a solution of Equation (11) with a two-stream method, delta scaling of the optical paramaters for cases of highly anisotropic phase functions is necessary to improve the accuracy of the solution, based on Joseph et al. (1976):

eq12

This is done to reduce the anisotropy caused by the strong forward-scattering peak characteristic of larger particles (aerosols and clouds).

For the two-stream method, the diffuse radiance is divided into up-welling and down-welling components, I+ and I-, respectively. This produces a pair of integral radiance equations from Eq. (11) (for detail derivation, see Meador and Weaver 1980):

eq13

where coefficients gammai's are determined by the approximation used. For the Eddington approximation, I+ and I- are assumed to have a simple angular distribution given as

eq14

This yields coefficients of the radiance equations as:

eq14a

The coefficients I0 and I1 are evaluated by solving the radiative equation for I+ and I- subject to the boundary conditions that no diffuse radiation is incident at the top of the atmosphere, while at the bottom of the atmosphere (the Earth's surface), radiation is reflected isotropically with a known albedo. Vertical inhomogeneity of the atmosphere is parameterized by subdivison into 120 layers, each being then taken as internally homogeneous with specified vertical optical depth, single scattering albedo, and asymmetry factor. To solve for the set of two-stream equations in this multi-layer inhomogeneous atmosphere, the method described in Toon et al. (1989) is used. The triadiagonal matrix solution of their method provides improved computational speed necessary for the large number of vertical layers in our model.

The solar actinic flux at the top of the atmosphere (q0=4piI0 ) used in the model is shown in Table 1 . The solar flux is in units of photons/cm2/s, with the size of the wavelength bin taken into account (q0(lambda)dellambda). Between the wavelength of 120 to 417 nm, the solar flux is the average of the flux measured by UARS SOLSTICE from Jan. 1st through March 29th, 1995 (representing long-term solar minimum condition) and from Jan 1st through March 31st, 1992 (representing solar maximum condition) (courtesy of G. Rottmann). Also listed in Table 1 is the fractional change of the solar flux from UARS day 169 to 1050 [(max-min)/(min)] representing the long-term solar variability. For wavelengths longer then 417 nm, the solar flux data are adapted from WMO (1986). To account for the ellipticity of the earth's orbit, the solar flux is adjusted with the annual variability : [1.-0.0342*cos(2pi*day/365)], with maximum flux during the winter season and minimum flux during the summer season.

The altitude variation of the optical parameters: extinction optical depth (from absorption and scattering), single scattering albedo, and asymmetry factor, are needed to obtain the solution of the solar radiative flux. Optical depth can be calculated knowing the absorption cross sections (sigma) and number densities ([n]) of absorbing particles. For a clear sky atmosphere, oxygen and ozone are the main absorbers of solar radiation in the middle atmosphere. Their wavelength dependent absorption cross sections are available from laboratory measurements, reported in DeMore et al. (1997) (Note: the particular case of the Schumann-Runge bands will be described later). For the Rayleigh scattering by air molecules, an empirical formula proposed by Nicolet (1984) is used for the cross section:

eq15

The single scattering albedo omega is then calculated from tauscat/(tauscat+tauabs). The asymmetry factor g for Rayleigh scattering is assumed to be zero (isotropic scattering).

To simulate aerosol scattering and absorption in the UV, some additional optical parameters are needed. Different approaches can be taken to obtain the optical parameters of scattering particles, depending on the conditions of interest. For background aerosol levels, Hitchman et al. (1994) compiled from SAGE and SAM measurements the climatology of aerosol extinction coefficient (extaer) at 1 µm over a period of roughly a decade. This aerosol extinction coefficient distribution can be used to calculate the extinction optical depth of aerosols by multiplying it with the thickness of the aerosol layer. To estimate its wavelength dependency, the simplest approach is to scale the extinction coefficient (for example, at 1 µ m) by the inverse of wavelength, which is a fair assumption according to Mie scattering calculations (Pinnick et al. 1980). An alternative is to use a precalculated spectral distribution of the extinction coefficient for background aerosol levels to scale the satellite-measured extinction coefficient at a given wavelength (for example, calculations shown in Fenn et al. 1985). As for the scattering parameters, the sulfate aerosols are known to be very effective scatterers with a single scattering albedo value close to 1, and an asymmetry factor of 0.7, which are approximately independent of wavelength (Michelangeli et al. 1989).

Optical properties for water cloud absorption and scattering in the UV and visible are calculated by the parameterization scheme of Hu and Stamnes (1993) based on Mie theory. The input parameters needed to be specified for this scheme are the liquid water content, effective droplet radius, cloud location, and vertical thickness. Following Kylling (1992), a typical water cloud effective radius of 10 mm and a liquid water content of 0.15 g/m3 are used in the model.

The solar zenith angles chosen for the radiative transfer equation are a function of both latitude and time of the year. As will be mentioned in the description of chemistry, for daytime conditions, the chemical equation is integrated in time with a timestep equal to one quarter of the length of the daylight part of the day. Thus, 2 photolysis calculations are made per day for each molecule subject to photodissociation. The zenith angle Z is calculated from the formula

eq16

where delta is the inclination angle which is a function of day of year only (the latitude at which the sun is directly over head with Z=0), and ha is the hour angle ( ha=0 at noon). The value of the hour angle at sunrise or sunset h1/2day as a function of day and latitude can be determined from Equation (16) by setting Z=90o. Once the hour angle at sunrise (or sunset) is determined, Z is estimated for two specific times: ha=1/8h1/2day, and ha=3/8h1/2day.

To account for the sphericity of the earth in spite of the plane-parallel approximation, the Chapman function replaces the secant of zenith angle in the calculation of the optical depth. The method described by Smith and Smith (1972) for zenith angle less than 90o is used to estimate the Chapman function (Ch(Z)):

eq17

where x=(a+z)/H . For zenith angles equal to or larger then 90o, night conditions are assumed.

a. Photolysis rates.

Once the solar flux at a given level is known from the solution of the radiative equation (Equation 11), the photodissociation coefficient for a given molecule can be calculated by:

eq18

This photolysis rate J (s-1) for a given species i is proportional to the solar actinic flux 4piI, the absorption cross section sigmai (cm2) and the quantum efficiency epsiloni. The photodissociations considered in the model are listed in Table 2 . The sources of reference for the absorption cross sections used in our model are listed in Table 3 . Molecules cross section that exhibit temperature dependencies include: ozone in the 175.8 to 347.5 nm range, N2O, NO2, HNO3, CFC12, CFC-11, CCl4 in the 200-250 nm range, ClONO2, N2O5 in the 200-281.7 nm range, H2O2, H2CO, PAN. For some source gases such a CH3CCl3, CH3Cl, CFC-113, HCFC-22, Ha-1211, Ha-1301, CH3Br, CFC-114, and CFC-115, temperature dependent cross section parameterizations are taken from Gillotay andSimon (1988) since data presented in DeMore et al. (1997) [JPL97] are for 298 K only. Cross checking of the cross sections from Gillotay and Simon (1988) and JPL97 showed that at 298 K, the values are close to identical. For OClO, cross sections are spectrally interpolated directly from Wahner et al. (1987) instead of from JPL97 because its data are take from band peak amplitude which tends to overestimate the cross sections. The quantum efficiency epsilon is equal to unity in most cases, except for O3, H2CO, CF2O, NO3, ClONO2, and HO2NO2. The sources of their epsilon values are listed in Table 3 .

The treatment of O2 absorption in the Schumann-Runge bands (SRB: 175-205 nm) requires special attention because of the cross section's highly variable spectral structure in this wavelength region. The parameterization of the O2 absorption developed by Kockarts (1994) is used in the model. The method uses recent spectroscopic data of O2, and accounts for the temperature dependence of the O2 cross sections. It remains valid even when the slant O2 optical depth becomes very large. The photolysis rates for any molecule (M) other than O2 in the SRB are computed by

eq19

and N is the slant total O2 overhead content, and q0 is the solar flux at the top of atmosphere. In the particular case of molecular oxygen (O2) , the photolysis frequency in the SRB is expressed by :

eq20

The SRB is divided into 17 (iv) wavelength subintervals, and each wavelength subinterval is assigned 12 coefficients (a1~a12, b1 ~ b12) (Table 1 and 2 in Kockarts, 1994).

The photolysis of NO in the middle atmosphere occurs primarily in the delta(0-0) band (190.9 nm) and delta(1-0) band (182.7 nm). Since thes bands coincide with the Schumann-Runge bands of oxygen, the determination of JNO requires a high spectral resolution treatment of the SRB absorption. Minschwaner and Siskind (1993) use the method of opacity distribution function to parameterize JNO with the latest spectroscopic information of both O2 and NO. Their scheme is used in SOCRATES to estimate the photolysis rates of NO. The photolysis frequency of NO coincident with a given SR band of O2 [(5-0), (9-0) and (10-0)] is determined by

eq21

where P(z) = 1 for the delta(1-0) band, and P(z)=1.65e9/(5.1e7+1.65e9+1.5e-9[N2(z)]) for delta(0-0). One given SR band is subdivided into six wavelength regions (i=1,6) and inside each of these six regions by two values (j=1, 2) of NO cross section and weighting factor (WNO). The values of these coefficients are listed in Table 1 of Minschwaner and Siskind (1993).

b. UV solar heating.

The thermal effect of the solar radiation can be calculated from the radiative transfer equation (11) for the mean radiance I. Solar photon absorption by O2 and O3 are the major sources of heat in the stratosphere and mesosphere. Heating in the UV by O2 and O3 knowing I (J/s/cm2) and the absorption cross section sigma of O2 and O3 (expressed in cm2) can be determined by

eq22

where mair is the molecular weight of air (28.8 g/mole), [M] is the air number density (in cm-3), [O2] is the O2 number density (cm-3), and [O3] the ozone number density (cm-3). When the effect of aerosol and cloud scattering is considered in the radiative transfer equation, solar heating due to aerosol (aer) and cloud (wc) heating is determined by:

eq23

in units of (K/s).

In Equation (22), it is assumed that the products of O2 and O3 photolysis recombine instantly as the solar energy absorbed, so that the heat release is equal to the solar energy is absorbed. Although this is generally valid for altitudes below 80 km, the recombination may not be instantaneous at altitudes above 80 km where the atmospheric density is low. Instead, at these higher altitudes, chemical combination may occur only after some time span comparable to the transport time scale. Under this circumstance, the solar energy absorbed less the dissociation energy of photolyzed molecules shall be called 'residual heating' (QR), and the heat released as a result of chemical recombination called 'chemical heating' (QC). Thus, the residual heating is calculated by:

eq24

where zeta is the energy required to break the O2 and O3 bonds. To estimate heating from exothermic chemical reactions, the recombination of odd oxygen species and odd hydrogen species are considered, following Brasseur and Offermann (1986):

		O+O+M ->  O2+M          k1, e1
		O+O2+M -> O3+M          k2, e2
		O+O3 -> 2O2             k3, e3
		H+O3 -> OH+O2           k4, e4
		OH+O -> H+O2            k5, e5
		HO2+O -> OH+O2          k6, e6
		H+O2+M -> HO2+M         k7, e7 
k and e are the reaction rate and the bond breaking energy, respectively. The rate of energy released by these chemical reactions is given by:

eq25

For reaction k4, the heating term is multiplied by a solar enegry efficiency factor kappa , whose value is given below. The unit of QC is converted from J/s/cm2 to K/s, by multiplying QC by NA/([M]maircp), where NA is Avogadro's number, mair is the air molecular mass and cp is the specfic heat at constant pressure.

In and above the mesosphere, solar energy absorbed by molecules in the form of chemical potential energy may be lost by airglow from excited photolysis products or by chemiluminescent emission from product species of exothermic chemical reactions before the energy is able to be converted to thermal energy. Mlynczak and Solomon (1993) evaluated the reduction of the solar heating efficiency by airglow processes for different absorption bands of O2 and O3, and by chemiluminescent loss for exothermic chemical reactions. They showed the heating efficiencies (eta ) to be less than one in the O3 Hartley band, the O2 Schumann-Runge continuum, the O2 Lyman alpha line, and in the case of H+O3 chemical reaction. These efficiencies are expressed by the following expressions:

eta = 1 for all other wavelengths. To obtain the residual solar heating rate, Equation (24) is multiplied by eta at each wavelength region mentioned above. Thus,

eq26

2.1.2 Infrared (IR) radiative regime ( l > 3 µm)

a. Troposphere and stratosphere

The infrared radiative transfer routine used in SOCRATES for altitudes between 0 to 50 km is the long wave band model described in Briegleb (1992). Briegleb's model computes broadband absorptivity and emissivity with a modified Malkmus random band model, corrected for overabsorption of the traditional Malkmus band transmission. The model covers a spectral region of 0 to 3000 cm-1 with a spectral width of 100 cm-1. With this wide spectral region and a relatively high spectral resolution, minor bands of CO2 and O3, and several other radiatively active trace gases can be included in the calculation. In SOCRATES, radiative contributions from CO2, H2O,O3, CH4, N2O, CFC11, CFC12, and aerosols are included.

The optical depths for all gases are computed first, then the broadband transmission is calculated as an exponential of the summed optical depths of all gases. The optical depth for CO2, O3, CH4 and N2O for each 100 cm-1 band takes the form

eq27

where ua is the column absorber amount (g/cm2), Sp and Sd are band parameters [see Appendix of Brigleb (1992)] calculated as a function of line strength and half-line width taken from the HITRAN database (Rothman et. al 1987), beta is the band dependent diffusivity factor to take into account the effect of solid angle integration, alpha is a band dependent parameter which can be estimated from reference calculations, and p* = p+delta , where delta is an asymptotic pressure required to simulate Voigt line effects. For H2O, a different transmission function is given due to the more gradual variation of absorption strength with wavenumber and its wide wavelength absorption span:

eq28

CFC11 and CFC12 weakly absorbs in the H2O window region between 800 to 1200 cm-1, and their optical depths are given as

eq29

where sigma is the absorption cross section in units of cm2/g, and beta(t) the diffusivity factor. The data for the absorption cross sections of the CFCs are from Massie et al. (1991). The diffusivity factor for the CFCs is taken from Ramanathan et al. (1985) as:

eq30

where tau is the total optical depth of all absorbing species. Similarly, the aerosol optical depth is calculated by multiplying its spectral extinction coefficient (provided either through observation or through theoretical calculations) by the aerosol column depth.

Knowing the transmission ( gamma=e-tau), the absorptivity (A) and emissivity (E) of the infrared radiative transfer equation is calculated by:

eq31

B(T) is the spectrally integrated Planck function, and is the transmission for all gases absorbing between pressure p and p' for abs, and p and pt (top boundary pressure) for ems, respectively. The longwave fluxes (upward + and downward -) at pressure p are calculated from the transfer equation by the formulation:

eq32

where subscript g denotes surface ground, subscript s denotes surface air, e is the broadband isotropic surface emissivity, and m is an empirical constant which partially corrects a bias in the surface atmosphere exchange. The upward and downward fluxes are differenced across layers to produce the longwave heating rates

eq33

b. Mesosphere and lower thermosphere (CO2 15 µm)

The infrared radiative heating above the stratosphere is dominated by the CO2 15 µm bands, while contributions from other trace gases are negligible (Kuhn and London, 1969). Above 60 km, where the local thermodynamic equilibrium (LTE) does not necessarily apply due to low air density, a different method for the longwave radiation is required for the CO2 15 µm band. Under non-LTE conditions, the populations of the excited vibrational levels are no longer controlled by collisional processes and the source function of CO2 is no longer equal to the Planck function. For this purpose, Formichev et al. (1993) developed a parameterization scheme for the 15 µm CO2 band radiative heating rate for the middle atmosphere, which is incorporated in SOCRATES to estimate the infrared radiative heating rate above 60 km.

The parameterization scheme of Formichev et al. (1993) for altitudes 15 to 85 km is largely based on the work of Akmaev and Shved (1982), but accounts for the line overlapping in the CO2 band and for the non-LTE effects, in order for the parameterization scheme to be applicable to a wide altitude range, including a non-LTE layer (70-85 km). The cooling rate is calculated by the following formula:

eq34

where j is altitude index set from -8 to 5, x=ln(1000/p(mb)) is the dimensionless height, and

eq34a.

The cooling rate accounts for the radiative exchange between a given level x0 and the other layers of the atmosphere denoted with subscript j. The coefficients aj, bj, and cj in Equation (34) are listed in Tables A, B, and C of Formichev et al. (1993).

For the non-LTE layer between 85 to 115 km, the cooling rate calculation is based on the recurrence formula of Kutepov and Fomichev (1993), which is expressed by

eq35

where XCO2,j is the mixing ratio of CO2 at level j. The quantum probability lambdaj is calculated by

eq36

where XN2,j, XO2,j, and XO,j are the mixing ratios of N2, O2, and O, respectively, at level j, and k0 is the collisional deactivation rate between O and CO2, which is not known with certainty. Currently, k0 is set at 5x107 s-1atm-1. The value of epsilon(xj) is determined by the recurrence formula

eq37

where Aj, and Dj are calculated from

eq38

The parameters dj are listed in Table D of Fromichev et al. (1993).