Introduction

If one shines a light on a snowball or a “bubbly” ice cube, a large fraction of the incident light is scattered out through the sides. Scattering is also seen in the various kinds of ice which form on northern lakes and rivers. If the ice contains very few bubbles, solar radiation is specularly reflected. In directions away from the sun, little light is reflected back to the observer, and the ice appears black. However, if the ice contains many bubbles, or if snow covers it, the solar beam is scattered in all directions, and the surface appears almost uniformly bright. Clearly, scattering should be important in distributing light in layers of snow and “bubbly” ice.

Quantitative reasons for considering multiple scattering

The qualitative observations above may be supplemented by more precise measurements. To describe the behavior of light, we introduce the specific intensity I (z, μ, ϕ). The basic definition of I may be found in a number of books, e.g. Reference ChandrasekharChandrasekhar (1950) or Reference GillhamGillham (1970). The arguments of I are as follows: z is the vertical distance from the top of the ice or snow, μ is cos θ, where (θ, ϕ) are polar coordinates describing the direction in which the light is travelling, and θ = 0 is in the downward (+z) direction.

The downward flux Φ↓(z) measured by a photocell or similar instrument imbedded in the ice or snow approximates the flux

(1)

Similarly, the total flux is

(2)

If a thick layer of ice contains no bubbles, the intensity of the reflected solar beam and the intensity of the beam in the ice are simply related to the incident solar intensity. Within this homogeneous ice layer, the specific intensity will obey the familiar Bouguer–Lambert law:

(3)

where a is the albedo of the ice, Φ ₀ is the solar constant above the ice, n is the index of refraction, κ is the absorption coefficient, ϕ ₀ is the azimuth of the sun, and μ ₀ is cos θ ₀. θ ₀ is related to the solar altitude by Snell’s law. The flux of solar radiation within the ice will be

(4)

The constants a and n are required to give conservation of energy at the surface of the ice. If the sensor within the ice is surrounded by a hemispherical bubble of air, then we must multiply the result in Equation (4) by n ² to derive the measured flux. Clearly, once we have measured the index of refraction and the extinction coefficient in the laboratory, we know the distribution of solar radiation in clear ice. Since n and κ are fairly independent of pressure and temperature over the range of these variables encountered in meteorological work, the radiation distribution in and the albedo of clear ice should also be independent of the weather.

This is not true of the properties of snow or “bubbly” ice. For example, snow albedos range from 0.90 for fine-grained, wind-packed snow to 0.45 for old, wet snow (Reference LiljequistLiljequist, 1956). Snow albedos are even reported to change from morning to afternoon, for the same solar altitude (Reference HubleyHubley, 1955; Reference LiljequistLiljequist, 1956). Large variations of albedo are to be expected for snow or “bubbly” ice if their scattering properties depend strongly upon temperature or pressure. For example, we might expect a considerable difference in the distribution of light scattered by a particle of ice covered with a thin layer of water, as compared with the distribution of light scattered by a similar dry particle. Also, the size and shape of ice grains and air bubbles should have a strong effect on the scattering properties of the medium.

There are other difficulties in applying Equation (3) or Equation (4) to snow or “bubbly’ ice. The extinction coefficients of these media are large near the surface and decrease to a constant value larger than that of clear ice, as may be seen in Reference Ambach and MockerAmbach and Mocker (1959) or Reference LyubomirovaLyubomirova (1962). If interpreted strictly, an extinction coefficient that depends on depth indicates that radiation is not extinguished exponentially. Attempts to explain this phenomenon on the basis of selective absorption, e.g. Reference LiljequistLiljequist (1956), are not entirely convincing. Reference LyubomirovaLyubomirova (1962) has suggested that multiple scattering will also lead to a non-exponential flux decrease near the surface and an exponential flux decrease deeper in the medium. As we will show, it should be possible to distinguish observationally between selective absorption and multiple scattering by the dependence of the flux extinction on solar altitude.

