Introduction

The unexpectedly low microwave brightness temperatures of the Greenland and Antarctic ice sheets, as measured by the electrically scanning microwave radiometer (ESMR) on Nimbus-5, have been the subject of several theoretical investigations in recent years. Reference Chang, Gloersen, Schmugge, Wilheit and ZwallyChang and others (1976) showed that, for a uniform snow medium, volume scattering is a dominant factor affecting the microwave emission and that the source of radiation can emanate from several meters within the medium. Subsequently, taking into account variations of snow properties with depth, Reference ZwallyZwally (1977) used an analytic approximation to solve the microwave emissivity of the polar firn. The approximation overestimated the effect of radiative scattering in the medium because contributions due to multiple scattering were neglected. The scattering coefficient ys had to be multiplied by an empirical parameter f = 0.12 to obtain agreement between calculated emissivities, e(model) and observed emissivities e(obs). Using a different numerical technique for solving the radiative transfer equation from the one used in this paper, Chang and others (1980) showed good quantitative agreement with their quoted brightness temperatures from Nimbus-5, but the quoted values differ substantially from other published values Reference Zwally and Gloersen(Zwally 1977, Zwally and Gloersen 1977).

This paper presents the result of a numerical solution of the radiative transfer equation similar to the method of Reference Chang, Gloersen, Schmugge, Wilheit and ZwallyChang and others (1976), but with depth-dependent parameters. When a value of f = 0.3 is used to adjust the scattering coefficient, the results agree with observed microwave emissivities and qualitatively with Zwally's results. In addition, the model is used to simulate the seasonal variation of the microwave emission. Comparison of modeled and observed seasonal amplitudes provides additional information on the validity of the radiative transfer model and, in particular, of the values of the scattering and absorption parameters used in the model.

Radiative transfer model

The intensity (energy flux per unit wavelength per unit solid angle) of radiation 1(θ) passing through a medium at an angle 6 with the vertical can be calculated from the equation of radiative transfer given by:

