Background

The discrete ordinates method of Reference ChandrasekharChandrasekhar (1950) gives solutions of any desired degree of accuracy to the equation of radiative transfer for a plane parallel, homogeneous, multiple scattering medium. The discrete ordinates, or streams, refer to the angles at which radiances are determined. Including more streams gives greater accuracy but also gives a greater degree of computational complexity. From Reference ChandrasekharChandrasekhar (1950) the equation of radiative transfer for a plane parallel medium is

(1)

where I is the radiance, S is the souce function, φ is the azimuth angle, and μ is the cosine of the zenith angle θ. The optical depth τ is defined in terms of the volume scattering coefficient σ and the absorption coefficient k by

(2)

The minus sign results from the choice that positive μs represent upward angles. For the case of sea ice at visible wavelengths, thermal emission is negligible and S describes the contribution of multiple scattering. For a plane parallel scattering atmosphere with a direct incident beam

Axial symmetry holds for the experimental results under study and Equation (1) gives

where E ₀ is the direct solar irradiance and μ ₀ is the absolute magnitude of the cosine of the solar zenith angle.

The scattering processes within the medium are jointly defined by the phase function p(μ, μ’) and the single scattering albedo ϖ ₀. The phase function describes the angular dependence of scattering and typically is normalized so that its integral over solid angle is equal to one. For use in the discrete ordinates method, it is expressed in terms of the following expansion in Legendre polynomials,

(3)

where N is the number of streams. The relative values of the ϖ₁ s in Equation (3) are obtained for a particular phase function from

(4)

which is derived by exploiting the orthogonality of the Legendre polynomials. These values must then be multiplied by the single-scattering albedo which describes the relative importance of scattering and absorbing processes and is defined as

(5)

and ranges from zero for a purely absorbing medium to one for a purely scattering medium.

Reference LiouLiou (1973, Reference Liou1974) adapted Chandrasekhar’s general solution to investigate radiative transfer within cloudy and hazy atmospheres. Analytic solutions were determined for the two- and four-stream cases with numerical solutions for eight and sixteen streams. In two-stream models the errors became prohibitively large for optically thin layers. The eight- and sixteen-stream models, while offering greater angular resolution of radiance and a more precise representation of complex phase functions, required a complicated iterative solution scheme. Thus, the four-stream case was chosen to take advantage of the mathematical convenience while retaining sufficient accuracy (Reference LiouLiou, 1974). For this model, Reference LiouLiou (1974) gives the solutions for the upwelling and downwelling irradiances in terms of a four-point Gaussian quadrature as follows:

(6)

and

(7)

where I(τ, ±μ) is the radiance in each of the streams, the a_j s are the Gaussian weighting functions (a ₁ = 0.6521452 and a ₂ = 0.3478548), and the μs are the cosines of the zenith angles for each of the four streams (μ ₁ = 0.3399810 and μ ₂ = 0.8611363).

Application to sea ice

In general the ice structure is assumed to be homogeneous and is represented by a single-layer model. It is also necessary to take into account that the medium is a collection of scattering centers embedded in ice rather than in air. Thus, refraction and specular reflection modify the incident radiation field as it enters the ice. To illustrate the importance of refraction and reflection, solutions are carried out both ignoring and including these effects for simple cases.

No refraction

In this case the boundary conditions are

and

where τ _M is the optical thickness of the ice slab. For isotropic incident radiation representing cloudy sky conditions, C is equal to one and there is no direct beam component (E ₀ = 0). To represent clear skies, C is assumed to be zero and the incident radiation field is defined by a direct beam of radiance E₀ at an angle arccos (μ ₀). The lower boundary condition is set by assuming that the underlying water does not scatter radiation back into the ice.

The albedo and transmittance are determined from α = F↑(0)/F ₀ and T = F↓(τ _M)/F ₀. F ₀ is the incident irradiance and is equal to π for isotropic incident radiation and μ ₀ πE ₀ for an incident direct beam. Equations (6) and (7) are used to calculate F↑(0) and F↓(τ_M).

Refraction

At the ice-air boundary, refraction must be considered in detail. From geometric considerations (Reference JerlovJerlov, 1976) or variable transformation theory (Reference PapoulisPapoulis, 1965), the redistribution of a diffuse incident radiation field by refraction is given by