In addition to a variable extinction coefficient, snow, and presumably “bubbly” ice, show a sizeable upward flux. Reference LiljequistLiljequist (1956) suggests that his measurements show that the upward flux is about 40% of the downward flux. It is difficult to explain this observation without resorting to multiple scattering.

Despite the apparent importance of multiple scattering there appears to be little work which has used modern techniques in radiative transfer. Reference Dunkle and GierDunkle and Gier (1953) and Reference Dunkle and BevansDunkle and Bevans (1956) have formulated the problem in terms of a two-stream model. The former find that in this approximation, the ratio of upward to downward flux is a constant determined by a distributed reflectivity and the absorption coefficient. Reference LiljequistLiljequist (1956) has modified Reference Dunkle and GierDunkle and Gier’s (1953) approach near the upper boundary by introducing an ad hoc surface reflectivity. Reference LyubomirovaLyubomirova (1962) assumes that the scattered light is proportional to the diminishing solar flux. She shows from observations that the scattered radiation builds up beneath the surface and then decreases. We will formulate the scattering problem in terms of conventional radiative transfer theory.

Mathematical formulation of the scattering problem

On a scale of a few tenths of a millimeter, both snow and “bubbly” ice are extremely heterogeneous materials. Thus, we might think of snow as small, irregular particles of ice imbedded in air, and of “bubbly” ice as small parcels of air imbedded in ice. For transfer calculations, we could compute the way in which light is scattered by an “average” particle or parcel, and then use the results to compute the effects of multiple scattering by a large number of similar particles. However, because the irregular shape and random orientation of either ice grains or air bubbles make such an approach very complicated, we will use a phenomenological approach to radiative transfer to gain some idea of the effects of multiple scattering.

In the phenomenological approach, we note that snow or “bubbly” ice often has homogeneous properties on the scale of a few millimeters. We can then follow Reference ChandrasekharChandrasekhar’s (1950) approach to the optical properties of the medium. We assume that the snow or ice can be characterized by a volumetric scattering coefficient σ (m⁻¹), From a beam of radiation incident normally on a thin slab of volume dV having a cross-section dΣ and a thickness ds, the power scattered is

(5)

dω is a small element of solid angle about the normal to the surface of the slab. The power scattered into an element of solid angle dω′ about a direction inclined at an angle Θ to the direction of incidence of dV is

(6)

where p (cos Θ) is the phase function. The total power within the beam scattered by dV is

(7)

so that

(8)

Similarly, the total power absorbed in dV from the beam in dω is

(9)

where κ (in m⁻¹) is the (true) absorption coefficient.

The scattering and absorption coefficients are assumed to depend on the physical condition of the ice and snow. However, κ, σ and p do not completely specify the scattering problem. In some cases, it is useful to have sources of radiation which add to the beam within a solid angle dω a power

(10)

For example, if we wish to consider only the diffuse radiation within the snow or ice, we can use a source q which describes the radiation which is directly scattered from the solar beam. Similarly, we could describe by q the radiation which is directly scattered from the “beam” of diffuse solar radiation that penetrates the medium. Finally, we must specify the boundary conditions on I.

If we assume that the snow or ice layer can usually be treated as a plane-parallel medium, then the changes in I are described by the equation of transfer:

(11)

Although Equation (11) can be written so that it applies to each frequency, we assume that I, κ, σ, p and q are independent of frequency. We do so because the sunlight reflected from snow or ice is essentially uncolored, i.e. they are grey media.

Our particular interest lies in the distribution of solar radiation in a thick layer of snow or “bubbly” ice. The properties of ice and snow are often macroscopically homogeneous over large depths. The assumption of homogeneity allows us to define a new depth variable—the optical depth—

(12)

where χ (also measured in m⁻¹) is called the extinction coefficient. The constant

(13)

is the single-scattering albedo. As we have already noted, little is known about the nature of scattering in snow or ice, i.e. p is unknown. We still therefore take the simplest possible choice:

(14)

