- Email: [email protected]
Mapping transformations for broadband atmospheric radiation calculations
Vol. 43,No. 3,pp.191-199, 1990 Printed in Great Britain. All rights reserved
0022-4073/90$3.00+ 0.00
J. Quanr. Spectrosc. Radiar. Transfer
Copyright Q 1990Pergamon Press plc
MAPPING TRANSFORMATIONS ATMOSPHERIC RADIATION
FOR BROADBAND CALCULATIONS and LUKE CHEN
ROBERT WEST, DAVID CRISP,
Earth and Space Sciences Division, Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, U.S.A. (Received 21 July 1989)
Abstract-We describe two methods of mapping the line-by-line absorption coefficient (k) spectrum onto a new variable to produce a set of k coefficients which can be used for accurate and efficient broadband radiative transfer calculations in vertically inhomogeneous, non-gray scattering atmospheres. These methods are intended for applications for which the less accurate correlated-k method is inadequate. We tested these methods for model atmospheres containing CO, and HzO, with pressure in the range 0.0549 bar. We obtained results that differed by no more than 1 or 2% from the line-by-line results if the number of k coefficients was about 100 times less than the number required for the line-by-line calculation. These algorithms can be used for gas mixtures and for a variety of pressure-temperature profiles.
1. INTRODUCTION Efficient broadband radiation algorithms are required for a variety of scientific problems, including remote sensing of planetary atmosphere properties, and radiative heating rate calculations for climate models. Both multiple scattering by aerosols and non-gray absorption by gases contribute to the radiation field at many wavelengths. The line-by-line method (hereafter lbl), which computes the radiation field on a spectral grid fine enough. to resolve individual lines, is computationally expensive. Band models cannot be used with multiple scattering, except in limiting circumstances.
Some investigators ‘~6 have adopted methods which sample the distribution of absorption coefficients (k) to reduce the number of calculations needed in these applications. The key idea is that one can group together spectral intervals which have similar absorption coefficients, and calculate the radiation field for the group as a whole rather than for each spectral interval. The development of these k-distribution methods can be traced back to work by Ambartzumian.’ By the early 197Os, a nomenclature had become established in the literature (Arking and Grossman,8 Domoto’), with the absorption coefficient distribution function denoted byf(k). In more recent work (Lacis and Oinas,” Goody and Yung,” Goody et al”), g(k) is used for the cumulative frequency distribution, i.e. g(k) =
‘f(F) s0
dk’.
(1)
The quantity g varies between 0 and 1. Equation (1) relating g to k can be inverted to yield k as a function of g. The absorption coefficient k increases monotonically with g. The quantity g replaces v in the radiation calculations, with dg = dv/Av (Av is the width of the spectral window). The cumulative k distribution can be calculated in at least two ways. Some investigators (e.g., Wang and Shi,13 Lacis and Oinas”) have fitted band-model parameters to k-distributions or vice versa. Goody et al’* (hereafter, GWCC) formed k(g) by sorting the k(v) produced by a lb1 code. We refer to both the Lacis and Oinas technique and that of GWCC as the correlated-k (hereafter c-k) method. It is useful to dispense withf(k), which is not used in the radiation calculations, and to think of the c-k method as a series of spectral mapping transformations. The c-k method performs a mapping for each atmospheric layer j of the frequency v onto variables gj such that k increases monotonically with g for each layer. It is assumed that the k are correlated, which is to say that QsRT 43,3--A
191
ROBERTWEST et al
192
the mapping is the same for all layers, and the individual gj can be replaced by a single variable g. This step is justified under restricted conditions.“*‘2 The accuracy of the c-k method for inhomogeneous atmospheres is limited by the degree to which the k coefficients are correlated over the relevant pressure-temperature regime. The mapping of v onto g as defined above differs for each atmospheric layer except in special cases.“<‘* Figures 1 and 2 illustrate the nature of the problem. In Fig. 2, the CO2 spectrum shown in Fig. 2 was mapped onto a variable g such that k increases monotonically for the bottom layer in our model atmosphere (which will be described in more detail). The same mapping was applied to the k spectrum of other layers in the atmosphere. The fundamental assumption in c-k that the k values are correlated means that a single mapping would produce a monotonically increasing k spectrum in every layer. Figure 2 shows that a mapping which produces a monotonically increasing k spectrum for one layer produces a k spectrum for other layers in which neighboring k values fluctuate by orders of magnitude (just how much depends on the k distribution). The assumption of correlation in the c-k method will introduce error into the radiation calculations. It is impossible to specify, a priori, the accuracy of the c-k approximation, and the only way to discover the magnitude of error is to compare the c-k results with those computed by the lb1 method. Goody et al” found that the errors (defined by the percent difference between c-k and lb1 calculations) were usually of order 1% but in some cases could not be reduced to less than a few tens of percent. Errors that are large are unacceptable for many applications. To address this shortcoming, we have developed two spectral mapping methods which preserve a one-to-one correspondence between the mapped variable, R, and v. These methods allow the user to achieve arbitrary accuracy, relative to the-lb1 result. The number of k coefficients needed to achieve this accuracy depends on the properties of the k distribution and the atmospheric path distribution. 2. TWO
MAPPING
ALGORITHMS
The c-k method and the two new algorithms described here begin with a k spectrum for each layer (e.g., Fig. 1). When calculated numerically, this spectrum is a list of frequency intervals, Vi,i and v2,,, and associated mean k values, k,,j. The subscripts 1 and 2 denote the beginning and ending frequency values for a discrete spectral interval i. The subscribt j refers to the layer number. The size of each interval, 6vi = v2.!- v’,~, is small enough (i.e., a fraction of a line width) to resolve individual lines. In the discussion that follows it is assumed that the user wishes to limit errors, defined as the difference between the broadband result (i.e., solar or thermal fluxes, flux divergences, transmission along a path, etc. integrated over the spectral interval Av) and the lb1 result, to a value smaller than L.
5000
5010
5020 Wavenumber
5c30
5040
5050
(cm -1)
Fig. 1. CO, absorption coefficient spectra for layers 1 (top curve), 4. and 7 of the model atmosphere (pressure and temperature are listed in Table 1). In this figure and those following, the coefficients for layer 4 have been multiplied by IO-“ and those for layer 7 by IO-‘.
-1OL
0
0.2
0.4
0.6
0.8
1.0
1
9
Fig. 2. CO* absorption coefficient spectra of Fig. I sorted in wavelength such that k increases monotonically for layer 7 (bottom curve). The wavelength was mapped onto g, as described in Sec. 2 for Method I.
Mapping transformations
for RT calculations
193
2.1. Method I In Method I we perform a series of mappings of spectral frequency onto a new variable g. Knowledge of the mapping is maintained until the final step in which the mean k value for each sub-interval is calculated. The procedure is defined by a series of steps as follows: (i) we choose a starting layer (e.g., layer A) and map v onto gA as is done for the c-k method, such that k&A) increases monotonically with g for this layer. We call this mapping TA. Following the c-k convention (this is not necessary but reinforces the connection with c-k), gA spans the interval O-l, and 6gi = 6v,/Av. (ii) We apply map TA to the k spectrum for all other layers in the atmosphere. Unless the k distribution is perfectly correlated in the sense assumed by the c-k method, the resulting k(gA) will not increase monotonically with gA for the other layers. The result of this operation for the CO2 spectrum of Fig. 1 is shown in Fig. 2, where the bottom layer was chosen as layer A. The mapping steps are fast compared to the calculation of k(v). We used the INDEXX routineI to perform the mappings. (iii) We establish subdivision boundaries in gA space for subsequent numerical integration based on the largest range of k, Ak,,j, which will provide results in the radiation calculation with errors no larger than L. The quantity Ak,,j can vary with each layer, and with each subinterval (I) of gA . The magnitude of Ak,,j for a given t will depend on the details of the k spectrum and on the path length distribution and gas abundance and cannot be specified here. The user must determine the appropriate step size. We are not aware of any systematic way to determine Ak. In our tests we derived Ak by trial and error. (iv) We examine the range of k within each g, subdivision for every other layer in the atmosphere. If, for every layer, and within each subdivision of g,, the range of k is no larger than Ak,,j defined in the previous step, then go to step (viii). (v) We choose another layer, layer B, and map g, onto a new variable, g,, , such that within each subdivision defined in step (iii), k,(g,,) increases monotonically. This mapping is performed with a sorting algorithm operating on one g, subdivision at a time. The mapping is denoted by TAB, since the k coefficients of layer A define the map TA for the interval 0
0
0.2
0.4
0.6
0.6
1.0
g
Fig. 3. layers
absorption coefficient spectra after sorting k for and 4 according to step (vi) of Method I (see the text).
CO? 7
0
0.2
0.4
0.6
0.8
1.C
g
Fig. 4. CO2 absorption coefficient spectra after sorting k for layers 7, 4, and I according Method I [cf. step (vii) of Method I].
ROBERT WEST et al
194
mean k coefficients, k,,j, for each subdivision and layer in the final g,,,,,, space. This step also requires insight by the user, and depends on which quantities are of interest (transmissions, fluxes, flux divergence, etc.). One way to formulate an expression for k is to require that the transmission through the layer along a mean path be identical to that given by lbl. If fi, the product of the mean path length and absorber amount, is known for each layer, the specified requirement implies -In k=
[i li
16g, exp( - k,C)/Ag
I,
(2)
where 6g, is the size of the g-interval associated with each ki value, and Ag is the size of the subdivision in gABC,., . Equation (2) can be used to define C, by minimizing errors in the desired quantity. Equation (2) reduces to the frequency- (or g-) weighted mean of k in the weak absorption limit. If it is known that the U are small (e.g., if one is only interested in vertical transmission, and the atmospheric layers have small k and small column abundance), the linear average of k is appropriate. In practice a distribution of paths and absorber amounts must be considered, and the width of the path distribution and the value E determine the size of Ak used in step (iii). For numerical tests, we used Eq. (2) by calculating 6 along a slant path with emission angle equal to 45”. The 90-term case developed in Figs. l-4 produced a sequence of & values shown in Fig. 5. 2.2. Method II Our second method uses a simple binning scheme to group together all 6v, within Av that have similar k values at all levels of the model atmosphere. The binning proceeds as follows: (i) for each layer j, we define a value for the maximum acceptable range Ak,,, of k to be combined into a given bin, 1. This quantity can vary as a function of the absorption coefficient. As for Method I, the user must investigate the appropriate Ak,,j by trial and error. For numerical tests, we choose Ak,,, such that 1k - &I < 0.5 6 ,.,. (ii) W e use the k values in the first 6v interval, k,,,, to initialize the width and mean absorption coefficients for the first bin: 6g, = 6v, /Av; k,,, = k,,,. (iii) We examine the k values for the next frequency interval, 6v,, and compare them to the k/,, in existing bins, If, for any layer, j, and for every bin, I, (k,., - I& > 0.5Ak,,j, we create a new bin, as in step (ii). If, for some bin, I, and for every layer, j, 1k,,j - k,,, 1< 0.5Ak ,,,, we add this interval to bin 1, set 6g, = 6g, + 6v,/Av and redefine the mean absorption coefficient for bin 1,j by adding in the contribution from k,.,. As in step (viii) of Method I, the user must investigate the appropriate averaging technique. We calculated a linear (6g-weighted) average of k for the numerical tests described in the next section. (iv) We repeat step (iii) until all k,, are gathered into bins. This simple scheme is extremely fast because it only involves differences and comparisons. It also requires very little storage because it operates on only one monochromatic spectral interval at a time. The total storage is L(2 + J) + J, where L is the number of bins and Jis the number of model atmosphere layers. This algorithm generated 123 terms for the spectrum shown in Fig. 1, with Ak = k. These are shown in Fig. 6 for three layers. 1
-101
0
0.2
0.4
0.6
0.6
1.0
9
Fig. 5. CO? mean absorption coefficients for layers 1,4, and 7, calculated from Eq. (2) by using the mapping produced by Method I [cf. step (viii) of Method I].
I
-101
0
0.2
0.4
0.6
0.6
1.0
9
Fig. 6. CO, mean absorption coefficients for layers 1,4, and 7 from the binning produced by Method II.
I
Mapping transformations
for RT calculations
3. NUMERICAL
195
TESTS
We performed a series of calculations to test the accuracy of our methods. To compare our methods with the c-k technique, we calculated radiation fields for atmospheric models used by GWCC. These calculations use COZ and HZ0 opacities in the 5000- to 5050-cm-’ region (Rothman et al15). The minimum step size in frequency is 0.25 times the Doppler width at the lowest temperature (layer 1). In the calculation of thermal emission we follow GWCC and set the Planck function to its mean value in the 50-cm-’ interval to isolate errors intrinsic to the sorting algorithms. We begin with calculations for a pure CO* atmosphere (Case A of GWCC), and follow this with calculations for mixtures of CO, and H,O (beginning with Case 0 of GWCC). Temperature, pressure, and gas abundances for each of the models are given in Table 1. Layers 2,4, and 6 contain isotropically scattering aerosols. Aerosol extinction optical depth is 10 in each layer. Aerosol single scattering albedos are 1.0 except for Case A, where they are 0.9999, 0.999, and 0.995 for layers 2, 4, and 6, respectively. Solar and thermal radiation fields were calculated with doubling and adding codes ” having 22 streams. We examined the same quantities (up and down solar and thermal fluxes, net fluxes, and flux divergences) that were discussed by GWCC for the c-k method. All methods, including c-k, produced results which agreed with the lb1 calculations to within about 1% for most quantities except the thermal flux divergence for layers 1 and 2, where the error in the c-k method was as large as 30%. We show in this paper only thermal flux divergence errors. 3.1. Single gas Figure 7 shows thermal flux divergences for Case A. Errors arising from the various methods are listed in Table 3. The cases in Tables 2 and 3 differ in the number of intervals for which E is computed. The c-k calculation examined by GWCC used 500 intervals, but accuracy does not improve for more than about 20 intervals. Two sources contribute to the errors. First, if the Ak are too large in any layer, the range of k will be too large to accommodate a wide range of path lengths. That effect imposes a minimum number of sub-intervals for each layer, of order 10 for Case A if the desired error is a few percent or less. It is the dominant source of error in Method I for layer 7 (Table 3). A second source of error is related to the range of k values for neighboring g values in the plot of k vs g in Fig. 2. This source is responsible for the large errors in layers 1 and 2 in the c-k method. It requires that several k coefficients for the top layer be chosen for a given k value in the bottom layer, and vice versa. It introduces a multiplicative factor in the number of terms each time a new mapping must be carried out according to step (v) of Method I. The same effect influences the number of terms in Method II. Case A required sorting on the bottom, middle, and top layers to achieve errors of less than a few percent using Method I. The case with the minimum number of intervals has 10 divisions for the bottom layer, times 3 for the middle, times 3 for the top, or 90 total. 1”~‘1”“J~!,‘1’1~~
Table 1. Model atmosphere layer parameters.
+
2x10-1' iti 5g\
0 ______~_____i___‘__~____._______.__________._____._.___
gg 4 E 2xd'
The
parameterfis the mass mixing ratio of CO,. The mixing ratio of H,O is 1-J The surface is a Lambert reflector with reflectivity and emissivity 0.5.
+ + + I,.,.IIIII1ll~.ll.IIIIIIII~IIIIIII~ 012345678
Layer Fig. 7. Thermal flux divergences for Case A.
ROBERTWESTet al
196 Table 2. Single gas test cases.
ays
Whod
-z-
c-k
90
3
135
3
150
5
160
3
225
3
250
5
300
3
123
-
?-
Akii
1
LaVEr
hlethoi
-K
c-k
2
3
T-
T-
-.-L-
4 __ < 0.10
-ii-
6
0.10
-K2
7 0.17
10
90
-0.21 -0.76 -0.51
0.39
0.63
15
135
-0.55 -0.76 -0.49 < 0.10
0.26
0.59
1.6
10
150
-0.46 -0.26 -0.20 < 0.10
0.37
0.61
2.4
20
180
-0.69 -0.69 -0.47 < 0.10
0.16
0.36
1.2
15
225
0.63
-0.39 -0.24 < 0.10
0.26
0.56
1.6
10
250
0.49
-0.17 -0.14
0.32
0.56
2.2
20
300
0.49
-0.41 -0.36 < 0.10
0.14
0.34
1.1
1.72
-0.13 -0.69 -0.21
-1.0
-0.49
-0.63
1.0
-
Table 3. Thermal flux divergence errors, single gas.
The total number of terms in the k distribution are listed in the second column, the number of subdivisions for each of the top, middle, and bottom layers for Method I are listed in the next three columns. The bin size criterion, Ak/k, for Method II is listed in the final column.
123
-
2.0
-
The error is 100 x (sorting algorithm - lbl)/lbl. The c-k results were obtained by GWCC.
The several cases for Method I listed in Tables 2 and 3 were designed to examine the errors as the number of terms increases. The error for layer 7 decreases as the number of divisions for that layer increases from 10 to 20. Errors in the other layers are all < 1% and do not show any systematic relationship to the number of subdivisions. It is difficult to know the sign and magnitude of the various contributions to error once the error becomes 5 1%. It is also difficult to justify the computational expense of achieving higher accuracy because errors or uncertainties in gas absorption coefficients and properties of the atmosphere impose greater limitations. Both Methods I and II produce comparable errors. 3.2. Gas mixtures We tested two approaches to the problem of gas mixtures. The first is to calculate the k spectrum for the gas mixture and treat it as a single-gas spectrum using the methods described in Sets. 2 and 3. We used gas abundances for Case 0 (cf. Table 1) to test this method. The CO* mixing ratio varies from 0.2 in layer 1 to 0.95 in layer 7. Calculations using our sorting algorithms produced results with errors of a few percent or less, listed in Table 4. This method requires that a new mapping be produced each time the mixing ratios change. We devised a second, more flexible, method for gas mixtures. It is a perturbation approach which begins by calculating the k spectrum for a model with gas mixing ratios having mean values in each layer. The sorting algorithm is then executed to map the combined-gas spectrum according to the procedures described in Sets. 2 or 3. Instead of forming L for the combined gas, we calculate lists of E for each gas separately and store them. When models with different mixing ratios are constructed, 6 values for each gas are combined as in the lb1 procedure. We tested this method with four models. The first, Case T is a pure CO, atmosphere (same as Case A). The second, Case U, is a pure Hz0 atmosphere. In the third, Case V, the mixing ratio of the minor constituent in each layer is increased by 30% over that for Case 0. The fourth, Case N, (this case was studied by GWCC) contains equal amounts of both gases in all layers. The lists of E for the separate gases were calculated using the mixing ratio profile of Case 0. The errors for each of these cases are listed in Table 5. Errors remain small when the mixing ratio profile differs by 30% from the profile used to compute E, as the results for Case V demonstrate. They can become quite large for large departures (Case T) because the sorting was performed on a combined-gas spectrum and is not optimum for a single gas. Errors for Case N, in the intermediate regime, are a little larger than those for Case V. The errors produced by the perturbation approach described above can become significant when the mixing ratio profile departs by more than a few tens of percent from the value used to calculate the E coefficients. In those cases one can use the first approach described in this section, but a more flexible method is desired.
Mapping transformations Table 4. Thermal flux divergence errors, Case 0. I
1
a”ar
for RT calculations
Table 5. Thermal flux divergence errors, Cases T, U, V, and N. C=e
Tems
layer 1
T
1
j
-0.17 -0.61 -0.66
/
-1.1
j
1.3
197
300
-74
1 2 1 3 1 4 1 5 1 6 1 7 ) -29 ) -14 1 -5.0 1-0.16 /co.101 4.2
j
The c-k results were obtained by GWCC.
If the spectra of the two gases are uncorrelated it is possible to generate a set of k coefficients for gas mixtures from the lists of Sg, and E,,jfor each gas. The individual Eare combined by creating a new set of intervals (denoted by g,,*,,). For each g interval in the listing for gas 1, 6g,.,, there are N 6g,,2,, intervals for the gas mixture corresponding to N 6g,,l intervals for gas 2. If the number of 6g,,, intervals is M, the number of 6g,,2,, intervals is M x N. The &j for the combined gas for each of these intervals is the average of the &,,j (for gas 1) and &,j (for gas 2) weighted by the mixing ratios of the gases in layer j, where 1= m x (n - 1) + n, and interval m for gas 1 was combined with interval n for gas 2 to produce interval 1 for the combined gas. The it4 x N intervals are then mapped and regrouped by Methods I or II to produce a much smaller number ( N the larger of M or N) of intervals for use in the radiation calculations. This method requires little storage capacity because only the E for each gas are stored. It is flexible because it is not restricted by the mixing ratio profile. Its accuracy is limited by the degree of correlation of k coefficients for the two gases. We have not tested this method for the mapping algorithms derived here, and we cannot predict how well it will perform in general. 4.
TEMPERATURE
VARIATIONS
It is cumbersome to recalculate the lb1 k spectrum for each new temperature-pressure profile. A much more efficient approach is to calculate k coefficients on a temperature grid and to interpolate each time a new set of k coefficients is needed. Goody et al” demonstrated the method for the c-k technique. Since the number of k values which are stored and used in the interpolation is small, the method is very efficient and does not require large storage capacity. The same technique can be used with the algorithms presented here. The same mapping must be used for each temperature to ensure that the spectral content of each of the bins does not change with temperature. We tested the interpolation method with a CO, hot band (where the temperature sensitivity is strong) in the region 595-600 cm-‘. We compared thermal flux divergences for a temperature profile appropriate for the Martian atmosphere (T varies from 140 K at altitudes above the 5-p bar level to 2 14 K at the 6.4-mbar level). The interpolation used k coefficients calculated at temperature intervals of 30 K. The largest difference in thermal flux divergence between the calculation using interpolated coefficients derived from Method I, and the non-interpolated lb1 calculation was 3%. We consider that difference to be acceptable since the lb1 absorption coefficients are not known to that accuracy. 5. DISCUSSION
We have introduced two spectral mapping algorithms for broadband radiation calculations in planetary atmospheres. It is appropriate to examine their relationship to each other and to another mapping algorithm, the c-k method. The algorithms described here are quite similar. Their goal (shared with the c-k method) is to reduce the number of radiation calculations by grouping together spectral regions with similar values of k. The two methods differ in the order in which the steps are executed. Method I performs a series of mapping transformations based on the k values of the various atmospheric layers, and then combines the k values at the end to form mean k values. Method II begins by defining criteria
198
ROBERT WFST et al
(the A/c,.~)which determine which group a given spectal interval will be mapped onto. Method I retains the full spectral information (list of v and k values) until the final step to compute 5 whereas Method II discards that information (unless it is needed for other mixing ratios or temperature profiles) after each spectral interval has been assigned a bin. The mapping transformations performed by Methods I and II are not necessarily equivalent. By equivalent we mean that the k,,j are formed from the same set of ki,j. If the grouping criteria were defined in the same way for both methods, the resulting radiation calculations would be identical for both. The two methods as described above are not equivalent because the criterion for bin membership in Method II depends on a moving average of k. A spectral interval which satisfied the criterion Ik,,, - &,,,I6 OSAk,,, when it was tested may not satisfy that criterion for the final value of E,,,. The moving average is not an essential feature of Method II. We also use a variant of Method II which has fixed Ak intervals and is equivalent to Method I. The mapping algorithms we introduced here differ fundamentally from the c-k technique whose validity is restricted to the weak absorption limit, or to cases where correlation assumption is valid (i.e., Elsasser models or spectra which obey a scaling law”). When the correlation assumption is strictly valid, all three mapping transformations are equivalent (except for the nonessential difference between Method I and one form of Method II noted above). The mapping algorithms described in the Sets. 2 and 3 can therefore be considered to be more general ones, whose validity is unrestricted apart from the requirement that the source radiation be nearly constant over the interval. Of course, the number of intervals needed to achieve arbitrary accuracy may be large. In the limit Ak -+ 0, the lb1 results are reproduced and the number of g intervals can be as large as the number of v intervals in the original spectrum. More work needs to be done to understand how to make optimum choices for the size of k intervals needed to achieve a desired accuracy. This problem was encountered earlier for the c-k method. It is more complex for these new methods since there are more degrees of freedom. Only the number and distribution of g-intervals was in question in the c-k method. In the new methods one must explore the trade-off between the Ak for one layer vs the Ak for another. We currently investigate these issues by trial and error. We have described two techniques to reduce the number of terms to a minimum number (cf. Sec. 2). Although these techniques are useful, a more systematic approach (e.g., cluster analysis) is desired.
6.
SUMMARY
We have described two methods of mapping the lb1 k spectrum to produce a set of k coefficients which can be used for accurate broadband radiative transfer calculations in non-gray, vertically inhomogeneous scattering atmospheres. These methods are intended for applications where the less accurate c-k method is inadequate. The algorithms described here allow the user to achieve results arbitrarily close to the lb1 results. They reduce the number of radiation calculations by about two orders of magnitude, compared to the lb1 technique, for the cases we tested. These algorithms can be used for mixtures of gases, and for a variety of pressure-temperature profiles. The two methods we introduced are similar, but differ primarily in the way in which they are executed. Method II requires little memory since it operates on only one entry in the lb1 spectrum at a time. Knowledge of the wavelength mapping is not retained, although it could be if desired. Knowledge of the mapping of wavelength onto g is preserved in Method I until the final step, when averages of k are formed. This requires significant memory during the mapping, although the mapping is fast. Techniques for combining gas absorption coefficients and interpolating on temperature grids are identical for both methods. These methods are essentially equivalent to the c-k technique for Elsasser band models and for cases where the absorption coefficient obeys a scaling law (e.g., the strong-line limit). The mapping transformation developed for Method I illustrates why the c-k method fails (Fig. 2). Our numerical tests underscore the need to obtain high-resolution spectroscopic measurements in the laboratory. Much of the high-frequency variation seen in the upper curve in Fig. 2 is produced in the near-wings of lines. Calculations which ignore this effect may produce errors of order several tens of percent.
Mapping transformations
for RT calculations
199
have benefited from a close collaboration with Professor Richard Goody, especially in the early stages of the work described here. The support given by Daniel McCleese was an essential ingredient in the success of this oroiect. This work was nerformed bv the Jet Prouulsion Laboratory, California Institute of Technology, under contract bith the National Aeronautics and Space Administration. The authors were supported under contracTs administered by the NASA Planetary Atmospheres program.
Acknowledgemenrs-We
REFERENCES 1. W. C. Yang, Y. L. Yung, A. A. Lacis, T. MO, and J. E. Hansen, Science 194, 685 (1976).
2. J. E. Hansen, W. C. Wang, and A. A. Lacis, Science 199, 1065 (1978). 3. A. Lacis, J. Hansen, P. Lee, T. Mitchell, and S. Lebedeff, Geophys. Res. Lett. 8, 1935 (1981). 4. J. Hansen, G. Russell, D. Rind, P. Stone, A. Lacis, S. Lebedeff, R. Ruedy, and L. Travis, Mon. Weath. Rev. 111, 609 (1983). 5. J. Hansen, 1. Fung, A. Lacis, D. Rind, S. Lebedeff, R. Ruedy, and G. Russell, J. Geophys. Res. 93,934l (1988). 6. D. Rind, R. Suozzo, N. K. Balachandran, A. Lacis, and G. Russell, J. Atrnos. Sci. 45, 329 (1988). 7. V. Ambartzumian, Publs Obs. Astr. Univ. Leningrad 6, 7 (1936). 8. A. A. Arking and K. Grossman, J. Atmos. Sci. 29, 937 (1972). 9. G. A. Domoto, JQSRT 14, 935 (1974). 10. A. A. Lacis and V. Oinas, submitted to J. Geophys. Res. (1990). Il. R. M. Goody and Y. L. Yung, Atmospheric Radiation: Theoretical Basis, Oxford University Press, Oxford (1989). 12. R. Goody, R. West, L. Chen, and D. Crisp, JQSRT 42, 539 (1989). 13. W.-C. Wang and G.-Y. Shi, JQSRT 39, 387 (1988). 14. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes, Cambridge University Press, Cambridge (1986). 15. L. S. Rothman, R. R. Gamache, A. Goldman, L. R. Borwn, R. A. Toth, H. M. Pickett, R. L. Poynter, J. M. Flaud, C. Camy-Peyret, A. Barbe, N. Husson, C. P. Rinsland, and M. A. H. Smith, Appl. Opt. 26, 4058 (1986).
0022-4073/90$3.00+ 0.00
J. Quanr. Spectrosc. Radiar. Transfer
Copyright Q 1990Pergamon Press plc
MAPPING TRANSFORMATIONS ATMOSPHERIC RADIATION
FOR BROADBAND CALCULATIONS and LUKE CHEN
ROBERT WEST, DAVID CRISP,
Earth and Space Sciences Division, Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, U.S.A. (Received 21 July 1989)
Abstract-We describe two methods of mapping the line-by-line absorption coefficient (k) spectrum onto a new variable to produce a set of k coefficients which can be used for accurate and efficient broadband radiative transfer calculations in vertically inhomogeneous, non-gray scattering atmospheres. These methods are intended for applications for which the less accurate correlated-k method is inadequate. We tested these methods for model atmospheres containing CO, and HzO, with pressure in the range 0.0549 bar. We obtained results that differed by no more than 1 or 2% from the line-by-line results if the number of k coefficients was about 100 times less than the number required for the line-by-line calculation. These algorithms can be used for gas mixtures and for a variety of pressure-temperature profiles.
1. INTRODUCTION Efficient broadband radiation algorithms are required for a variety of scientific problems, including remote sensing of planetary atmosphere properties, and radiative heating rate calculations for climate models. Both multiple scattering by aerosols and non-gray absorption by gases contribute to the radiation field at many wavelengths. The line-by-line method (hereafter lbl), which computes the radiation field on a spectral grid fine enough. to resolve individual lines, is computationally expensive. Band models cannot be used with multiple scattering, except in limiting circumstances.
Some investigators ‘~6 have adopted methods which sample the distribution of absorption coefficients (k) to reduce the number of calculations needed in these applications. The key idea is that one can group together spectral intervals which have similar absorption coefficients, and calculate the radiation field for the group as a whole rather than for each spectral interval. The development of these k-distribution methods can be traced back to work by Ambartzumian.’ By the early 197Os, a nomenclature had become established in the literature (Arking and Grossman,8 Domoto’), with the absorption coefficient distribution function denoted byf(k). In more recent work (Lacis and Oinas,” Goody and Yung,” Goody et al”), g(k) is used for the cumulative frequency distribution, i.e. g(k) =
‘f(F) s0
dk’.
(1)
The quantity g varies between 0 and 1. Equation (1) relating g to k can be inverted to yield k as a function of g. The absorption coefficient k increases monotonically with g. The quantity g replaces v in the radiation calculations, with dg = dv/Av (Av is the width of the spectral window). The cumulative k distribution can be calculated in at least two ways. Some investigators (e.g., Wang and Shi,13 Lacis and Oinas”) have fitted band-model parameters to k-distributions or vice versa. Goody et al’* (hereafter, GWCC) formed k(g) by sorting the k(v) produced by a lb1 code. We refer to both the Lacis and Oinas technique and that of GWCC as the correlated-k (hereafter c-k) method. It is useful to dispense withf(k), which is not used in the radiation calculations, and to think of the c-k method as a series of spectral mapping transformations. The c-k method performs a mapping for each atmospheric layer j of the frequency v onto variables gj such that k increases monotonically with g for each layer. It is assumed that the k are correlated, which is to say that QsRT 43,3--A
191
ROBERTWEST et al
192
the mapping is the same for all layers, and the individual gj can be replaced by a single variable g. This step is justified under restricted conditions.“*‘2 The accuracy of the c-k method for inhomogeneous atmospheres is limited by the degree to which the k coefficients are correlated over the relevant pressure-temperature regime. The mapping of v onto g as defined above differs for each atmospheric layer except in special cases.“<‘* Figures 1 and 2 illustrate the nature of the problem. In Fig. 2, the CO2 spectrum shown in Fig. 2 was mapped onto a variable g such that k increases monotonically for the bottom layer in our model atmosphere (which will be described in more detail). The same mapping was applied to the k spectrum of other layers in the atmosphere. The fundamental assumption in c-k that the k values are correlated means that a single mapping would produce a monotonically increasing k spectrum in every layer. Figure 2 shows that a mapping which produces a monotonically increasing k spectrum for one layer produces a k spectrum for other layers in which neighboring k values fluctuate by orders of magnitude (just how much depends on the k distribution). The assumption of correlation in the c-k method will introduce error into the radiation calculations. It is impossible to specify, a priori, the accuracy of the c-k approximation, and the only way to discover the magnitude of error is to compare the c-k results with those computed by the lb1 method. Goody et al” found that the errors (defined by the percent difference between c-k and lb1 calculations) were usually of order 1% but in some cases could not be reduced to less than a few tens of percent. Errors that are large are unacceptable for many applications. To address this shortcoming, we have developed two spectral mapping methods which preserve a one-to-one correspondence between the mapped variable, R, and v. These methods allow the user to achieve arbitrary accuracy, relative to the-lb1 result. The number of k coefficients needed to achieve this accuracy depends on the properties of the k distribution and the atmospheric path distribution. 2. TWO
MAPPING
ALGORITHMS
The c-k method and the two new algorithms described here begin with a k spectrum for each layer (e.g., Fig. 1). When calculated numerically, this spectrum is a list of frequency intervals, Vi,i and v2,,, and associated mean k values, k,,j. The subscripts 1 and 2 denote the beginning and ending frequency values for a discrete spectral interval i. The subscribt j refers to the layer number. The size of each interval, 6vi = v2.!- v’,~, is small enough (i.e., a fraction of a line width) to resolve individual lines. In the discussion that follows it is assumed that the user wishes to limit errors, defined as the difference between the broadband result (i.e., solar or thermal fluxes, flux divergences, transmission along a path, etc. integrated over the spectral interval Av) and the lb1 result, to a value smaller than L.
5000
5010
5020 Wavenumber
5c30
5040
5050
(cm -1)
Fig. 1. CO, absorption coefficient spectra for layers 1 (top curve), 4. and 7 of the model atmosphere (pressure and temperature are listed in Table 1). In this figure and those following, the coefficients for layer 4 have been multiplied by IO-“ and those for layer 7 by IO-‘.
-1OL
0
0.2
0.4
0.6
0.8
1.0
1
9
Fig. 2. CO* absorption coefficient spectra of Fig. I sorted in wavelength such that k increases monotonically for layer 7 (bottom curve). The wavelength was mapped onto g, as described in Sec. 2 for Method I.
Mapping transformations
for RT calculations
193
2.1. Method I In Method I we perform a series of mappings of spectral frequency onto a new variable g. Knowledge of the mapping is maintained until the final step in which the mean k value for each sub-interval is calculated. The procedure is defined by a series of steps as follows: (i) we choose a starting layer (e.g., layer A) and map v onto gA as is done for the c-k method, such that k&A) increases monotonically with g for this layer. We call this mapping TA. Following the c-k convention (this is not necessary but reinforces the connection with c-k), gA spans the interval O-l, and 6gi = 6v,/Av. (ii) We apply map TA to the k spectrum for all other layers in the atmosphere. Unless the k distribution is perfectly correlated in the sense assumed by the c-k method, the resulting k(gA) will not increase monotonically with gA for the other layers. The result of this operation for the CO2 spectrum of Fig. 1 is shown in Fig. 2, where the bottom layer was chosen as layer A. The mapping steps are fast compared to the calculation of k(v). We used the INDEXX routineI to perform the mappings. (iii) We establish subdivision boundaries in gA space for subsequent numerical integration based on the largest range of k, Ak,,j, which will provide results in the radiation calculation with errors no larger than L. The quantity Ak,,j can vary with each layer, and with each subinterval (I) of gA . The magnitude of Ak,,j for a given t will depend on the details of the k spectrum and on the path length distribution and gas abundance and cannot be specified here. The user must determine the appropriate step size. We are not aware of any systematic way to determine Ak. In our tests we derived Ak by trial and error. (iv) We examine the range of k within each g, subdivision for every other layer in the atmosphere. If, for every layer, and within each subdivision of g,, the range of k is no larger than Ak,,j defined in the previous step, then go to step (viii). (v) We choose another layer, layer B, and map g, onto a new variable, g,, , such that within each subdivision defined in step (iii), k,(g,,) increases monotonically. This mapping is performed with a sorting algorithm operating on one g, subdivision at a time. The mapping is denoted by TAB, since the k coefficients of layer A define the map TA for the interval 0
0
0.2
0.4
0.6
0.6
1.0
g
Fig. 3. layers
absorption coefficient spectra after sorting k for and 4 according to step (vi) of Method I (see the text).
CO? 7
0
0.2
0.4
0.6
0.8
1.C
g
Fig. 4. CO2 absorption coefficient spectra after sorting k for layers 7, 4, and I according Method I [cf. step (vii) of Method I].
ROBERT WEST et al
194
mean k coefficients, k,,j, for each subdivision and layer in the final g,,,,,, space. This step also requires insight by the user, and depends on which quantities are of interest (transmissions, fluxes, flux divergence, etc.). One way to formulate an expression for k is to require that the transmission through the layer along a mean path be identical to that given by lbl. If fi, the product of the mean path length and absorber amount, is known for each layer, the specified requirement implies -In k=
[i li
16g, exp( - k,C)/Ag
I,
(2)
where 6g, is the size of the g-interval associated with each ki value, and Ag is the size of the subdivision in gABC,., . Equation (2) can be used to define C, by minimizing errors in the desired quantity. Equation (2) reduces to the frequency- (or g-) weighted mean of k in the weak absorption limit. If it is known that the U are small (e.g., if one is only interested in vertical transmission, and the atmospheric layers have small k and small column abundance), the linear average of k is appropriate. In practice a distribution of paths and absorber amounts must be considered, and the width of the path distribution and the value E determine the size of Ak used in step (iii). For numerical tests, we used Eq. (2) by calculating 6 along a slant path with emission angle equal to 45”. The 90-term case developed in Figs. l-4 produced a sequence of & values shown in Fig. 5. 2.2. Method II Our second method uses a simple binning scheme to group together all 6v, within Av that have similar k values at all levels of the model atmosphere. The binning proceeds as follows: (i) for each layer j, we define a value for the maximum acceptable range Ak,,, of k to be combined into a given bin, 1. This quantity can vary as a function of the absorption coefficient. As for Method I, the user must investigate the appropriate Ak,,j by trial and error. For numerical tests, we choose Ak,,, such that 1k - &I < 0.5 6 ,.,. (ii) W e use the k values in the first 6v interval, k,,,, to initialize the width and mean absorption coefficients for the first bin: 6g, = 6v, /Av; k,,, = k,,,. (iii) We examine the k values for the next frequency interval, 6v,, and compare them to the k/,, in existing bins, If, for any layer, j, and for every bin, I, (k,., - I& > 0.5Ak,,j, we create a new bin, as in step (ii). If, for some bin, I, and for every layer, j, 1k,,j - k,,, 1< 0.5Ak ,,,, we add this interval to bin 1, set 6g, = 6g, + 6v,/Av and redefine the mean absorption coefficient for bin 1,j by adding in the contribution from k,.,. As in step (viii) of Method I, the user must investigate the appropriate averaging technique. We calculated a linear (6g-weighted) average of k for the numerical tests described in the next section. (iv) We repeat step (iii) until all k,, are gathered into bins. This simple scheme is extremely fast because it only involves differences and comparisons. It also requires very little storage because it operates on only one monochromatic spectral interval at a time. The total storage is L(2 + J) + J, where L is the number of bins and Jis the number of model atmosphere layers. This algorithm generated 123 terms for the spectrum shown in Fig. 1, with Ak = k. These are shown in Fig. 6 for three layers. 1
-101
0
0.2
0.4
0.6
0.6
1.0
9
Fig. 5. CO? mean absorption coefficients for layers 1,4, and 7, calculated from Eq. (2) by using the mapping produced by Method I [cf. step (viii) of Method I].
I
-101
0
0.2
0.4
0.6
0.6
1.0
9
Fig. 6. CO, mean absorption coefficients for layers 1,4, and 7 from the binning produced by Method II.
I
Mapping transformations
for RT calculations
3. NUMERICAL
195
TESTS
We performed a series of calculations to test the accuracy of our methods. To compare our methods with the c-k technique, we calculated radiation fields for atmospheric models used by GWCC. These calculations use COZ and HZ0 opacities in the 5000- to 5050-cm-’ region (Rothman et al15). The minimum step size in frequency is 0.25 times the Doppler width at the lowest temperature (layer 1). In the calculation of thermal emission we follow GWCC and set the Planck function to its mean value in the 50-cm-’ interval to isolate errors intrinsic to the sorting algorithms. We begin with calculations for a pure CO* atmosphere (Case A of GWCC), and follow this with calculations for mixtures of CO, and H,O (beginning with Case 0 of GWCC). Temperature, pressure, and gas abundances for each of the models are given in Table 1. Layers 2,4, and 6 contain isotropically scattering aerosols. Aerosol extinction optical depth is 10 in each layer. Aerosol single scattering albedos are 1.0 except for Case A, where they are 0.9999, 0.999, and 0.995 for layers 2, 4, and 6, respectively. Solar and thermal radiation fields were calculated with doubling and adding codes ” having 22 streams. We examined the same quantities (up and down solar and thermal fluxes, net fluxes, and flux divergences) that were discussed by GWCC for the c-k method. All methods, including c-k, produced results which agreed with the lb1 calculations to within about 1% for most quantities except the thermal flux divergence for layers 1 and 2, where the error in the c-k method was as large as 30%. We show in this paper only thermal flux divergence errors. 3.1. Single gas Figure 7 shows thermal flux divergences for Case A. Errors arising from the various methods are listed in Table 3. The cases in Tables 2 and 3 differ in the number of intervals for which E is computed. The c-k calculation examined by GWCC used 500 intervals, but accuracy does not improve for more than about 20 intervals. Two sources contribute to the errors. First, if the Ak are too large in any layer, the range of k will be too large to accommodate a wide range of path lengths. That effect imposes a minimum number of sub-intervals for each layer, of order 10 for Case A if the desired error is a few percent or less. It is the dominant source of error in Method I for layer 7 (Table 3). A second source of error is related to the range of k values for neighboring g values in the plot of k vs g in Fig. 2. This source is responsible for the large errors in layers 1 and 2 in the c-k method. It requires that several k coefficients for the top layer be chosen for a given k value in the bottom layer, and vice versa. It introduces a multiplicative factor in the number of terms each time a new mapping must be carried out according to step (v) of Method I. The same effect influences the number of terms in Method II. Case A required sorting on the bottom, middle, and top layers to achieve errors of less than a few percent using Method I. The case with the minimum number of intervals has 10 divisions for the bottom layer, times 3 for the middle, times 3 for the top, or 90 total. 1”~‘1”“J~!,‘1’1~~
Table 1. Model atmosphere layer parameters.
+
2x10-1' iti 5g\
0 ______~_____i___‘__~____._______.__________._____._.___
gg 4 E 2xd'
The
parameterfis the mass mixing ratio of CO,. The mixing ratio of H,O is 1-J The surface is a Lambert reflector with reflectivity and emissivity 0.5.
+ + + I,.,.IIIII1ll~.ll.IIIIIIII~IIIIIII~ 012345678
Layer Fig. 7. Thermal flux divergences for Case A.
ROBERTWESTet al
196 Table 2. Single gas test cases.
ays
Whod
-z-
c-k
90
3
135
3
150
5
160
3
225
3
250
5
300
3
123
-
?-
Akii
1
LaVEr
hlethoi
-K
c-k
2
3
T-
T-
-.-L-
4 __ < 0.10
-ii-
6
0.10
-K2
7 0.17
10
90
-0.21 -0.76 -0.51
0.39
0.63
15
135
-0.55 -0.76 -0.49 < 0.10
0.26
0.59
1.6
10
150
-0.46 -0.26 -0.20 < 0.10
0.37
0.61
2.4
20
180
-0.69 -0.69 -0.47 < 0.10
0.16
0.36
1.2
15
225
0.63
-0.39 -0.24 < 0.10
0.26
0.56
1.6
10
250
0.49
-0.17 -0.14
0.32
0.56
2.2
20
300
0.49
-0.41 -0.36 < 0.10
0.14
0.34
1.1
1.72
-0.13 -0.69 -0.21
-1.0
-0.49
-0.63
1.0
-
Table 3. Thermal flux divergence errors, single gas.
The total number of terms in the k distribution are listed in the second column, the number of subdivisions for each of the top, middle, and bottom layers for Method I are listed in the next three columns. The bin size criterion, Ak/k, for Method II is listed in the final column.
123
-
2.0
-
The error is 100 x (sorting algorithm - lbl)/lbl. The c-k results were obtained by GWCC.
The several cases for Method I listed in Tables 2 and 3 were designed to examine the errors as the number of terms increases. The error for layer 7 decreases as the number of divisions for that layer increases from 10 to 20. Errors in the other layers are all < 1% and do not show any systematic relationship to the number of subdivisions. It is difficult to know the sign and magnitude of the various contributions to error once the error becomes 5 1%. It is also difficult to justify the computational expense of achieving higher accuracy because errors or uncertainties in gas absorption coefficients and properties of the atmosphere impose greater limitations. Both Methods I and II produce comparable errors. 3.2. Gas mixtures We tested two approaches to the problem of gas mixtures. The first is to calculate the k spectrum for the gas mixture and treat it as a single-gas spectrum using the methods described in Sets. 2 and 3. We used gas abundances for Case 0 (cf. Table 1) to test this method. The CO* mixing ratio varies from 0.2 in layer 1 to 0.95 in layer 7. Calculations using our sorting algorithms produced results with errors of a few percent or less, listed in Table 4. This method requires that a new mapping be produced each time the mixing ratios change. We devised a second, more flexible, method for gas mixtures. It is a perturbation approach which begins by calculating the k spectrum for a model with gas mixing ratios having mean values in each layer. The sorting algorithm is then executed to map the combined-gas spectrum according to the procedures described in Sets. 2 or 3. Instead of forming L for the combined gas, we calculate lists of E for each gas separately and store them. When models with different mixing ratios are constructed, 6 values for each gas are combined as in the lb1 procedure. We tested this method with four models. The first, Case T is a pure CO, atmosphere (same as Case A). The second, Case U, is a pure Hz0 atmosphere. In the third, Case V, the mixing ratio of the minor constituent in each layer is increased by 30% over that for Case 0. The fourth, Case N, (this case was studied by GWCC) contains equal amounts of both gases in all layers. The lists of E for the separate gases were calculated using the mixing ratio profile of Case 0. The errors for each of these cases are listed in Table 5. Errors remain small when the mixing ratio profile differs by 30% from the profile used to compute E, as the results for Case V demonstrate. They can become quite large for large departures (Case T) because the sorting was performed on a combined-gas spectrum and is not optimum for a single gas. Errors for Case N, in the intermediate regime, are a little larger than those for Case V. The errors produced by the perturbation approach described above can become significant when the mixing ratio profile departs by more than a few tens of percent from the value used to calculate the E coefficients. In those cases one can use the first approach described in this section, but a more flexible method is desired.
Mapping transformations Table 4. Thermal flux divergence errors, Case 0. I
1
a”ar
for RT calculations
Table 5. Thermal flux divergence errors, Cases T, U, V, and N. C=e
Tems
layer 1
T
1
j
-0.17 -0.61 -0.66
/
-1.1
j
1.3
197
300
-74
1 2 1 3 1 4 1 5 1 6 1 7 ) -29 ) -14 1 -5.0 1-0.16 /co.101 4.2
j
The c-k results were obtained by GWCC.
If the spectra of the two gases are uncorrelated it is possible to generate a set of k coefficients for gas mixtures from the lists of Sg, and E,,jfor each gas. The individual Eare combined by creating a new set of intervals (denoted by g,,*,,). For each g interval in the listing for gas 1, 6g,.,, there are N 6g,,2,, intervals for the gas mixture corresponding to N 6g,,l intervals for gas 2. If the number of 6g,,, intervals is M, the number of 6g,,2,, intervals is M x N. The &j for the combined gas for each of these intervals is the average of the &,,j (for gas 1) and &,j (for gas 2) weighted by the mixing ratios of the gases in layer j, where 1= m x (n - 1) + n, and interval m for gas 1 was combined with interval n for gas 2 to produce interval 1 for the combined gas. The it4 x N intervals are then mapped and regrouped by Methods I or II to produce a much smaller number ( N the larger of M or N) of intervals for use in the radiation calculations. This method requires little storage capacity because only the E for each gas are stored. It is flexible because it is not restricted by the mixing ratio profile. Its accuracy is limited by the degree of correlation of k coefficients for the two gases. We have not tested this method for the mapping algorithms derived here, and we cannot predict how well it will perform in general. 4.
TEMPERATURE
VARIATIONS
It is cumbersome to recalculate the lb1 k spectrum for each new temperature-pressure profile. A much more efficient approach is to calculate k coefficients on a temperature grid and to interpolate each time a new set of k coefficients is needed. Goody et al” demonstrated the method for the c-k technique. Since the number of k values which are stored and used in the interpolation is small, the method is very efficient and does not require large storage capacity. The same technique can be used with the algorithms presented here. The same mapping must be used for each temperature to ensure that the spectral content of each of the bins does not change with temperature. We tested the interpolation method with a CO, hot band (where the temperature sensitivity is strong) in the region 595-600 cm-‘. We compared thermal flux divergences for a temperature profile appropriate for the Martian atmosphere (T varies from 140 K at altitudes above the 5-p bar level to 2 14 K at the 6.4-mbar level). The interpolation used k coefficients calculated at temperature intervals of 30 K. The largest difference in thermal flux divergence between the calculation using interpolated coefficients derived from Method I, and the non-interpolated lb1 calculation was 3%. We consider that difference to be acceptable since the lb1 absorption coefficients are not known to that accuracy. 5. DISCUSSION
We have introduced two spectral mapping algorithms for broadband radiation calculations in planetary atmospheres. It is appropriate to examine their relationship to each other and to another mapping algorithm, the c-k method. The algorithms described here are quite similar. Their goal (shared with the c-k method) is to reduce the number of radiation calculations by grouping together spectral regions with similar values of k. The two methods differ in the order in which the steps are executed. Method I performs a series of mapping transformations based on the k values of the various atmospheric layers, and then combines the k values at the end to form mean k values. Method II begins by defining criteria
198
ROBERT WFST et al
(the A/c,.~)which determine which group a given spectal interval will be mapped onto. Method I retains the full spectral information (list of v and k values) until the final step to compute 5 whereas Method II discards that information (unless it is needed for other mixing ratios or temperature profiles) after each spectral interval has been assigned a bin. The mapping transformations performed by Methods I and II are not necessarily equivalent. By equivalent we mean that the k,,j are formed from the same set of ki,j. If the grouping criteria were defined in the same way for both methods, the resulting radiation calculations would be identical for both. The two methods as described above are not equivalent because the criterion for bin membership in Method II depends on a moving average of k. A spectral interval which satisfied the criterion Ik,,, - &,,,I6 OSAk,,, when it was tested may not satisfy that criterion for the final value of E,,,. The moving average is not an essential feature of Method II. We also use a variant of Method II which has fixed Ak intervals and is equivalent to Method I. The mapping algorithms we introduced here differ fundamentally from the c-k technique whose validity is restricted to the weak absorption limit, or to cases where correlation assumption is valid (i.e., Elsasser models or spectra which obey a scaling law”). When the correlation assumption is strictly valid, all three mapping transformations are equivalent (except for the nonessential difference between Method I and one form of Method II noted above). The mapping algorithms described in the Sets. 2 and 3 can therefore be considered to be more general ones, whose validity is unrestricted apart from the requirement that the source radiation be nearly constant over the interval. Of course, the number of intervals needed to achieve arbitrary accuracy may be large. In the limit Ak -+ 0, the lb1 results are reproduced and the number of g intervals can be as large as the number of v intervals in the original spectrum. More work needs to be done to understand how to make optimum choices for the size of k intervals needed to achieve a desired accuracy. This problem was encountered earlier for the c-k method. It is more complex for these new methods since there are more degrees of freedom. Only the number and distribution of g-intervals was in question in the c-k method. In the new methods one must explore the trade-off between the Ak for one layer vs the Ak for another. We currently investigate these issues by trial and error. We have described two techniques to reduce the number of terms to a minimum number (cf. Sec. 2). Although these techniques are useful, a more systematic approach (e.g., cluster analysis) is desired.
6.
SUMMARY
We have described two methods of mapping the lb1 k spectrum to produce a set of k coefficients which can be used for accurate broadband radiative transfer calculations in non-gray, vertically inhomogeneous scattering atmospheres. These methods are intended for applications where the less accurate c-k method is inadequate. The algorithms described here allow the user to achieve results arbitrarily close to the lb1 results. They reduce the number of radiation calculations by about two orders of magnitude, compared to the lb1 technique, for the cases we tested. These algorithms can be used for mixtures of gases, and for a variety of pressure-temperature profiles. The two methods we introduced are similar, but differ primarily in the way in which they are executed. Method II requires little memory since it operates on only one entry in the lb1 spectrum at a time. Knowledge of the wavelength mapping is not retained, although it could be if desired. Knowledge of the mapping of wavelength onto g is preserved in Method I until the final step, when averages of k are formed. This requires significant memory during the mapping, although the mapping is fast. Techniques for combining gas absorption coefficients and interpolating on temperature grids are identical for both methods. These methods are essentially equivalent to the c-k technique for Elsasser band models and for cases where the absorption coefficient obeys a scaling law (e.g., the strong-line limit). The mapping transformation developed for Method I illustrates why the c-k method fails (Fig. 2). Our numerical tests underscore the need to obtain high-resolution spectroscopic measurements in the laboratory. Much of the high-frequency variation seen in the upper curve in Fig. 2 is produced in the near-wings of lines. Calculations which ignore this effect may produce errors of order several tens of percent.
Mapping transformations
for RT calculations
199
have benefited from a close collaboration with Professor Richard Goody, especially in the early stages of the work described here. The support given by Daniel McCleese was an essential ingredient in the success of this oroiect. This work was nerformed bv the Jet Prouulsion Laboratory, California Institute of Technology, under contract bith the National Aeronautics and Space Administration. The authors were supported under contracTs administered by the NASA Planetary Atmospheres program.
Acknowledgemenrs-We
REFERENCES 1. W. C. Yang, Y. L. Yung, A. A. Lacis, T. MO, and J. E. Hansen, Science 194, 685 (1976).
2. J. E. Hansen, W. C. Wang, and A. A. Lacis, Science 199, 1065 (1978). 3. A. Lacis, J. Hansen, P. Lee, T. Mitchell, and S. Lebedeff, Geophys. Res. Lett. 8, 1935 (1981). 4. J. Hansen, G. Russell, D. Rind, P. Stone, A. Lacis, S. Lebedeff, R. Ruedy, and L. Travis, Mon. Weath. Rev. 111, 609 (1983). 5. J. Hansen, 1. Fung, A. Lacis, D. Rind, S. Lebedeff, R. Ruedy, and G. Russell, J. Geophys. Res. 93,934l (1988). 6. D. Rind, R. Suozzo, N. K. Balachandran, A. Lacis, and G. Russell, J. Atrnos. Sci. 45, 329 (1988). 7. V. Ambartzumian, Publs Obs. Astr. Univ. Leningrad 6, 7 (1936). 8. A. A. Arking and K. Grossman, J. Atmos. Sci. 29, 937 (1972). 9. G. A. Domoto, JQSRT 14, 935 (1974). 10. A. A. Lacis and V. Oinas, submitted to J. Geophys. Res. (1990). Il. R. M. Goody and Y. L. Yung, Atmospheric Radiation: Theoretical Basis, Oxford University Press, Oxford (1989). 12. R. Goody, R. West, L. Chen, and D. Crisp, JQSRT 42, 539 (1989). 13. W.-C. Wang and G.-Y. Shi, JQSRT 39, 387 (1988). 14. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes, Cambridge University Press, Cambridge (1986). 15. L. S. Rothman, R. R. Gamache, A. Goldman, L. R. Borwn, R. A. Toth, H. M. Pickett, R. L. Poynter, J. M. Flaud, C. Camy-Peyret, A. Barbe, N. Husson, C. P. Rinsland, and M. A. H. Smith, Appl. Opt. 26, 4058 (1986).