where R(−μ ^a) is the Fresnel reflection coefficient for natural light and μ ^a and μ are corresponding direction cosines in the air and the ice respectively. For isotropic incident radiation I _air(0, −μ ^a) = 1 and μ ^a ranges from −1 to 0, but μ extends from only −1 to −0.646. Figure 2 displays the redistribution of this radiation field at the ice-air interface. Before refraction, the radiance is one at all angles. Specular reflection removes much of the radiation near grazing incidence, and refraction concentrates the remaining downwelling radiance between 0 and 49.8°.

Fig. 1.

Schematic representation of the four-stream model for a floating ice slab including refraction and specular reflection at the upper boundary.

Fig. 2.

The angular distribution of isotropic incident radiation undergoing refraction: (1) isotropic incident radiation; (2) net incident radiation; and (3) radiation distribution in the ice after refraction (θ = 49.5° is the angle of total internal reflection).

After refraction, I _ice(0, −μ ₂) = 1.67 while I _ice(0, −μ ₁) = 0. As a result the quadrature in Equation (7) has only one non-zero term remaining which introduces considerable inaccuracy in the calculation of downwelling irradiance. In order to conserve energy, the value of I_ice(0, −μ ₂) is adjusted to make Equation (5) consistent with the formula (Reference ChandrasekharChandrasekhar, 1950).

Numerical integration gives F↓_ice(0) = 0.937π so that I _ice(0, −μ ₂) = 1.56 which is used in Equation (7) in place of 1.67.

The upwelling radiance at τ = 0 is similarly influenced by refraction. None of the radiation striking the boundary at angles exceeding the critical angle for total internal reflection (θ _c = 49.8°) emerges from the ice, including I_ice(0, μ ₁). In addition, 2.2% of I _ice(0, μ ₂) is reflected downward by the interface. For diffuse incident radiation then, the upper boundary conditions are

(8)

and

(9)

If the incident radiation field is a direct beam, E ₀(μ _in), the magnitude of the refracted beam is given by

where μ _in is the cosine of the zenith angle of the direct incident beam in air and μ ₀ is cosine of the corresponding angle after refraction (Reference Born and WolfBorn and Wolf, 1964). Since this beam is included explicitly in the equation of transfer as a volume source term, the boundary conditions are the same as for the diffuse case except that the 1.56 does not appear in Equation (9).

At the bottom interface, the small difference in the refractive indices of ice and water gives a very small amount of specular reflection. For example, only 0.26% and 0.01% of the radiation reaching the lower boundary at angles −μ ₁ and −μ ₂ is reflected back into the ice. Consequently, the effect can be ignored, and the lower boundary conditions remain the same as for the no-refraction case.

The formula for total transmittance retains the same form, T = F↓(τ _M)/F ₀, but calculating the albedo is more involved. The general expression for albedo is

(10)

For the direct-beam case, the first integral in the numerator is equal to πμ _in R(μ _in)E ₀(μ _in). For diffuse incidence, numerical integration gives a value of 0.065π. The second integral would normally be evaluated by Gaussian quadrature, but because I(0, μ ₁) is totally reflected, only one non-zero term would be present with a resulting loss in accuracy. Therefore the upwelling radiance is fitted to a function of the form I _ice(μ)=Aμ+B (as is the Eddington assumption (Reference Joseph, Joseph, Wiscombe and WeinmanJoseph and others, 1976)), where A = [I _ice(0, μ ₁)−I _ice(0, μ ₂)]/(μ ₁−μ ₂) and B = I _ice(0, μ ₁)−μ ₁ A. This allows the value of I _ice(0, μ ₁) to contribute to α. Then by assuming R(μ)≈0.065, the second integral in the numerator of Equation (10) can be evaluated analytically giving

where R ₀ = 0.065 for isotropic incident radiance and R(μ _in) for a direct incident beam.

In some cases, a two-layer model is needed. To extend the model, four additional boundary conditions must be satisfied at the interface. These are given by

where the index of refraction is the same in both layers.

Model parameters

To complete the model, the following quantities are specified: the phase function, the absorption coefficient, the ice thickness, and either the scattering coefficient or the single-scattering albedo. The total optical depth and ϖ ₀ or σ are then given by Equations (2) and (5).