where T is the optical depth of the medium, ε is the emittance, u0 is the single scattering albedo, B is the Rayleigh-Jeans approximation to the Planck function, μ = cos θ, P(θ) is the scattering phase function, and T' and 8' are the integration variables for optical depth and scattering angle respectively Reference Swihart(cf. Swihart 1971). The first term on the right-hand side of the equation is the contribution of radiation incident on the medium from the bottom and attenuated by the medium before it reaches the surface. The second term, the non-scattering term, takes into account emission and absorption within the medium, whereas the third term is the contribution to the radiation due to scattering. The brightness temperature T[3 is related to intensity by:

where c is the speed of light, k is the Boltzman constant, and A is the wavelength of the radiation.

The integral equation was solved numerically by dividing the medium into 1000 layers of unequal optical depths and assuming that the depth-dependent properties do not change within any one layer. The layering structure is a numerical technique only and is not intended to relate to any physical properties of the medium. The layer thickness increases from 0.001 optical depth per layer for the top 200 layers to 0.002, 0.005, 0.01, and 0.02, respectively, for increments of 200 layers. This procedure permits more accurate computation near the surface, where a large amount of the outgoing radiation originates. The calculations were not extended to the bottom of the ice sheet, but rather to 7.5 optical depths (corresponding to an average depth of 25 m) because contributions to the outgoing radiation from further down were negligible.

Angular dependence of the intensity is controlled by the phase function, which, in this model, is the Rayleigh phase function given by:

where θ is the angle between the scattered radiation and the direction of incidence. The angular variable was divided into 72 equal divisions, each of which is associated with an intensity that, in turn, is updated from layer to layer to provide the angular distribution of the radiant energy.

The numerical procedure is similar to the Gauss- Seidel iterative technique employed by Reference Herman and BrowningHerman and Browning (1965) in atmospheric studies. Assuming isotropic up-welling radiation at the bottom of the medium, the calculation proceeded upward layer by layer. At each layer, the contributions of the upwelling components for each angle were summed, and the down-welling components were stored for a subsequent iteration. At the top layer, it was assumed that no downward radiation is received from outside the medium. The down-welling radiation was then calculated layer by layer, similar to the previous calculation but with summation of down-welling components and storage of up-welling values. This procedure was repeated until the percentage increment in net radiation at the surface was less than 0.05%. The average number of iterations that was made before convergence was about 10. Various tests were made to check the accuracy of the numerical solution. One such test case is the calculation of emission as a function of angle from a pure scattering medium with Rayleigh phase function. The results from this model agreed very well with the analytic calculation of Chandrasekhar (1960) as shown inFigure 1. Although the model calculates brightness temperature as a function of emission angle, the Tg discussed below are nadir values to match the electrically scanning microwave radiometer (ESMR) observations.

Fig.1. Comparison of model calculations with the analytic results of Reference ChandrasekharChandrasekar (1960) for a pure scattering medium with a Rayleigh phase function.

Depth-dependent parameters

To evaluate the radiative transfer equation, the emittance, the single scattering albedo, and the physical temperature of the medium in each layer must be known. The emittance and the albedo at _a particular optical depth T are functions of the absorption coefficient y_a and the scattering coefficient ys of the particle as follows: where y_e = y _a + y _s is the extinction coefficient. The optical depth is given by T (Z) = s^z re(_z)dz, where z is the depth in meters. For ice particles at _a particular temperature, Mie scattering calculations show that y _a depends greatly on the imaginary part of the refractive index, ", and is insensitive to crystal size for radii <1 mm Reference Zwally(Zwally 1977). Also, although Reference CummingCumming's (1952) measurements of the loss tangent at 3.2 cm indicate _a strong dependence of the loss tangent on temperature from about 258 to about 273 K, the loss tangent approaches an almost constant value at temperatures less than 258 K. Because the temperatures relevant to this study are mostly less than 258 K, the same value of y a was therefore used at all seven stations and was kept constant as a function of depth in most of the calculations presented.

Nevertheless, both the magnitude of " and its temperature dependence at microwave wavelengths are not well known. A compilation by Reference EvansEvans (1965)shows some variations of the loss tangent with temperature down to about 213 K. However, the measurements from different investigators in the compilation were not consistent with each other. The Debye equation Reference Hobb(Hobbs1974:81-199) also shows a temperature dependence of the dielectric constant of ice, but, near the ESMR wavelength, It does not agree with Evans.Figure 2compares interpolated y_a at 1.55 cm wavelength based on Evans1 compilation, those derived from the Debye equation using an activation energy of 0.45 eV, and the constant value used in our model calculations. The absorption coefficient is considerably more sensitive to temperature in the Debye equation than in Evans' compilation. The Debye equation could also give a wide range of values for a particular temperature, depending on the activation energy, for which measured values for firn and ice range from 0.45 to 0.62 eV Reference Hobb(Hobbs 1974: 81-99). The effect of a possible temperature-dependent Y_a on model results is discussed below.

Fig.2. Temperature dependence of the coefficient of absorption as derived from the Debye equation and Evans (1965) at 1.55 cm. A constant value was used in the model.

The scattering coefficient Ys was determined from the size of the particles and the dielectric constant of the medium, assuming Rayleigh scattering (Van de Hulst 1957: 63-84). In this model, following Reference Zwally and GloersenZwally (1977), we used the depth-dependent expression of Y_s=(l-8 r)³ with r³ = r₀ + az based on Rayleigh scattering cross-sections and on the crystal size data of Reference GowGow (1969) for different snow depths and for various locations in Greenland and Antarctica. To simulate the seasonal dependence of brightness temperature, a time-dependent temperature profile appropriate for the firn at Maudheim Reference Dalrymple, Lettau, Wollaston and Rubin(Dalrymple and others 1966) was used:

where T_m is the mean annual surface temperature, t is the time, and (2A) is the amplitude at the surface. Although the temperature profile varies from one station to another, Equation (6) was used with appropriate values of T _m and A for each station.

Discussion of results

The model was tested for seven different locations in Greenland and Antarctica where depthdependent data are available. To compare results with the data, the mean bulk emissivity e =T_B/T_m Reference Zwally and Gloersen(Zwally 1977) was determined for both observed and calculated brightness temperatures. The Tg is a yearly average. The observed emissivity e(obs) was computed by using measured Tm and an average of Nimbus-5 ESMR brightness temperature data from June 1973 to May 1975. Table I lists the observed e(obs) values, which differ slightly from those reported by Reference ZwallyZwally (1977) because they were determined from an improved and larger data set of observed brightness temperatures.

For e(model), Tm was also used and T was obtained from a solution of Equation (1). In a test case, Tg was calculated for 12 months, using Equation (6) for the temperature profile. The average of the 12 monthly values was the same as that of the Tg obtained by using a constant Tm in Equation (6) rather than a time-dependent profile. Therefore, to calculate the mean emissivity, it was not necessary to use a temperature profile, except for the case in which Ya was made temperature-dependent. Density is another depth dependent parameter that was set at a constant value 480 kg m-3, as was done by Reference Zwally and GloersenZwally (1977). In a test case, the firn density was varied ±20%, and the calculated variation in brightness temperature, and therefore emissivity, was <2%. This insensitivity to density occurs because both Ya and Ys are proportional to density so that their ratio Is independent of density. Therefore, the variations in density from station to station and the variations in density with depth were neglected in the calculations of Tg

The single most important variable in this model that affects the mean emissivities Is the grain-size distribution that was used to calculate Ys at eacn location. Of the firn parameters considered, grain size is the parameter that varies most with location. The initial calculations using Ya = 0.15 m"1 (based on Cumming's (1952) value of ") and Ys = (l«8r)3 gave calculated emissivities e(model) that differed significantly from e(obs) for most stations. Because there is substantial uncertainty in both Ya ar,d Ys» these coefficients were adjusted as follows. First, a value was chosen for the absorption coefficient Ya and an empirical factor f was used to vary the scattering coefficient (Y$ = f Ys) until the differences between E(model) and e(obs) at four selected stations (South Pole, Plateau, South Ice, and Site 2) were minimized. The YS was then kept constant, and Ya was varied until a minimum rms deviation was obtained.Figure 3 shows the effect of varying Ya and f at all stations on the calculated emissivities. Varying Ya W1tn constant ys changes the slope and the magnitude of the E(model) versus e(obs) relationship, whereas varying YC with constant Ya primarily changes the magnitude of the e(model). At larger values of Ya» the required value of f increases but, contrary to the observations, the e(model) tends to be the same for all locations. This tendency occurs at larger values of Ya because the calculated radiation comes from nearer the surface, causing the sensitivity to variations In grain sizes at lower depths to decrease.

