Theoretical study of solar light reflectance from vertical snow surfaces

The influence of horizontal and vertical inhomogeneity of snow surfaces on solar light reflectance is studied using the radiative transfer theory (RTT). We compared 1-D RTT and 2-D RTT and found that large errors are produced if the 1-D RTT is used for the calculation of the snow reflection function (and, therefore, also in the retrievals of the snow grain radii) in 2-D measurement geometries. Such 2-D geometries are common in the procedures for the determination of the effective snow grain radii using near-infrared photography and spectroscopy of vertical snow walls. In particular, we have considered three cases for the numerical calculations: (1) the case with no black film; (2) the case with a black film at the pit’s bottom; (3) the case with a black film at the pit’s bottom and also at one of the vertical snow walls.


Introduction
Optical measurements are commonly used to derive snow microphysical parameters from plane-parallel snow layers (Kokhanovsky et al., 2011).In particular, snow grain size is obtained from near-infrared (NIR) measurements (in the spectral range 865-1240 nm) of intensity of solar light reflected from flat snow layers.The corresponding retrieval algorithms are based upon the physical phenomenon of the enhancement of light absorption by larger ice grains (and as a consequence, a smaller light reflectance for snow layers with larger grains).The main problem with such a method is that only upper snow layers can be observed.The information on the snow microphysical parameters and snow pollution in deeper layers cannot be retrieved because of high absorption of NIR radiation by snow grains.As a matter of fact, NIR radiation does not penetrate deep into a snowpack and, therefore, does not contain information on the properties of snow from the depths above 1-5cm depending on the size of particles and the wavelength (Kokhanovsky and Rozanov, 2012).To avoid this problem, recently measurements along vertical snow walls have become popular (see, e.g.Fig. 1 in Matzl and Schneebeli, 2006;and Fig. 2 in Painter et al., 2007).Also measurements along the length of cylindrical holes in snow are used (Barker and Korolev, 2010;Arnaud et al., 2011).
In most of cases (see e.g.Kokhanovsky et al., 2011) the 1-D transfer theory valid for plane-parallel slabs is used for the interpretation of optical measurements and determination of snow grain sizes.Although there could be some influences of 3-D effects (e.g.shadowing from the snow walls, enhancement of brightness, etc.) on corresponding measurements.For the measurements involving 2-D and 3-D geometries (e.g.along snow walls), the approach based on the correlation of the reflectance and the snow grain size or the snow specific surface area (see e.g.Matzl and Schneebeli, 2006) is used.This is because more quantitative approaches based on the solution of radiative transfer equation in 2-D and 3-D geometries have not been developed in applications relevant to optics of vertical snow walls.
The aim of this work is twofold.Firstly, we develop software, which can be used for studies of 3-D effects in snow and, secondly, we study the corresponding effects using numerical simulations.Only optical snow parameters (extinction coefficient, single scattering albedo, and phase function) are used in the analysis.This makes it possible to perform a study with greater generality.The link to the snow grain size and concentration of pollutants is obvious (Kokhanovsky et al., 2011).The paper is structured as follows.In the next section we introduce the radiative transfer equation and boundary conditions relevant to the studies of light propagation in snow.The numerical algorithm developed for the solution of the corresponding integro-differential radiative transfer equation in the 2-D geometry is described in Sect.3. The results of numerical experiments are reported in Sect. 4. The present study can be used to design and interpret real-world experiments relying on the spectrometry of vertical snow walls.