The ice thickness is measured directly, and absorption coefficients for ice are taken from the results of Reference SaubererSauberer (1950) as reported by Reference GoodrichGoodrich (1970). For the phase function of sea ice, the experimental data of T. C. Grenfell and D. Hedrick (personal communication) are used. Curve 1 in Figure 3 is a smoothed fit of their preliminary results. A large backscattering component is present compared with that due to air bubbles in water (Reference DavisDavis, 1955) or ice. For mathematical convenience, a Henyey–Greenstein phase function (Reference Henyey and Greenstein.Henyey and Greenstein, 1941) is used to represent p(Θ). To include the backscattering component, a three-parameter function of the form

is adopted, where the second term contributes the back scattering (g ₂ < 0), and γ is the relative weighting factor between two terms. A non-linear regression analysis was used to determine the coefficients. Best-fit values were found to be γ = 0.279, g ₁ = 0.71, and g ₂ = −0.25, and the resulting phase function is plotted as curve 2 in Figure 3. To show how accurately p(Θ) is represented by the four-term expansion, the associated ϖ_l s were determined from Equation (4) by numerical integration and substituted into Equation (3). This derived phase function, curve 3 in Figure 3, deviates from the shape of the original curve most strongly at large scattering angles. But since most of the scattering in the four-stream model involves deflections of 120° or less, the approximation is quite good, and errors in the irradiance are not large (Reference LiouLiou, 1974).

Fig. 3.

Representations of the phase function used for young sea ice. Curve 1 is the experimental results of Grenfell and Hedrick. Curve 2 is the best fit of curve 1 using a three-parameter Henyey–Greenstein phase function. Curve 3 displays the four-stream approximation to curve 2.

Response of the model to variations in the basic parameters

A range of cases was examined to investigate the effect of the single-scattering albedo and the phase function on the calculated albedo and transmittance. For a particular ice thickness, increasing ϖ ₀ gives an increase in albedo and a decrease in transmittance. As ϖ ₀ approaches one, albedo and transmittance become very sensitive to small changes in ϖ ₀. Calculations for a variety of phase functions indicate that albedo decreases and transmittance increases as p(Θ) becomes more forward-peaked.

Albedos and transmittances, both neglecting and considering refraction at the air–ice interface, were calculated for direct-beam and diffuse incident radiation. The zenith angle of the direct beam was taken to be 65°, a representative value for Arctic latitudes in the summer. The results are given in Table I. They show that refraction causes a decrease in albedo and an increase in transmittance. This results from the following effects: A portion of the upwelling radiance at the surface undergoes total internal reflection back into the ice, and much of this radiation is then absorbed by the ice or transmitted to the water. Also, refraction converges the incident radiation into a more normal direction, and, because the scattering is primarily forward-directed, a larger fraction of the radiation is transmitted to the water.

Table I.

Effects of refraction on albedo and transmittance

Comparison of theoretical and experimental albedos and transmittances

Comparisons of four-stream calculations with experimental results were performed primarily at 650 nm. This wavelength was nearly optimum for the spectrophotometer used in the observations (Reference Roulet, Roulet, Maykut and GrenfellRoulet and others, 1974) since it combines a large signal-to-noise ratio with high spectral resolution. This also made it possible to use the phase functions of T. C. Grenfell and D. Hedrick (personal communication), which were measured with a helium neon laser. Since an optical diffuser was used during the laboratory experiments, an isotropic incident radiation field was assumed for all calculations pertaining to the observational results.

For a particular experimental case, ϖ ₀ was selected so that the albedo calculated using the four-stream model would best match the experimental albedo. In order to compensate for the effects of the spectrophotometer’s bandwidth, a corrected absorption coefficient of 0.42 m⁻¹ was used rather than the standard value at 650 nm of 0.32 m⁻¹ (Reference SaubererSauberer, 1950). The corrected value was determined by convolving a representative spectral transmission curve with the bandpass of the instrument.

This technique was applied to data from ten experimental cases and the results summarized in Table II. Albedos for cases 3 through 10 were matched with values of ϖ ₀ ranging from 0.955 to 0.988. The difference between theoretical and experimental transmittances ranges from 0.6 to 41% and in all cases but one, the theoretical transmittance is larger than the experimental. The differences are smallest for warm, slowly grown ice and largest for the very cold, rapidly grown ice. Best agreement is achieved for cases 7, 8, and 10, where the ice was most nearly homogeneous.