Fig.3. Observed versus modeled emissivities for various values of (a) absorption coefficients Ya and (b) scattering coefficients Ya = f TV

The values chosen for f and Ya a"*e 0.30 and

Table I Comparison of observed and calculated emissivities

0.038m⁻ⁱ, andFigure 4 shows the resulting correlation. The correlation coefficient is 0.986, and the standard deviation is 4.6%. This value of f is roughly three times greater than that of Zwally's analytic solution, which required a small f because of its overestimate of the effect of scattering. However, an f less than unity means that the amount of scattering is still over-estimated in the model, despite the fact that multiple scattering was considered. Therefore, the calculated y_s apparently does not reflect the true scattering cross-sections of the snow grains in the firn. This suggests a limitation in the Rayleigh and Mie scattering models, which use a plane-wave approximation and assume independent spherical-scattering centers. Although the absorption coefficient is a factor of four lower than the estimate made from the 3.2 cm data of Reference CummingCumming (1952), it agrees well with the interpolated value from Evans (1965) at 1.55 cm wavelength and 232 K.

Fig.4. Comparison of observed and calculated mean emissivities at seven locations.

As a further check on the reliability of the model and the choice of Y_a and Y_s tne seasonal distribution of the brightness temperature was modeled. The results for Plateau and South Ice stations are shown in Figure 5. Plotted with the model results are the ESMR seasonal brightness temperature data for 1973, Reference Hobb1974, 1975, and 1976. These seasonal variations are useful because a large difference occurs in the depth-dependent temperature profile from summer to winter.