The boundary conditions for the problem are complicated by refraction. We included refraction in the treatment of solar radiation in clear ice. However, we have not included it explicitly in our derivation of Equation (11) because it is almost impossible to measure the index in translucent media. The path of a light beam is very difficulty to follow far into the medium, so Snell’s Law cannot be used. We shall continue to neglect effects of refraction, even at the surface of the medium. No special surface reflectivity is to be invoked. The grains of snow on the surface scatter light in the same way as the grains deeper in the layer; the bubbles which pit the surface of a homogeneous layer of “bubbly” ice scatter light in essentially the same way as the bubbles deeper in the layer. In other words, light is not specularly reflected from the surface of the ice or the snow. We assume that the direct solar beam entirely penetrates the medium. The outward flux comes from the scattered radiation. Since the diffuse solar radiation under clear skies usually contributes ten to twenty per cent of the global solar flux (Reference LiljequistLiljequist, 1956; Reference RusinRusin, 1961), we neglect the contribution of the diffuse solar radiation to the distribution of solar radiation in the snow or ice. Finally, if the snow or ice layer is fairly thick, the solar radiation is almost completely extinguished at the bottom of the layer. This allows us to neglect the light absorbed or reflected by the underlying medium, and to treat a deep layer of snow or “bubbly” ice as a semi-infinite, plane-parallel “atmosphere”.

These abstractions allow us to simplify the transport problem considerably. The radiation scattered from the direct solar beam will obey the equation

(15)

where I is subject to the boundary conditions

(16a)

and

(16b)

Equation (15) agrees with equation (126) in Chapter 1 of Reference ChandrasekharChandrasekhar (1950), except that our μ. is measured in the opposite direction from his.

If we replace the incident solar beam in Equation (15) by the axi-symmetric incident beam (Φ ₀/2π) δ(μ – μ ₀), the diffuse solar radiation in the scattering medium still obeys Equations (15) and (16a, b). In addition, the flux and mean intensity of the new incident beam are identical with those of the original beam. Therefore, the solar radiation in snow and “bubbly” ice should be essentially the same as the distribution of the specific intensity

(17)

which satisfies the homogeneous transfer equation

(18)

and has the boundary conditions

(19a)

and

(19b)

Solution for the solar radiation field

The problem set down in Equations (18) and (19a–b) can be solved by using the singular eigenfunction method. The solution may be found in Reference Case and ZweifelCase and Zweifel (1967) in the form

(20)

where the expansion coefficients are

(21)

and

(22)

with

(22)

P being the mnemonic symbol to remind us to take the principal value of the integral in Equation (20),

and

X(–v) is the X-function tabulated in Appendix L of Reference Case and ZweifelCase and Zweifel (1967) (after Mendelson, unpublished) for various values of ϖ. The functions ϕ ₀₊(μ) and ϕ_ν (μ) are the eigenfunctions for the homogeneous transfer equation with respective eigenvalues v ₀ and v. ϕ ₀₊ and ϕ _v are normalized in the sense

(23)

The total solar flux within the snow or ice is related to the specific intensity by

(24)

Integrating Equation (18) over μ, and multiplying by 2π, we find

(25)

Using the normalization condition, Equation (23), the solar flux divergence can therefore be written

(26)

Properties of the solution

The solution to the transport problem has several interesting features. First, the specific intensity decreases non-exponentially near the surface, and exponentially deep in the medium. This property may be seen by noting that v ₀ ⁻¹ ⩽ 1 ⩽ v ⁻¹ for 0 ⩽v ⩽ 1, Thus, we can consider the integral in Equation (20) to be a boundary transient, and the discrete term to a characteristic or asymptotic solution.

The asymptotic solution might be considered a modified “Bouguer–Lambert Law”. However, the distribution of light deep in the scattering medium does not depend inversely upon the sine of the solar altitude. If we follow the procedure described in Reference LyubomirovaLyubomirova (1962) for measuring the “extinction coefficient” K _c, clear ice should give K _c = κ/μ ₀ in accordance with Equation (4). However, the “asymptotic extinction coefficient” in Equation (20) would be K _s = (κ + σ)/v ₀. Clearly, K _s is independent of solar altitude. So far as the author knows, there are no observations to decide whether or not the extinction coefficient in snow or translucent ice depends upon μ ₀.

