Introduction

In a previous paper by one of the authors (Reference BarkstromBarkstrom, 1972, hereafter referred to as paper I) it was shown that multiple scattering in snow could be modelled with considerable success by solving the phenomenological equation of transfer with isotropic scattering. However, as Reference Middleton and MungallMiddleton and Mungall (1952) and Reference Salomonson and MarlattSalomonson and Marlatt (1968) have shown, the light “reflected” from snow is directed in the forward scattering direction, particularly for low solar altitudes. This phenomenon cannot be explained by isotropic scattering. The purpose of this paper is to refine the results in I by modelling the optical properties of snow including the effects of anisotropic scattering in finite layers. In this model the transport equation is solved by an algorithm that gives the optical properties of a multiply scattering layer composed of two stacked layers of known properties. With this algorithm it is possible to infer a phase function that will result in a layer whose surface albedo and angular reflectance is in good accord with experimental measurements. The modelling procedure is sufficiently general that vertically inhomogeneous layers can be constructed which closely simulate stratified layers with different scattering properties. Moreover, the reflectance properties of the “ground” underlying the “snow” layer can be realistically included.

Numerical results of the modelling procedure are given that compare well with Reference Middleton and MungallMiddleton and Mungall’s (1952) measurements on “new snow fallen in calm”. The computations indicate that inclusion of anisotropy causes the flux divergence to maximize deeper in the snow than for isotropic scattering (in terms of optical depth) and that layers of optical depth less than five will exhibit a flux “extinction coefficient” three to five times that observed deep within “semi-infinite” layers (for the same extinction coefficient). The success of the anisotropically scattering model in predicting layer albedo and angular surface reflectance for the above snow suggests that it would be both possible and useful to characterize other types of snow layers by “effective” phase functions. Although no simple analytic procedures are found for deducing these phase functions, they are readily formed by comparing observations of surface albedo and angular reflectance with the results of calculations. A physical theory for the phase function will be found elsewhere (Reference Bohren and BarkstromBohren and Barkstrom, 1974).

Formulation of the Problem

In paper I the equation of transfer, equation (I-11) was solved with an isotropic phase function, equation (I-14), to obtain an analytic expression for specific intensity (spectral radiance) I(τ, μ, φ) within and at the surface of a vertically semi-infinite medium. In this paper we use an alternate, but equivalent, method of solving the transfer equation to obtain the specific intensities within and emerging from the top and bottom of a laterally infinite, but vertically finite, plane-parallel slab containing anisotropic scatterers. The method of solution is sufficiently general that a wide variety of scatterers can be incorporated in the slab, so that realistic, vertically inhomogeneous slabs such as snow packs can be modelled fairly easily. The solution algorithm is well known in the study of stellar and planetary atmospheres, but because it is less well known elsewhere and because it appears to be a useful tool in characterizing the optical properties of snow, the algorithm is discussed in some detail in the text and in the Appendix.

As in paper I, I(τ, μ, φ) is the specific intensity of radiation in a laterally infinite, plane-parallel slab at an optical depth τ travelling in the direction specified by the azimuth φ and the cosine of the polar angle μ = cos θ where θ is measured from the downward normal to the slab. The optical depth is the dimensionless product τ = χd of the volume extinction coefficient χ and the vertical distance d. It is measured downward from the top of the slab. It is convenient to set μ = |cos θ| and explicitly prefix μ with the proper sign so that +μ denotes radiation travelling downward and –μ denotes radiation travelling upward. It is usual (but not necessary) to assume that no sources are present within the slab so that the only radiation sources are the intensities I(0, +μ, φ) and Ι(τ ₀, —μ, φ) respectively incident at the top and bottom of a slab of total optical depth τ ₀.