Fig.5. Observed and modeled seasonal variation of brightness temperatures at Plateau and South Ice stations. Model (dotted line) is compared with 4 a of ESMR TB data.

Because the amount of volume scattering controls the transmission of radiation from within the medium, the seasonal variation is an independent check on the Y_a and y_s used in the model. For example, higher Y_a gives higher seasonal amplitude because the radiation in that case comes from nearer the surface, where substantial seasonal variability occurs in the firn temperature. Except for September, October, and November, the agreement between theory and observed seasonal distribution for South Ice is good. The agreement is also good for Plateau except January, February, and December for which the data showed significantly higher values than the model. A better agreement could probably be obtained for Plateau by using a time-dependent temperature function approp-Mate to the East Antarctic plateau. On the other hand, increasing the y_a to 0.15 m⁻¹ for example, or using the Debye formula substantially increased the amplitude of the seasonal cycle above observed values.

A better correlation for the emissivities in Figure 4 could be obtained by using a smaller absorption coefficient Y_a = 0.02 m⁻¹ with y's = 0.12 Y_S. However, this set of values would have a small seasonal variation in brightness temperature, and the correlation with observed data would not be as good. Furthermore, Ya= 0.02 m⁻¹ is consistently lower than both measurements from Reference CummingCumming (1952) and Evans (1965), and Ys = 0.12 Y_S is even further from calculated values.

A set of calculations was made using a temperature-dependent ya that varied with station temperature. The temperature dependence of Y_a in the Debye equation is quite strong, and using it in the model gave results that were inconsistent with observed emissivities and seasonal amplitudes of brightness temperature. The interpolation of dielectric constant from Evans (1965), shown in Figure 2, was also used to obtain a value of Y_a based on the mean annual temperature for each of the seven stations. These values for Y_a were also used in the model, and the emissivities obtained are shown in Table I. The results show that a good correlation at all seven stations is not possible with this temperaturedependent Y_a. This fact suggests that the temperature dependence of Y_a is negligible as was assumed.

Conclusion

These results show that both the abnormally low brightness temperatures and the lack of simple correlation of these temperatures with physical temperatures can be successfully modeled by using a numerical solution of the radiative transfer equation with Rayleigh scattering. A consistently good correlation between the calculated and observed bulk emissivities at seven different stations can be obtained if the absorption coefficient, Y_a iS assumed to be identical at all these stations. This indicates that in the temperature range applicable for the ice sheet the temperature dependence of Y_a is either very small or negligible, which is consistent with the loss tangent measurements of Reference CummingCumming (1952). The results also show that modeling of the seasonal dependence of brightness temperature can be a useful aid in selecting values of the scattering and absorption coefficients. In particular, the seasonal variation of brightness temperature is strongly dependent on emission depth, which, in turn, is controlled by the values of the scattering and absorption coefficients. Although the value chosen for Y_a is quite consistent with some measurements of the index of refraction, the value for y_s obtained is substantially smaller (~70%) than Rayleigh scattering calculations based on measured grain size. This is not unexpected because the use of Rayleigh and Mie scattering models have some known limitations when applied to snow grains as scatterers. In particular, because the Rayleigh/Mie models use plane-wave approximations and assume independent spherical-scattering centers, exact quantitative agreement should not be expected.

Acknowledgement

This research was supported by the National Aeronautics and Space Administration (NASA) Climate Program.