Table II.

Comparison of one-layer model with experimental results

The assumption of homogeneous ice is probably the largest source of error in the theoretical calculations. Temperature gradients, changes in growth rate with ice thickness, and differences in crystal structure all cause vertical variations in the physical state of the ice which can affect the optical properties. In addition, changes in the ice surface strongly influence the albedo and transmittance. This will be discussed in detail in the next section.

For cases I and 2, which were most poorly represented by the one-layer analysis, a two-layer model was applied. These cases had the largest temperature gradients, and the surface temperature of the ice was at or below the eutectic point. In case 1, for example, the top two centimeters of ice were below the eutectic point for NaCl. This resulted in a region containing precipitated salts whose optical properties were distinctly different from those of the underlying ice.

The two-layer calculations which best represent cases 1 and 2 are summarized in Table III. While the agreement between experimental and theoretical results is improved, substantial differences still remain suggesting that for cold ice a multilayer model may be necessary. In addition, for case 1 it was not possible to fit the albedo precisely. This was probably due to the presence of patches of highly scattering NaCl crystals distributed unevenly across the surface of the ice. By assuming refraction for these patches their enhancement of the albedo is underestimated. To match the observed albedo, considerable reformulation of the model would be needed. However, a more detailed analysis of the experimental results is not warranted until more comprehensive data concerning the relationships of the single-scattering albedo and the phase function to the state and structure of the ice are available.

Table III.

Comparison of two-layer model with experimental results for cold ice.

For several cases spectral albedos and transmittances were calculated at 50 nm intervals between 400 and 1000 nm. The spectral dependence of k _λ for pure ice was determined from the results of Reference SaubererSauberer (1950) in the same fashion as k ₆₅₀. Since scattering inhomogeneitics in the ice are much larger than the wavelengths involved, the limit of geometrical optics applies for scattering by sea ice. Reference Ellison and PeetzEllison and Peetz (1959) show that in this limit the wavelength dependence of scattering by spheres is determined by the index of refraction and is not influenced by particle size. This result also applies to nonspherical scatterers. Since the index of refraction is a weak function of wavelength from 400 to 1000 nm, the scattering coefficient and the phase function are assumed to be independent of wavelength in this range.

Ice 21.1 cm thick grown at an air temperature of −10° C (case 7) was selected as the most suitable case for comparison since the agreement at 650 nm was best. A scattering coefficient of 11.9 m⁻¹ was determined from the value of ϖ ₀ at 650 nm given in Table II. Theoretical spectral albedos are plotted together with the experimental values in Figure 4. The agreement between the two curves is quite good at all wavelengths; however, the slope of the theoretical curve is slightly greater, and calculated albedos lie above observed values from 400 to 650 nm and below them from 650 to 1000 nm. A comparison of the corresponding transmittances is shown in Figure 5. Due to the previously mentioned instrumental limitations, the observations are only given from 500 to 800 nm. In this region the agreement with observations is acceptable, but again the theoretical curve drops off more rapidly with wavelength, crossing the observational curve at 650 nm. The reason for this behavior is not clear at present, but may be due in part to uncertainties in the absorption coefficients of pure ice. The same analysis was performed for cases 8 and 10. As indicated in Table II. the absolute agreement was not as good, but the results show the same relative behavior for both albedos and transmittances.

Fig. 4.

Comparison of spectral albedos for 0.21 m thick ice at an air temperature of —10° C:______________experimental,-----theoretical. Theoretical calculations used σ = 11.9 m⁻¹ and a diffuse incident radiation field.

Fig. 5.

Comparison of spectral transmittances for 0.21 m thick ice at an air temperature of –10° C: _______experimental,-----theoretical. Theoretical calculations used σ = 11.9m⁻¹ and a diffuse incident radiation field.

Further applications of the model

The effect of the angle of incidence on albedo and transmittance of a 25 cm thick ice slab for clear sky conditions is investigated. Clear sky conditions are approximated by representing the incident radiation field as a direct beam at θ _in. From 0 to 60° albedo and transmittance vary gradually with zenith angle, but from 60° to 90° they are very sensitive to changes in θ_in. This is a result of the increase in the Fresnel reflection of the direct beam and the forward peak of the phase function. For comparison α and T were computed for cloudy skies. At approximately 54° results for direct-beam incident radiation are equal to those for diffuse incident radiation in all three cases (Fig. 6).