The direct problem is to find the intensities I(τ,±μ, φ) within the slab and the emergent intensities I(0, —μ, φ) and I(τ ₀, +μ, φ) given the properties of the scatterers and the incident intensities. The algorithm used to find the internal and emergent intensities is variously known as the method of addition of layers (Reference SobolevSobolev, 1963), and doubling method (Reference Van de HulstVan de Hulst, 1963), and the star-product algorithm (Redheffer, 1963; Reference Grant and HuntGrant and Hunt, 1969[a],[b]), but it is formally identical to that given by Reference StokesStokes (1862) for the transmission and reflection from a pile of glass plates. The essence of the method is that if the reflection and transmission properties of one or more scattering layers are known, then the same properties of a stack of two such layers can be found in terms of those of the individual layers. Thus it is possible to compute the properties of a multiply scattering layer by repeated doubling of an initial, optically thin, single-scattering layer of known properties. It is also possible to combine successively layers of dissimilar properties to obtain a vertically inhomogeneous stack with known properties. The algorithm also makes it possible to add arbitrary reflecting boundaries within and at the surface of the slab.

The related inverse problem is to infer the properties of the scatterers given the incident and emergent intensities. The solution to the inverse problem is not as straightforward as that of the direct problem. At present there is no well established procedure for inferring scattering properties when intensities incident upon and emergent from a scattering slab are given. The numerical example given later used an iterative procedure which proceeded from initial estimates of the single scattering properties of the slab with isotropic scattering to a model which produced emergent intensities which compared well with experimental measurements. Fortunately the procedure converged quickly to the results given below.

Method of solution: The Addition Algorithm

The addition algorithm assumes that there exist reflection and transmission operators which linearly transform intensities incident from above and below into intensities emerging from the slab. Assume that the slab is of optical depth τ, and that it is illuminated by either collimated or diffuse radiation from above and below. For brevity let the intensity incident on the top of the slab I(0, +μ, φ) be given by I ₀↓ where the down arrow denotes the direction of travel and the subscript denotes the level τ = 0. Let the intensity incident from below be I ₁↑ where the subscript denotes τ = τ _I,. The emergent intensities are

where R, T, R ^∗ and T^∗ are the reflection and transmission operators. When the slab is homogeneous or symmetric about the median plane, R = R ^* and T = T ^∗. When the slab is inhomogeneous or when polarization properties are included, R ≠ R ^∗ and T = T ^∗. In the text which follows, but not in the appendix, it will be assumed that R ^∗ = R and T ^∗ = T.

Suppose that reflection and transmission operators R ₁, and T ₁, are known for a slab of optical depth τ ₁ and that the corresponding operators R ₂ and T ₂ are known for a slab of optical depth τ ₂. If the two slabs are stacked with the slab of optical depth τ, on top, there will be operators R and T for the stack. In terms of these operators the intensities emerging from the stack are

where the subscript 2 denotes the bottom of the stack. These intensities can also be related to the intensities at the (fictitious) interface between the top and bottom slabs

The intensities at the interface are related to one another and to I ₀↓ by

The equations may be solved for explicit expressions for the interface intensities

The inverse operator [I —R ₁ R ₂] is assumed to exis T, and, as is made clear in the appendix, it may be represented by a simple matrix inversion. Substitution of Equations (9) and (10) into Equations (5) and (6) and comparison with Equations (3) and (4) gives the expressions for R and T:

The R and T operators also may be defined in terms of analogous scattering and transmission functions S ⁺ and T ⁺ given by Reference ChandrasekharChandrasekhar (1950). For diffuse incident radiation they are defined by