Radiative transfer equation and boundary conditions
It is assumed that the surface of the snow is flat (no sastrugi, no microstructures on the snow surface).Let us assume that there is a pit with the width D, the length L and the depth H in the snowpack with the width 2X and the height H , see Fig. 1.The pit is covered by a sheer film to convert the direct solar light to diffuse one.Reflected radiation is registered along the line AB on the vertical wall of the pit and along the line CS at the top of the pit.
To find radiation intensity in this region, we introduce the coordinate system with the origin at the centre point O at the pit's bottom, see Fig. 1.We assume that the length L (10-15 m) is larger than the typical width D (about 1-3 m).Therefore, the radiation intensity near the central plane y = 0 of the pit depends only on the spatial coordinates x and z.The dependence on the third coordinate can be neglected.So the problem can be considered in the 2-D framework in the plane y =0.This reduces calculations as compared to the 3-D modelling.Moreover, the region is symmetrical with respect to the plane x = 0, see Fig. 1, hence only the half-region The radiative transfer equation (RTE) for the monochromatic radiation intensity I takes the following form in the case under consideration: The function I (x, z, θ, ϕ) depends on the spatial coordinates x, z and angles θ, φ, defining direction of the radiation transfer, see Fig. 2. The first term in Eq. ( 1) gives the change of intensity in the direction and the second term represents the extinction of radiation by the medium.The integral describes the re-radiation of scattered light, here incident light has the direction (θ , ϕ ) and re-radiated light has the direction (θ, ϕ).
The pit [0, R] × [0, H ], R = D/2, is filled by air, whereas the medium out of the pit is snow.Then it follows that for the extinction coefficient: for the single scattering albedo: for the scattering phase function: Note, the single scattering albedo ω snow 0 (z) is considered as a piecewise function of depth, which describes a layered snowpack; the special case ω snow 0 (z) = const corresponds to a homogeneous snow layer.
One needs to define the boundary conditions for Eq. ( 1).The intensity of the radiation entering the region is defined on each boundary at the corresponding angular intervals, see    ( , , , , , ) with albedo A(x).Therefore, it follows that: Such a surface reflects incident radiation uniformly into all possible directions.If the bottom is covered by a black film absorbing all incident radiation, the albedo A(x) is made to be equal to zero and the bottom boundary is called "black".
It is assumed that the radiation does not enter the region via the right boundary x = X: Reflecting condition with albedo A s is defined on the left boundary x = 0: = A s I (0, z, θ, π − ϕ).
When A s = 1 Eq.( 7) gives the condition of symmetry of the whole region [−X, X] × [0, H ] with respect to the plane x = 0. Actually, the solution of the RTE (Eq. 1) in the whole region [−X, X] × [0, H ] is symmetrical: I (x, z, θ, ϕ) = I (−x, z, θ, π −ϕ).Therefore, one can consider the half-region [0, X] × [0, H ] under the boundary condition (Eq.7).When A s = 0, the boundary x = 0 is black.It means that this boundary is covered by a black film absorbing all incident radiation.
The diffuse source on the top boundary is imposed.Therefore, it follows that: where S 0 is the incident light irradiance.
We will study three problems depending on albedo A s at the left snow wall and the albedo A(x) at the bottom.First, under A s = 1, A(x) = A snow > 0 one has the pit with no black film.Second, under A s = 1, A(x) = 0 as x ≤ R, A snow else, the black film lies only at the pit's bottom.Third, under A s = 0 and A(x) = 0 as x ≤ R, A snow else, the black film covers both the bottom and the vertical plane x = 0.
Relative intensities of reflected radiation at the snow surface in the directions * and * * : The function Î (x) can be used to retrieve the optical properties of the upper layer of the snow (up to 5cm in depth).This function can be approximated by the piece-constant function: where the values I snow (H, * ) and I air (H, * ) are obtained via the two 1-D radiative transfer models as shown here (Chandrasekhar, 1950): (11) The boundary conditions are ( see the definition of angles in Fig. 3) Equation ( 10) is used to find the solution over areas, where snow is present.Otherwise, Eq. ( 11   ( , , , , , ) Eq. ( 10) is solved along the line CO passing through the centre of the pit, whereas the problem formulated in Eq. ( 11) is defined along the line TV passing through the centre of the snowpack, see Fig. 1.
The function Ĩ (z) is found by measuring light reflectance from snow walls and can be used to retrieve the optical and microphysical properties of the snow layers at any depth.This is not possible if, for example, the measurements of the light reflectance from the snow top (Kokhanovsky et al., 2011) are analysed.This is due to weak dependence of the snow reflectance in the UV and, also in the visible, on the size of particles and small penetration depths of IR radiation (sensitive to the snow microstructure) into snowpack.The retrieval algorithms are based upon the assumption that the registered radiation intensity is a constant function of the spatial coordinates in each homogeneous sub-region (layer) and this constant value does not depend on neighboring sub-regions (layers).It is actually believed that each sub-region (layer) of the snowpack can be considered separately from others.The constant value Ĩ (z) in each homogeneous layer is often believed equal to the value I snow (H, * ) for the homogeneous snowpack under the same optical properties.Such an approach is "the horizontal 1-D transfer model".The difference of the 2-D solution and the 1D solutions is termed as "2-D effects".We will check the accuracy of the 1D radiative transfer models using exact solutions of the 2-D problem (see Eqs. 1-8).We do not use the term "3-D effects" because the 2-D problem is under consideration.

Numerical algorithm
Below follows an outline of the numerical method used by us for the solution of the above radiative transfer problem.We introduce a quadrature with nodes m {θ m , φ m }, see Fig. 4a,  and weights m , m = 1, . . .M. For this purpose we use the mesh over angle θ : and the mesh over angle ϕ for each interval θ −1/2 , θ +1 / /2 of the mesh over θ:  Points θ +1/2 and ϕ n+1/2, have been chosen so that squares of all cells ϕ n−1/2, , ϕ n+1/2, × θ −1/2 , θ +1/2 are identical.We define the node {ϕ n , θ } in each cell, see Fig. 4b, and renumber all nodes with a single index m.Further we define the weight m as a square of the corresponding cell m .Then we approximate the continuous function I (x, z, θ, φ) with functions I m (x, z) = I (x, z, θ m , φ m ) and replace the scattering integral in Eq. ( 1) with the quadrature sum: where The coefficients ρ nm (x, z) correspond to the light scattering event from the direction n {θ n , φ n } to the direction m {θ m , φ m }.Therefore, they are the integrals of a complicated forward-peaked phase function ρ(x, z, θ m , ϕ m , θ , ϕ ) over the cell n under fixed values of the angles {θ m , ϕ m }.To find the integral of a forward-peaked phase function one introduces the additional quadrature in the cell n by the nodes j, n {θ j,n , ϕ j, n } and the weights j, n , j = 1, . . ., L n , here L n is the number of the additional nodes.This quadrature is refined in the subregions of the cell n , where the integrand function has a great gradient, see Fig. 5, where the additional nodes j, n {θ j, n , ϕ j, n } are designated by the black circles.The following equality is always kept: Then the coefficient ρ nm (x, z) can be found by a quadrature sum: So Eq. ( 1) for the functions I m (x, z) takes the form (with account for Eq.10): where the derivative ∂I ∂ m is written as: The values are projections of the unit vector m onto coordinates axes x and z, see Fig. 4a.To solve the system of differential equations (Eq.16) for the functions I m (x, z), we introduce a regular mesh over spatial variables x, z: is considered as a homogeneous one.Integrating Eq. ( 16) over this cell, one obtains the exact algebraic relation: Here the values σ k,j , ω 0, k, j , ρ nm, k, j correspond to the cell The values I m, k, j , I m, k±1/2, j , I m, k, j ±1/2 are averaged light intensities defined by the following integrals: The boundary conditions for Eq. ( 19) follow from Eqs. (5-8): To close Eqs.( 19) and (23-25) one needs additional relations.They are taken as where the function s(ξ ) = sgn(ξ ) = 1 as ξ > 0, −1 as ξ < 0 and the parameters v m, k, j , u m, k, j are defined on the interval [0, 1].Therefore, the piece-linear approximation to the solution is sought in the spatial cell, see Fig. 6.Here the parameters v m, k, j , u m, k, j define the variation of the solution in the cell (Carlson, 1972).
Let the solution in the nodes n , n = 1, . . ., m-1 and n = m + 1, . . .M be known.Then the solution of the system (Eqs.21 and 25-29) for the node m can be found by the so-called sweep procedure in a following way.
If the angles θ m , ϕ m are from other intervals, then indices are to be sorted in yet another order.Generally, the index k increases as s(ξ m ) > 0 and decreases as s(ξ m ) < 0, the index j increases as s(γ m ) > 0 and decreases as s(γ m ) < 0.
To solve the system (Eqs.21 and 25-29) for all nodes m , the iterative Seidel's method (Saad, 2000) is used.In this method the already obtained values I n, k, j are used to calculate the right side of Eq. ( 21) and further the values I m, k, j at other nodes.x The previous version of the presented algorithm was outlined by Sokoletsky et al. (2009), where it was applied to the calculation of solar light reflectance by natural sea waters.There the scattering phase functions were defined by their values in nodes of a very refined mesh over the interval [−1, 1] and approximated by piecewise linear functions.
Here the scattering phase functions are given by their Legendre coefficients.Furthermore, we apply the adaptive method of choosing additional meshes j, m to calculation of integrals (Eq.15).

Results of numerical experiments
All computations were done by the code RADUGA-6 (Nikolaeva et al., 2005;Sokoletsky et al., 2009) on the hybrid cluster k100 (http://www.kiam.ru/MVS/resourses/k100.html) assuming the following parameters:   6. snow albedo A snow = 0.8; 7. the diffuse source, when both a snowpack and a pit are covered by a sheer film.
We have selected a typical snow phase function as suggested by Kokhanovsky et al. (2011).The phase function does not depend strongly either on the wavelength (in the optical range) or on the size of ice grains.The extinction coefficient 1 mm −1 and the values of snow grain albedos in the range 0.98-1.0are typical for snow.
Both homogeneous and heterogeneous snowpack were under consideration.A homogeneous snowpack is defined by the constant single scattering albedo ω snow 0 .A heterogeneous snowpack contains a polluted layer, see Fig. 8.It was assumed that: Here parameter t is thickness of a polluted layer, ω0 is single scattering albedo of clean snow.
We define three experimental conditions (see Fig. 1): 1. no black film: 2. a black film is only on the pit's bottom EB: 3. a black film is on the pit's bottom EB and left boundary If the snow surface is covered by the black film, radiation is not reflected at this surface and does not influence the registered radiation intensity on the line AB.So the registered radiation intensity depends only on properties of the snow on the vertical wall AB.If there is no black film, one registers the radiation reflected by both the bottom and two walls of the pit.Then the registered radiation depends on optical properties of the snow at all walls of the pit.
The following parameters are used in the numerical calculations.
1. N = 800 is the number of the Legendre polynomials to represent both phase functions; 2. M = 360 is the number of nodes of the quadrature; one needs a dense quadrature to approximate the strongly anisotropic solution I (x, z, θ, φ); 3. K = 468, J = 1610 are numbers of cells of the spatial meshes.The mesh over z is refined in the vicinity of the top boundary z = H , where the intensity I (x, z, θ, φ) has a large gradient.The mesh over x is refined near snowpack wall AB, see Fig. 2, for the same reason.
Let us consider relative radiation intensity Ĩ (z) given by Eq. ( 9) on the vertical wall AB of a snowpack, see Fig. 2, in the direction * * , which is perpendicular to the wall AB and at the top boundary CS of the system in the zenith direction * .
The relative intensity at the horizontal line CS in the zenith direction * for homogeneous snowpack is given in Fig. 9.One can see that the intensity of reflected radiation has extrema near the air/snow boundary; similar effects are observed in clouds illuminated by direct solar light (Nikolaeva et al., 2005).In the problem under study a maximum of radiation intensity in the snow near the air/snow boundary is formed by radiation penetrating in the snowpack and only weakly absorbed near this boundary; the maximum enhances as snow absorption enhances.In a similar way, a minimum of radiation intensity arises outside of the snow near the air/snow boundary due to absorption of radiation by the snow.Thereby the extrema in the radiation intensity in Fig. 9 arise due to the neighbourhood of two different media (snow and air).Note the 1-D vertical transfer model leads to the  This deviation shows whether it is possible to consider the intensity Ĩ (z) as a constant function far from the upper and lower boundaries of the pit.In other words it shows whether the 1-D model is applicable to process measurement data along vertical walls of snowpack.
It follows from Fig. 11 that the 1-D transfer model is not applicable for experiments under the black film because in this case the deviation r(z) is less than 10 % only near the central point z = H /2. Actually, the size of the sub-region, where the deviation r(z) is less than 10 %, is equal to 7 cm -if both bottom and opposite sides are covered by the black film -and about 17 cm -if only the bottom is covered the black film.The size of this sub-region for the case without black film is about 55 cm.
The results for the pit without black film are shown in Figs.12-15.It should be stressed that in this case the radiation registered on the wall AB is reflected by the bottom and the opposite walls of the pit and depends on the optical properties of the whole surface of the pit.
Let us consider homogeneous snowpack.Here the deviation r(z) decreases as absorption decreases (see Fig. 12) and the width D of the pit increases (see Fig. 14).The function r(z) only weakly depends on the depth H (see Fig. 13).
At small values of the probability of photon absorption β = 1 − ω 0 and in broad pit, the deviation r(z) is less than the threshold value 10 % far from bottom and upper edges of the pit; here the 1-D model can be used.At the same time this deviation is large near the bottom and upper edges (boundary effects), where the 1-D model is not applicable.
The influence of heterogeneity of a snowpack on relative radiation intensity is presented in Figs.12-15.The thin polluted layer in the centre of the pure snowpack, see Fig. 8, leads to a minimum in reflected radiation intensity in the vicinity of the layer (the shadow of the minimum is spread over the whole wall, if absorption in snow is weak enough).Let us define the width of the spread of the optical influence of the polluted layer as the size of the sub-region, where the relative intensity of the polluted layer differs more than by threshold value b% from the relative intensity of the homo-   geneous snowpack: Here the point z = H /2 is the central point of the whole snowpack and the polluted layer and the function p(z) is defined by the relation: The function Ĩ0 (z) is the relative intensity in the homogeneous snowpack, Ĩ (z) is the relative intensity in the snowpack with the polluted layer.The widths of the spread of polluted layer optical influences t * are presented in Table 1 for the different single scattering albedo ω snow 0 outside of the polluted layer, the width of the polluted layer t and threshold value b.One can see that t * is always larger than the real width of the polluted layer t, especially when the outer snow is clean.The error in width of the polluted layer (when defined via the value t * ) can reach 200-400 % (see Table 1), especially if the polluted layer is thin.At the same time the value of the minimum of the relative intensity in the polluted   layer depends on the width of this layer and the single scattering albedo outside of this layer, see Fig. 15 and Table 1.Note that the minimum decreases as the width of the layer decreases and the albedo of surrounding medium increases.

Conclusions
We have presented the 2-D radiative transfer problem related to the reflection of solar light by a rectangular wide pit in a thick snow layer.Simulation (by the parallel code RADUGA-6) is based upon the mesh technique of the discrete ordinate method when peaked scattering phase functions of snow are exactly taken into account.A diffuse radiation source, produced by a sheer film covering a snowpack, is assumed.Such source models are close to those for real ground measurements.
We have checked whether the 1-D model, when the reflected radiation intensity is considered as constant function of the spatial coordinate in each homogeneous subregion of a snowpack, is applicable to describe the real measurements.We found that the 2-D effects (brightening and shadowing) on the top boundary of a snowpack near the vertical wall of the pit are significant in spite of a diffuse radiation source.
The 2-D effects are significant on the vertical wall of the pit in a homogeneous snowpack, especially near the upper boundary.At the same time, 2-D effects are less evident at large values of the pit's width far from its bottom and top boundaries, when snow is almost clean.
Additional 2-D effects arise in layered snowpack.Although minimum in intensity on a vertical wall of a pit is localized near a polluted layer, intensity out of minimum can be influenced by this polluted layer.
One can conclude that 1-D models can lead to large errors in the simulation of the measured radiation intensity on vertical walls of snow pits.The retrieval algorithms should, therefore, be based upon the 2-D and 3-D radiative transfer models.
, H, * /S 0 (9) are of interest.Here the function Î (x) defines the radiation intensity exiting from the top boundary in the zenith direction * .The function Ĩ (z) corresponds to radiation intensity reflected by the vertical wall AB of the snowpack (see Fig. 2) in the direction * * , perpendicular to the wall AB.

Fig. 6 .
Fig. 6.The piece-linear approximation to the RTE solution over the spatial variable x under ξ m > 0.
1. the region height H = 0.5 m, 0.6 m, 0.7 m, the region semi-width X = 5 m, seeFig.air (aerosol)  scattering phase function ρ air is obtained via Mie theory, the snow phase function ρ snow (see Eq. 4) is found by geometrical optics theory as described byKokhanovsky et al. (2011), see Fig.7.

Fig. 7 .
Fig. 7.The scattering phase functions.The molecular scattering is ignored and the air phase function is assumed to be equal to that of atmospheric aerosol.

Fig. 8 .
Fig. 8.The 3-D geometry of the region with the central polluted layer.

Fig. 9 .
Fig.9 air snow Fig. 9. Relative intensity Î (x) in the zenith direction * at the top boundary CS of the homogeneous snowpack.Width D = 1 m, depth H = 0.7 m and no black film, for the different single scattering albedos ω snow 0 .

Fig. 11 .
Fig. 11.Relative intensity Ĩ (z) in the direction * * (a) and the deviation r(z) (b) at the vertical wall AB of the homogeneous snowpack.The snow single scattering albedo ω snow 0 = 0.98, the width D = 1 m, and depth H = 0.7 m , with and without black film. Fig.12

Fig. 12 .
Fig. 12. Relative intensity Ĩ (z) in the direction * * (a) and the deviation r(z) (b) on the vertical wall AB of the heterogeneous snowpack.Single scattering albedo ω0 = 1 out of the inserted layer, width D = 1 m, depth H = 0.7 m, and no black film, for the different values of the thickness t of the inserted polluted layer.30

Fig. 13 .
Fig. 13.Relative intensity Ĩ (z) in the direction * * (a) and the deviation r(z) (b) on the vertical wall AB of the homogeneous snowpack.The snow single scattering albedo ω snow 0 = 0.98, the width D = 1 m, and no black film, for different depths H . Fig.13

Fig. 14 .
Fig. 14.Relative intensity Ĩ (z) in the direction * * (a) and the deviation r(z); (b) on the vertical wall AB of the homogeneous snowpack.The snow single scattering albedo ω snow 0 = 0.98, the depth H = 0.7 m for, no black film, different widths D. Fig.15

Fig. 15 .
Fig. 15.Relative intensity Ĩ (z) in the direction * * on the vertical wall AB of the heterogeneous snowpack.Width D = 1 m, depth H = 0.7 m, thickness of the inserted polluted layer t = 5 cm, and no black film, for different values of the single scattering albedo ωsnow 0

Table 1 .
The value of the minimum of the relative intensity Ĩmin and the width of the optical influence spread of the polluted layer t * (cm) for different values of the single scattering albedo ω snow 0 outside of the polluted layer at different widths of the polluted layer t.