Fig. 6.

Albedo_____and transmittance----as a function of the zenith angle of a direct incident beam for: (1) ϖ ₀ = 0.97, (2) ϖ ₀ = 0.98, and (3) ϖo = 0.99. The values for isotropic incident radiation are given by the large dots, •, and are plotted on the curves to give the effective zenith angle.

The magnitude and spectral dependence of albedo and transmittance are markedly influenced by ice thickness. Using a representative scattering coefficient for growing ice at −20° C of 15.7 m⁻¹, spectral albedos and transmittances have been calculated for ice ranging in thickness from 5 to 25 cm and are shown in Figures 7 and 8. As the ice thickens, the increase in albedo is greatest at shorter wavelengths where ϖ ₀ is largest. Near 1000 nm no increase occurs for ice thicker than 5 to 7 cm since the ice is already optically thick. Radiation which penetrates below these depths and is then backscattered is almost completely absorbed before it can emerge from the ice. As the thickness approaches zero, the albedo approaches the limit of Fresnel reflection for water (0.065). The transmittance is strongly dependent on wavelength, decreasing sharply in the infrared. As the ice thickness increases, the radiation near 1000 nm is quickly absorbed by the ice and when the ice is 15 cm thick, the transmission is essentially zero. In the visible the attenuation is much more gradual for 25 cm ice and the transmittance is reduced by only a factor of 2.

Fig. 7.

Spectral albedos calculated for ice ranging in thickness from 0.025 to 0.25 m using a scattering coefficient of 15.7 m⁻¹ for a diffuse incident radiation field. The curve labels indicate the ice thickness in centimeters.

Fig. 8.

Spectral transmittances calculated for ice ranging in thickness from 0.025 to 0.25 m using a scattering coefficient of 15.7 m⁻¹ for diffuse incident radiation. The curve labels indicate the ice thickness in centimeters.

During the laboratory experiments, the ice frequently developed a thin surface layer which caused a visible change in its appearance. Two important cases of this were encountered, one occurred when the temperature of the surface dropped below the eutectic point precipitating out solid salts within the ice, and the other when the ice was warm and water began to collect on the surface. To give a quantitative estimate of the effects of these layers on the spectral albedo, two layer models are constructed, again assuming σ = 15.7 m⁻¹ for the underlying ice. In the surface layer, scattering coefficients of 400 m⁻¹ and 6.8 ⁻¹are chosen based on cases 1 and 10 in Table II for very cold ice and melting ice respectively. Spectral albedos are then calculated for ice 25 cm thick including 1 and 5 cm surface layers of both ice types and for homogeneous ice without a surface layer. The results are plotted in Figure 9.

Fig. 9.

Effects of different surface layers on spectral albedos for diffuse incident radiation. The total ice thickness is 0.25 m and the scattering coefficient for the underlying ice is 15.7 m⁻¹: (1) melting surface layer 50 mm thick (σ = 6.8 m⁻¹); (2) melting surface layer 10 mm thick (σ= 6.8 m⁻¹); (3) no surface layer; (4) highly scattering surface layer 10 mm thick (σ = 400m⁻¹); and (5) highly scattering surface layer 50 mm thick (σ = 400 m⁻¹). Diffuse incident radiation was assumed for all calculations.

The presence of a cold and highly scattering layer increases the albedo considerably and alters its spectral distribution. For a 1 cm surface layer (curve 4) the increase is about 29% at 400 nm rising to 250% at 1000 nm. The albedo is nearly constant across the visible region and the ice appears white. A 5 cm layer (curve 5) further increases the albedo at visible wavelengths by 25% almost uniformly from 400 to 650 nm. In the infrared, however, the increase is smaller and at 1000 nm amounts to only about 8%. This saturation effect at long wavelengths, also observed in Figure 7, arises because the thickness of ice contributing significantly to the albedo decreases as k _λ , and hence wavelength, increases.

The presence of a thin melting layer has a much smaller influence. A 1 cm layer decreases the albedo uniformly at all wavelengths (curve 2) by about 0.01. Increasing the layer thickness from 1 to 5 cm (curve 1) further decreases α_λ by 0.025 but does not significantly alter the shape of the curve.