If we take the asymptotic values of the measured extinction coefficients, we can derive a value for χ which will allow us to relate τ to z. Using Reference LiljequistLiljequist’s (1956) mean value for K _s of 0. 145 cm⁻¹ we find

(27)

In deriving this number, we used a value for ν ₀ which corresponds to ϖ = 0.990. This high single scattering albedo is based on the observed albedo, in a manner which will be discussed shortly. From Equation (27), we see that τ = 1 corresponds to about 0.85 cm, τ = 3 to about 2.55 cm. The absorption coefficient is related to χ according to

(28)

The interesting point about this value of κ is that it is slightly smaller than the value of 2 m⁻¹ observed for normally incident light on glass-like river ice (Reference LyubomirovaLyubomirova, 1962). The theoretical value of κ may be slightly larger if ϖ is smaller. In view of the possible experimental error, the agreement between theory and experiment is reasonable. This result suggests that the absorption occurs in the ice of the scattering media.

In addition to having an extinction coefficient that is independent of depth, the asymptotic distribution of solar intensity is independent of depth. It is given by

(29)

Clearly, I _AS becomes more isotropic as v ₀ (or ϖ) increases. This leads us to expect an increase in the ratio of upward to downward flux as ϖ increases. From Equation (29) we can show that

(30)

Again, this ratio is independent of τ and μ ₀. Table I shows a a _AS for several values of ϖ. As expected, a _AS increases with v ₀. Reference LiljequistLiljequist’s (1956) observation that a _AS in Antarctic snow is about 0.43 implies that ϖ is approximately 0.9. As we will show, such a value of ϖ is too low to give good agreement with the observed albedo. It is not clear whether the discrepancy is due to the theory or observational error. Apparently, the observations are not easy to make.

Table I. Asymptotic Ratio of Upward to Downward Flux as a Function of Single Scattering Albedo

Since we have an expression for the solar radiation within the snow or ice, we can find the amount of radiation which emerges from the upper surface. Thus, we expect that the albedo of snow or “bubbly” ice in direct sunlight is

(31)

From certain identities involving the X-function, Reference Case and ZweifelCase and Zweifel (1967) obtain an expression equivalent to

(32)

By using the identity

(33)

and the related fact from Reference McCormick and KusčerMcCormick and Kusčer (1965) that

(34)

it is possible to show that

(35)

Clearly, a is a function of μ ₀ and ϖ. From values of v ₀ ⁻¹ given by Reference CaseCase and others (1953) and the tables of the X-function given in Reference Case and ZweifelCase and Zweifel (1967) we can find the albedo as a function of μ ₀ and ϖ. The results are shown in Table II.

Table II. Albedo as a Function of the Angle of Incidence and the Single Scattering Albedo

X (–μ ₀) depends very little upon ϖ for the values of ϖ shown in Table II. By rearranging Equation (35), we find

(36)

Thus, if we have observations of the albedo for snow or “bubbly” ice under clear skies, we can compute approximate values of v ₀ by using the tables of X( – /μ ₀) for an approximate value of ϖ. Some Antarctic expeditions have this kind of data, e.g. Reference LiljequistLiljequist (1956); Reference RusinRusin (1961). Figure 1 shows the observed and computed values of the albedo as functions of μ ₀ (or sine of the solar altitude). Reference RusinRusin’s (1961) observations are those for station Mirny (lat. 66° 33′ S., long. 93° 01′ E.). His albedos are given to two significant figures, and this uncertainty is reflected in the bars of Figure 1. The observations of Reference LiljequistLiljequist (1956) are averages for clear skies in November 1951, and in December 1951–January 1952, at Maudheim (lat. 71° 03′ S., long. 10^o 56′ W.). The November observations are about 0.015 higher than those of December–January. No uncertainties are shown for Liljequisi’s data because his table gives the albedos to three figures.