These operators give only the diffuse emergent radiation. Accordingly a term exp (—τ ₀/μ) I(0, μ ^’, φ ^’) δ _{μ μ} , δ _{φ φ} , must be added to Equation (14) to give that portion of I(0, μ, φ') already travelling in the (+μ, φ) direction which is transmitted and attenuated by the slab. The Kronecker deltas δ _{μ μ} , δ _φ , serve as a reminder that the term is added only for the μ = μ ^’, φ = φ ^’ component of I(0, μ ^’, φ ^’). The definitions for R and T thus become

The analogy between Equations (15) and (16) and Equations (3) and (4) becomes somewhat clearer when the incident radiation is collimated. In that case the incident intensity in the direction (+μ ₀, φ ₀) is

where φ ⁰ is the net flux (spectral irradiance) and δ(μ—μ ₀) and δ(φ—φ ₀) are Dirac delta functions. The emergent intensities are

and the expression for R and T, Equations (15) and (16) remain the same without the integrals.

Although the fundamental procedure for finding R and T from given R ₁, R ₂, T ₁, and T ₂ has been given, explicit equations for the operators R ¹, and T, or R ₂ and T ₂ have not been given. Initial values for R and T operators may be obtained in either of two ways. Since the R and T operators are simply related to the single-scattering phase function in the limit τ → 0. optically thin layers (τ ≈ 10^-6), can be doubled repetitively to obtain R and T operators for thicker layers. Alternatively, integro-differential equations given by Reference ChandrasekharChandrasekhar (1950) can be solved numerically with reasonable ease to reach depths of τ≈ 10⁻³ in one integration step with doubling or addition used thereafter. The method used is partly a matter of taste, but should more properly be dictated by the computational economics of the problem to be solved.

Starting values for R and T are doubled from the phase function p using Reference ChandrasekharChandrasekhar’s (1950) small τ limits

These equations reduce to

for small τ ₀. The factor Ѓ F(τ ₀) is

The procedure for converting S ⁺ and T ⁺ functions for small τ ₀ to the corresponding R and T is evident from Equations (15) and (16), but it is discussed in detail in the appendix. The discussion of the initialization via integro-differential equations is also given in the appendix.

The effects of the albedo of the ground underlying a slab of snow can also be accommodated in the addition algorithm by treating the reflecting layer as if it were a scattering layer. For example, the scattering function S ⁺ for a Lambert surface (a Lambert surface has a diffuse reflectance which is equal in all directions) with albedo λ is 4μμ ₀ λ so that the reflection operator is μ ₀ λ/π. Fresnel reflections a T, for example, an ice-snow interface can be included by using the usual Fresnel reflection factors with the appropriate change in the sign of μ, i.e. μ → — μ upon reflection. Since reflections are not included in the numerical example no further mention of them is made.

Reflectance of a snow cover

Reference Middleton and MungallMiddleton and Mungall (1952) have given measurements of the angular reflectance for several types of snow layers. The measurements share the property that the layer reflectance increasingly diverges from that of a Lambert surface with decreasing solar elevation angles. This is consistent with the anisotropy present in scattering of light by large particles. We suggest that this behavior can be successfully modelled by multiple scattering in a plane parallel layer containing scatterers with an effective (or fictitious) single-scattering albedo

and phase function p.

The single-scattering albedo is the fraction of photons incident on a particle which are scattered rather than absorbed. The angular dependence of the scattering is given by the phase function p (cos θ) which explicitly depends on the scattering angle θ. When the phase function has the normalization

the product mp(cos Θ) dΩ/4π is the probability that a photon interacting with the particle will be scattered through an angle θ into an element of solid angle dΩ = d(cos Θ) dø. The anisotropy of the scattering is characterized by the asymmetry parameter g given by

It is convenient to expand the phase function in a Legendre series

whose coefficients are

in which the P_i(cos Θ) are Legendre polynomials. Equations (25)-(28) give f₀ ≡ 1 and g ≡ f₁/3.

We have attempted to model the angular reflectance data of Reference Middleton and MungallMiddleton and Mungall (1952) shown replotted in Figure 1 for “new snow fallen in calm”. The single-scattering albedo largely determines the level of the curves and is expected (see paper I) to be about 0.99 for isotropic scattering. The shape of the curves is determined by the phase function, with g playing the major role. The angular reflectance for near grazing incidence and reflection (μ, μ ₀ ≈ 0) is sensitive to the higher terms. However, the properties important for the energy balance, the albedo and the flux extinction, are extremely insensitive to f ₂, f ₃., ... (Bohren and Barkstrom, in press; Reference Van de HulstVan de Huls T, 1970). Accordingly, we have searched for

values of

and g that give reasonable agreement between the data and the transfer computations. The higher terms f₂, f₃, …, have been given arbitrary, but physically plausible values.

Fig. 1.

Plot of the angular dependent of the reflected intensity I from snow new-fallen in calm weather. Data are taken from Reference Middleton and MungallMiddleton and Mungall (1952). The azimuthal angle is denoted by ø, with ø = 0° being the direction in which the solar radiation is travelling, ø = 180° being the direction in which light back-scattered from the surface is travelling, μ is the sine oftht altitude at which Ike light is travelling or, equivalently, the cosine of the zenith angle, μ ₀ denotes the sine of the solar altitude.

Our early attempts to reproduce Figure 1 used a Henyey-Greenstein phase function (see for example Reference Van de HulstVan de Huls T, 1970)

which has an infinite Legendre expansion with coefficients

which depend only on l and g. For values of g greater than 0.4 or 0.5 these coefficients decrease very slowly with increasing l and it is necessary to retain many terms in the Legendre expansion of Equation (29) to avoid truncation problems. This makes it necessary to use large matrices for R and T (see the Appendix) to preserve numerical stability and accuracy. It is also possible to ignore truncation problems and use too small an order of quadrature to estimate an appropriate value of g. This procedure suggested that

= 0.94 and g = 0.5 should be approximately correct. A search then began for an alternate, finite phase function with g ≈ 0.5.

Fig. 2.

Plot of the phase function p for single scattering, as a function of cos θ, where θ is the angle through which the light is scattered.

The alternate phase function, shown in Figure 2, has g = 0.506 and was found by a Mie calculation for a polydispersion of spheres (Reference Querfeld, C.W., M. and J. P.Querfeld and others, 1972) distributed according to a zeroth order logarithmic distribution (Reference KerkerKerker, 1969, p. 356) with modal diameter 0.245 μm' width parameter σ₀ = 0.07, and real refractive index n = 1.20 at an incident wavelength of 0.62 μm. The significant Legendre coefficients for this phase function are g.iven in Table I. This phase function with

= 0.996 was used to compute the intensities emerging upward from a slab with τ ₀ = 128, shown in Figure 3. The lower slab boundary was assumed to be completely absorbing, but this has no effect on the intensities at the upper boundary. Apart from small differences in incident angles the computed intensities in Figure 3 are directly comparable to the measured intensities in Figure 1. Comparison of the two figures shows that the calculated intensities agree remarkably well in spite of differences in detail. The agreement is the more remarkable since the phase function in Figure 2 is for a polystyrene latex of very uniform properties and the microstructure of a snow cover is probably chaotic by comparison. The integration of phase functions of scatterers of quite diverse sizes, shapes, refractivities, and orientations yields an effective albedo and phase function of a relatively simple character. It would be interesting to characterize this and other types of snow with g.reater care and perhaps g.reater attention to higher order terms in the phase function. At this point the authors can only ask for more high-quality observations.

Fig. 3.

Plot of the angular dependence of the light intensity I reflected from a very thick layer of anisotropic scatterer with a phase function as shown in Figure 2. Notation is as in Figure 1.

Table I.

Legenfire expansion coefficients of the phase function

The emergent upward intensities in Figure 3 can be integrated (see the Appendix) to obtain a surface albedo A(μ ₀). Figure 4 shows the albedo A(μ _α) computed with the phase function in Table I and Figure 2 using

— 0.996. The isotropic scattering computations from paper I are also included with the Antarctic data of Reference LiljequistLiljequist (1956) and Reference RusinRusin (1961) as well. The anisotropic albedo lies above the isotropic curve because the higher single scattering albedo used (oo = 0.996 compared with 0.99), but is otherwise very similar. Surface albedos calculated with a variety of anisotropic phase functions with near 0.99 reproduce the variation of a with μ ₀ for isotropic scattering to about 0.1%. If measurements of the surface albedo alone are available, it seems reasonable to estimate using isotropic scattering theory. If it is important one can estimate the effects of anisotropy on the asymptotic intensity extinction and flux divergence, noting that an eigenvalue about twice the isotropic discrete eigenvalue will g.ive a reasonable first approximation.

Fig. 4.

Plot of clear-sky albedo A versus the sine of the solar altitude μ ₀. Data are from Antarctic expeditions (Reference LiljequistLiljequis T, 1956; Reference RusinRusin, 1961), as noted in paper I. Solid line is the albedo from a very thick layer of anisotropic scatterer with the phase function shown in Figure 2; dashed line is the albedo from a semi-infinite layer of isotropic scattertr.

Behavior of the Radiation within the snow Layer

It was shown in paper I that the flux divergence in an isotropic scattering atmosphere has a maximum beneath the surface for certain angles of incidence. The flux divergence in an “atmosphere” with a forward-scattering phase function, such as snow, exhibits a similar behavior, except that the maximum of the flux divergence lies deeper in the snow than the maximum in the isotropic scattering case when the distance is measured in photon mean free paths (i.e. in terms of optical depth τ = χz). This behavior is shown in Figure 5. Again, we see that for small solar altitudes, the Sun heats the top few centimeters of the snow with respect to the snow further down in the layer. However, the flux divergence is much flatter than was the flux divergence for isotropic scattering.

Fig. 5.

Plot of the flux divergence —d{ΦΦ₀}/dτ versus the optical depth τ for various solar altitudes, using the anisotropic scatterer with a phase function shown in Figure 2

We also observe that the non-exponential behavior of the flux extends deeper into the snow than it did in the case of isotropic scattering. Using Figure 5, we find that the flux divergence becomes exponential below τ≈ 5, Thereafter

where we measure v_a ^’ to be 13.1 from the g.raph. This in turn implies that the extinction coefficient is now x= 190 m_-1 so that τ = 1 corresponds to about 0.526 cm. Since τ = 1 corresponds roughly to the optical depth at which one can distinguish objects through a layer, the new value of the extinction coefficient g.ives a more accurate idea of the thickness a layer of snow must have to be “opaque” to an underlying surface. The estimate of the physical depth corresponding to τ = 1 was incorrectly g.iven in paper I as 0.851 cm, and should have been x⁻¹ = (0.851 cm⁻¹)⁻¹ = 1.175 cm. This would require a rather thick layer to blanket the g.round so that it could not be seen. The new value of χ appears to bring the theoretical predictions of dΦ/dz into better agreement with the observations of Reference Ambach and MockerAmbach and Mocker (1959). The volume melting rate g.iven in Equation (42) of the first paper is changed as a result of anisotropic scattering to 0.569 kg cm⁻¹ S ⁻¹. However, the melting rate per unit area is unchanged. We have also calculated the asymptotic ratio of upward to downward flux deep in the medium. We find φ↑_as|φ↓_as = 0.815 for the phase function g.iven by Table I.

Fig. 6.

Plot of the angular dependence of the intensity of light I transmitted through slabs of various optical depths τ for neat normal incidence (μ ₀ = 0.995) on the top of the slab.

Transmission of Light through a Finite slab of snow

The transmission of light through a finite layer of snow is of considerable importance because it is one of the ways in which the “extinction coefficient” of snow is determined. Figure 6 shows the intensity transmitted by a thin layer of snow under nearly normal incidence for the phase function g.iven by Table I. It is clear that thin layers of snow (0.025 ^m thick or less) may g.ive an entirely erroneous impression of the character of the flux extinction. Such an effect was reported by Reference Giddings and LaChapelleGiddings and LaChapelle (1961), and indeed, was predicted by them on the basis of a diffusion theory for radiation. The steep portion of the curve of extinction coefficient versus optical depth is caused by the extinction of light transmitted directly through the slab without being scattered. It is clear that as the slab becomes thicker, the relative contrast between the directly transmitted light and the diffusely transmitted light decreases, until by τ = 7, only the diffuse light is coming through the slab. The strong directionality of the transmitted light for slabs up to about 3 cm thickness means that one may have to make special care in measuring the flux extinction, since detectors are often strongly direction dependent (Reference GillhamGillham, 1970).

Conclusions

We have extended the computations done in paper I to include anisotropic scattering by finite layers of snow. In doing so, the single scattering albedo and the phase function were adjusted until the calculated angular reflectance of a deep snow layer g.ives reasonable agreement with the observed reflectances. It has been shown that the flux divergence may be expected to peak slightly further beneath the surface than was the case for isotropic scattering (as measured in photon mean free paths). Finally, it has been suggested that measurements of the flux extinction coefficient using the light transmitted through finite layers of snow (under 2.5 cm thick) may lead to serious errors if extrapolated to effectively semi-infinite layers below the boundary layer in the latter.

Acknowledgements

The authors wish to thank Dr Lewis House and Dr Thomas Schlatter of NCAR for a critical reading of this paper. One of us (B.R.B.) would also like to g.ive especial thanks to Dr Thomas Schlatter of NCAR, Dr Kevin Pang of LASP, and Dr Roger Berry of INSTAR for their help in suggesting useful and pertinent references.

MS. received 20 August 1973 and in revised form 2 July 1974