Fig. 1. Plot of the clear-sky albedo (a) versus the sine of the solar altitude (μ ₀). Data are from Antarctic expeditions: and Liljequist and Reference RusinRusin (1961). Upper data (Δ)—average clear-sky albedo for November 1951. Middle data (□)— average clear-sky albedo for December 1951/January 1952. Theoretical curves were computed from Equation (35) using the X-function for ϖ = 0.990, and values of ν ₀ corresponding to ϖ = 0.990 (–––––––––), 0.998 (–––– – – ), and 0.983 (– – – – –).

A value of v ₀ was determined for each solar altitude, using the data of Rusin and Liljequist and interpolating in μ ₀ from the table of X( – μ ₀) for ϖ = 0.990. The mean values are v ₀ = 5.7±0.3 for the November data, v ₀ = 5.3±0.2 for the December–January data, and v ₀ = 4.5±0.1 for Rusin’s data. The uncertainties listed are standard deviations. Rusin’s albedo at 5° gave a value of v ₀ more than three standard deviations from the mean, so it was not included in computing the mean. The theoretical curves shown in Figure 1 were computed using Equation (35) and X(–μ ₀) for ϖ = 0.99. The values of v ₀ used in Figure 1 were 5.797, 5.289 and 4.407, in descending order. The corresponding values of ϖ are 0.990, 0.988 and 0.983. The theoretical curve for v ₀ = 5.797 was used rather than that for v ₀ = 5.713, because the former corresponds to ϖ = 0.99, for which the tabulated X-function is exact.

The observed albedos are for clear skies, so that direct solar radiation should be the predominant source of radiation in the snow. However, some light is scattered by the atmosphere, contaminating the direct solar beam with light scattered from altitudes lower than the sun’s. Since the albedo is highest for the most obliquely incident light, the scattered light tends to increase the albedo. This effect becomes more important at lower solar altitudes because there is more scattered light (Reference LiljequistLiljequist, 1956; Reference RusinRusin, 1961). In addition, there may be observational errors connected with measuring the direct solar flux (Reference FritzFritz, 1958). Thus, the observed albedos at low solar altitudes are more uncertain than those at high solar altitudes. Nonetheless, the computed albedos agree with the observed albedos to within 1%, except at μ ₀ = 0.08.

The increase of albedo with decrease of solar altitude also appears to be in qualitative agreement with Reference LanglebenLangleben’s (1968) ice albedo measurements and Reference HubleyHubley’s (1955) glacier albedo measurements. The latter suggests that powder snow has an albedo very near 1.00 independent of solar altitude. This is very interesting in the light of Equation (35). It would appear that the very small grains of powder snow have a single scattering albedo approaching 1.0 (v ₀ → ∞). Reference HubleyHubley (1955) also suggests that the albedo of snow on glaciers increases with the increase of light coming from small altitudes. In general, then, our theory seems to be in accord with albedo observations.

The solar flux divergence

Once we know the value of ϖ which applies to the medium in which we are interested, we can compute the radiation field within it. Two of the most important moments of the field are the flux, Φ(τ) and the mean intensity:

(37)

For the solar radiation in the snow or “bubbly” ice, these two moments are related according to a variant of Equation (25):

(38)

Therefore, if we have Φ, we can find J, and vice versa.

The flux divergence is also important because it is directly related to the heating of the medium by solar radiation. If we neglect conduction, convection, and latent-heat transfer within a snow or ice layer, then the temperature in the layer will change according to the equation

(39)

where T(z, t) is the temperature at a time t and a depth z, ρ is the density (kg m⁻³), c is the specific heat of the material (kJ kg⁻¹ deg⁻¹), Φ is the total flux of solar radiation (kW m⁻²), Φ _L, is the total flux of long-wave radiation (kW m⁻²), l is the latent heat of melting (334 kJ kg⁻¹), and m is the rate of melting (kg m⁻³ s⁻¹).

m is zero as long as T is less than 0° C, and the temperature in the layer then changes in accordance with Equation (39). Once the melting point is reached, the temperature remains at 0° C, and we have the steady-state condition

(40)

Near the top of the layer, the long-wave flux may be important. For example, if the temperaure of the air above the layer is below 0° C, long-wave radiation will be emitted by the top few centimeters of the layer in an attempt to heat the air, and the melting rate will be decreased. If we neglect the long-wave flux, for purposes of simplification, then the melting rate will be given by

(41)

Using Φ ₀ = 1.000 kW m⁻², χ = 85.1 m⁻¹, and l = 334 KJ kg⁻¹, we find that

(42)

The total melting rate (kg m⁻² s⁻¹) for any finite layer may be found by integrating m over the proper limits of depth. Thus, if we wish to find the total down-wasting of the layer, we must integrate m over the whole layer. Since the flux vanishes at the lower boundary, the number of grams of water per square meter per second produced by the direct solar radiation is (1 – a) μ ₀ Φ ₀ l ⁻¹.

Because the flux divergence is so important, we have computed it for a scattering medium with ϖ = 0.990. The results are shown in Figure 2. A definite region of non-exponential flux divergence is evident near the surface. If the melting rate is proportional to the flux divergence, it is clear from Figure 2 that the snow or ice will melt fastest just beneath the surface when the sun is high in the sky. As the sun sets, the fastest melting moves to the surface. The calculations show that the flux divergence is exponential beyond τ ≈ 3. On die basis of the measurements quoted previously, this depth corresponds to about 2.6 cm. Clearly, all measurements of the asymptotic flux extinction coefficient should be made well below this level.

Fig. 2. Plot of the flux divergence ( − d{Φ/Φ ₀}/dτ) versus the optical depth (τ) for various solar altitudes. The theoretical values of the flux divergence were computed using ϖ = 0.990. τ = 2.0 corresponds to about 1.7 cm for Liljequist’s (1956) data. This plot is also proportional to the melting rate under some conditions (see Equation (42)).

In some cases we may wish to melt snow more rapidly than would occur naturally. One way to accomplish this end is to scatter an absorber on top of the snow, decreasing the albedo, the flux divergence and albedo calculations suggest that a similar result could be achieved by scattering the absorber while the snow is falling. Since the absorption is then dispersed throughout the layer, less melt water should be evaporated than is the case when the absorber is on top of the layer.

Suggestions for further investigation

As we mentioned earlier, anisotropic scattering needs to be included in a definitive theory of multiple scattering in snow and “bubbly” ice. The ice grains in snow and the air bubbles in ice probably scatter radiation in a manner very similar to that of water droplets in air. The latter are known to scatter radiation very anisotropically. The observed albedo could be produced by any of a large number of combinations of single scattering albedo and phase function, so our determination of ϖ is not unique. Constraints might be placed on such combinations by either directly measuring or computing ϖ and p. In some cases of anisotropic scattering, the singular eigenfunction solutions to the transport equation can be extended (Reference Case and ZweifelCase and Zweifel, 1967; Reference MikaMika, 1961). Depending on the degrees of anisotropy, there are one or more discrete eigenvalues. We still anticipate that the flux divergence will be non-exponential near the surface, and exponential deep in the medium. The theoretical results will probably also need to be extended to include finite layers of snow.

There also appear to be a large number of observations which could improve our understanding of multiple scattering in snow and ice. First, it would be very helpful to have observations of the azimuthal distribution of scattered light, rather than just albedo measurements. As shown by the observations of Reference Salomonson and MarlattSalomonson and Marlatt (1968), snow tends to forward scatter. With an extension of the theory, such azimuthal observations could be used to determine the phase function. Second, measurements of the flux extinction in different colors would be very helpful in deciding what role is played by selective extinction. The theory presented here could lend some help in interpreting the results. Third, the relation between ϖ and ρ could be determined observationally for different kinds of snow and ice. Hopefully, this relationship would simplify the relationship between the physical properties of the media and their albedo. Finally, measurements of the ratio of upward to downward flux would be of great interest in deciding whether or not the single scattering albedo measured by a _AS agrees with the single scattering albedo deduced from the albedo. However, even at this stage, multiple scattering in snow and ice can tie together many disparate observations.