In recent studies (Hu et al., 2022; Lu et al., 2022), a simple relationship, < L > = 2 H , between snow depth (H) and the average traveling distance between a photon’s entry and exit through snow (<L>) was discovered through Monte Carlo simulations. This simple relationship is valid for snow with various densities and scattering phase functions. It has been shown that the directly retrieved snow depths achieved by applying this simple relationship to the multiple scattering signals of ICESat-2 lidar measurements agree reasonably well with aircraft-based snow depth measurements (Lu et al., 2022).

Analyzing the snow bidirectional reflectance of a simple two-stream radiative transfer solution, the objective of this paper is to derive a simple, analytical relationship between snow depth and the average distance photons travel inside snow, to derive the expression for the snow extinction coefficient as an analytical function of < L > and < L 2 > , and to derive the solution for snow grain size as an analytical function of lidar-measured snow surface reflectance at the 180° backscatter direction. Then, snow density can be obtained analytically from the snow extinction coefficient and snow grain size. Thus, a multi-wavelength lidar is capable of measuring not only snow depth but also the snow extinction coefficient, snow grain size, and snow density. We also present an alternative method to measure snow density directly from lidars using the vibrational Raman scattering signals from snow.

When a space-based lidar receiver’s ground footprint size is a few meters larger than the laser spot size, the reflectance of the laser beam by a layer of snow is equivalent to the bidirectional reflectance of sunlight due to reciprocity. This is because snow extinction coefficients are greater than 100 1/m and thus the lidar receiver collects all the multiple scattering of the laser light. Light propagation in the snow can be described by the following one-dimensional radiative transfer equation, μ d I τ , μ d τ = I τ , μ − ω 2 ∫ − 1 1 P μ , μ ′ I τ , μ ′ d μ ′ − ω 4 π P μ , μ 0 F 0 ⁡ exp ⁡ ⁡ − τ μ 0 (1) where I is radiance, µ the cosine of viewing zenith angle, τ the optical depth, ω the single-scattering albedo, P μ , μ ′ the scattering phase function, and F ₀ the solar irradiance (normal to the beam).

The diffuse component of the radiance can be computed with a simple two-stream solution and a modified Eddington approximation that discretizes the differential-integral equation as follows: d I + d τ = r 1 I + − r 2 I − − r 3 ω π μ 0 F 0 ⁡ exp − τ μ 0 (2) d I − d τ = − r 1 I − + r 2 I + + 1 − r 3 ω π μ 0 F 0 ⁡ exp ⁡ ⁡ − τ μ 0 (3)

In the above equations, I + and I − are the upward and downward diffuse radiances with cosine of viewing zenith angles equal to −1 and 1, r 1 = 1 4 7 − ω 4 + 3 g , r 2 = − 1 4 1 − ω 4 − 3 g , r 3 = 1 4 1 − 3 g μ 0 , and g is the asymmetry factor of the scattering phase function P cos ⁡ θ , defined as g = 1 2 ∫ − 1 1 P x x d x . The parameter x = cos ⁡ θ , θ is the scattering angle.

This modified Eddington approximation is optimized for calculating nadir-viewing azimuth independent diffuse radiance measurements of snow in order to properly capture its dependence on absorption, which is view angle and solar zenith angle dependent. Thus, for this modified Eddington approximation of diffuse radiance calculation, the source terms is different from the standard Eddington approximation for hemispheric radiative flux calculations, and r 3 = 1 4 P ( μ 0 , μ ′ = − 1 ) = 1 4 1 − 3 g μ 0 , with P μ 0 , − μ ′ being the first two terms of the Legendre expansion of the Henyey-Greenstein scattering phase function, P μ 0 , − μ ′ = ∑ l = 0 1 2 l + 1 g l P l − μ ′ P l μ 0 . This is different from the Eddington approximation, where the solution is for hemispheric diffuse flux measurements, thus, r 3 = ∫ 0 1 P μ 0 , − μ ′ d μ ′ = 1 2 ∫ 0 1 1 − 3 g μ ′ μ 0 d μ ′ = 1 4 2 − 3 g μ 0 . For a snow-like media with weak absorption, the modified Eddington approximation agrees with Monte Carlo simulations of space-based lidar measurements (red line in Figure 1), while absorptions using the standard Eddington approximation and quadrature method (Meador and Weaver, 1980) are significantly less than those from Monte Carlo simulations (blue and black lines).

The eigenvalues of the above two equations are k = r 1 2 − r 2 2 = 1 4 7 − ω 4 + 3 g 2 − 1 − ω 4 − 3 g 2 = 1 4 48 1 − ω − g ω + g ω 2 = 3 1 − ω 1 − g ω (4)

The bidirectional reflectance at the top of the layer (Meador and Weaver, 1980) is R μ 0 , μ = ω 1 − k μ 0 a + k r 3 e k τ − 1 + k μ 0 a − k r 3 e − k τ − 2 k r 3 − a μ 0 e − τ / μ 0 1 − k 2 μ 0 2 k + r 1 e k τ + k − r 1 e − k τ . (5)

Here, a = r 1 + ( r 2 − r 1 ) 1 − r 3 and τ is the total optical depth of the snow layer.

For nadir-pointing lidar ( μ 0 = 1 , μ = 1 ) measurements of conservative scattering ( ω = 1 ) in an optically thick medium with the lidar receiver’s footprint covering nearly all multiple scattering of laser light, the vertically integrated attenuated lidar backscatter is equivalent to the bidirectional reflectance R μ 0 = 1 , μ = 1 , ω = 1 , with r 1 = r 2 = a = 3 4 1 − g , k = 0 , and τ ≫ 1 , e − τ / μ 0 ≈ 0 , R μ 0 = 1 , μ = 1 , ω = 1 = r 1 τ + r 3 − r 1 μ 0 1 + r 1 τ = 3 1 − g τ − 2 3 1 − g τ + 4 . (6)

This equation is similar but not identical to Eq. (7.87) in Stamnes et al. (2017) for g = 0 (isotropic scattering).

For a snow measurement lidar with its receiver footprint diameter a few meters greater than that of the laser spot, the lidar receiver can capture nearly all multiple scattering signals in the 180-degree backscatter direction. Increasing the footprint size will not change the lidar measurements. Thus, the lidar measurements of snow can be considered as 1) a laser beam with a divergence angle equivalent to the field-of-view angle of a passive sensor pointing at nadir; 2) a receiver as far as the Sun at solar zenith angle = 0° with a footprint near infinity. Due to reciprocity, the bidirectional reflectance, R μ 0 = 1 , μ = 1 , can also be expressed as a function of the vertical-integrated attenuated backscattering profile, I L , of the lidar measurement, R μ 0 = 1 , μ = 1 = ∫ 0 ∞ I L d L . (7)

Here, L is the distance of photons traveled within the medium. L = c t e x i t − t e n t r y , of which c is the speed of light and t e x i t and t e n t r y are the time the photons exit and enter the snow, respectively.

In previous studies (Hu et al., 2022; Lu et al., 2022), a simple relationship between the snow depth, H, and the averaged distance of photons traveling inside the non-absorbing medium, <L>, is H = 1 2 < L > = 1 2 ∫ 0 ∞ I 0 L , ω = 1 L d L ∫ 0 ∞ I 0 L , ω = 1 d L . (8)

Here, I( L , ω = 1 ) is the lidar backscatter profile measurements of optically thick conservative media, I 0 L , ω = 1 = I L , ω ≠ 1 e σ a L (9) where I ( L , ω ≠ 1 ) is the backscatter profile of snow, σ a is the absorption coefficient of snow, and c is speed of light. The bidirectional reflectance of an extremely weakly absorbing medium can be expressed as a function of <L> and the bidirectional reflectance of the non-absorbing medium, R μ 0 = 1 , μ = 1 , 1 − ω → 0 = ∫ 0 ∞ I L , ω d L = ∫ 0 ∞ I 0 L , ω = 1 e − σ a b s L d L ≈ ∫ 0 ∞ I 0 L , ω = 1 1 − σ a b s L + o σ a b s 2 L 2 d L ≈ ∫ 0 ∞ I 0 L , ω = 1 d L ∗ 1 − σ a b s ∫ 0 ∞ I 0 L , ω = 1 L d L ∫ 0 ∞ I 0 L , ω = 1 d L = ∫ 0 ∞ I 0 L , ω = 1 d L ∗ 1 − σ a b s < L > (10) R μ 0 = 1 , μ = 1 , 1 − ω → 0 ≈ R μ 0 = 1 , μ = 1 , ω = 1 ∗ 1 − σ a b s < L > (11) < L > = 1 − 1 σ a b s R μ 0 = 1 , μ = 1 , 1 − ω → 0 R μ 0 = 1 , μ = 1 , ω = 1 . (12)

Similarly, we can also use the higher order Taylor expansion, e − σ a b s L ≈ 1 − σ a b s L + 1 2 σ a b s 2 L 2 , to derive R μ 0 = 1 , μ = 1 , 1 − ω → 0 ≈ R μ 0 = 1 , μ = 1 , ω = 1 ∗ 1 − σ a b s < L > + 1 2 σ a b s 2 < L 2 > (13) < L 2 > = 2 σ a b s 2 R μ 0 = 1 , μ = 1 , 1 − ω → 0 R μ 0 = 1 , μ = 1 , ω = 1 − 1 − σ a b s < L > . (14)

Here, σ a b s is the absorption coefficient. From the two-stream solution, it is found here that R μ 0 = 1 , 1 − ω → 0 ≈ R μ 0 = 1 , ω = 1 ∗ 1 − 2 σ a b s H , and thus proves that <L>=2H. The following is the derivation:

For the limit of an extremely weakly absorbing optically thick medium where k τ < 1 , 1 − ω ≪ 1 , e ± k τ ≈ 1 ± k τ + k τ 2 2 ± k τ 3 6 , and τ ≫ 1 , e − τ μ 0 ≈ 0 .

Using Eq. 5 and the above relation, we can rewrite R as R μ 0 = 1 , μ = 1 , 1 − ω → 0 ≈ ω 1 − k μ 0 a + k r 3 1 + k τ + k τ 2 2 + k τ 3 6 − 1 + k μ 0 a − k r 3 1 − k τ + k τ 2 2 − k τ 3 6 1 − k 2 μ 0 2 k + r 1 e k τ + k − r 1 e − k τ ≈ a τ 1 + k 2 τ 2 6 + r 3 − a 1 + k 2 τ 2 2 1 + k 2 τ 2 2 + r 1 τ 1 + k 2 τ 2 6 = − 2 + 3 1 − g τ 1 + k 2 τ 2 6 − 4 k 2 τ 2 6 4 + 3 1 − g τ 1 + k 2 τ 2 6 + 8 k 2 τ 2 6 ≈ − 2 + 3 1 − g τ 1 + k 2 τ 2 6 − 4 k 2 τ 2 6 ∗ 3 1 − g τ 4 + 3 1 − g τ 1 + k 2 τ 2 6 + 8 k 2 τ 2 6 ∗ 3 1 − g τ ≈ − 2 + 3 1 − g τ 4 + 3 1 − g τ 1 − 2 k 2 τ 2 3 1 − g τ + O k 4 τ 4 = − 2 + 3 1 − g τ 4 + 3 1 − g τ 1 − 2 1 − ω τ = R μ 0 = 1 , μ = 1 , ω = 1 1 − 2 σ a b s H . (15)

Thus, when the absorption optical depth 1 − ω τ approaches zero, the reduction of snow bidirectional reflectance R due to absorption is proportional to twice the absorption optical depth of the snow layer. 2 H = 1 − 1 σ a b s R μ 0 = 1 , μ = 1 , ω → 1 R μ 0 = 1 , μ = 1 , ω = 1 (16)

From Eq. 12, < L > = 1 − 1 σ a b s R μ 0 = 1 , μ = 1 , ω → 1 R μ 0 = 1 μ = 1 , , ω = 1 , thus, < L > = 2 H . (17)

A Appendix A provides a different derivation of Eq. 17. It is important to note that the above equation is valid for weakly absorbing or non-absorbing media, of which σ a b s τ < 1 . Absorption reduces the backscatter signals, I L , ω < 1 = I L , ω = 1 e − σ a b s L . To apply it to snow measurements, an absorption correction is needed when <L> is computed, < L > = ∫ 0 ∞ I L , ω = 1 L d L ∫ 0 ∞ I L , ω = 1 d L = ∫ 0 ∞ I L , ω < 1 e σ a b s L L d L ∫ 0 ∞ I L , ω < 1 e σ a b s L d L . (18)

At the top of the snow, the Raman scattering of snow is proportional to snow density. Thus, lidar can provide the most direct measurements of snow density by adding vibrational Raman channels (around 3,300 cm⁻¹ away from the laser frequency). The concept of Raman scattering for snow density measurement is similar to the ice water content measurements estimates from the Raman scattering of ice clouds (e.g., Wang et al., 2004). Assuming the Raman backscatter cross section of each ice molecule in snow as C R a m a n , the peak Raman backscatter of snow near the top of the snow from a lidar measurement is C R a m a n N m . Here, N m is the number of molecules the laser light will interact with in a unit volume. N m = s ρ s n o w 18 ∗ 6.022 ∗ 10 23 , where ρ s n o w is the snow density and s is a scaling factor that is related to the increased path of lights reflecting inside snow particles. Assuming snow particles are spherical, the scaling factor can be estimated using a ray tracing technique (Bohren and Barkstrom, 1974), s ≈ 0.84 ∗   π d 3 4 / π d 3 6 ≈ 1.26 . Thus, the peak Raman backscatter cross coefficient at the top of the snow, β R a m a n , is a function of snow density, β R a m a n ≈ 1.26   ∗   C R a m a n   ∗   ρ s n o w 18   ∗   6.022   ∗   10 23 (19) where C R a m a n ≈ 4.2 ∗ 10 − 33 m²sr⁻¹ for UV laser at 355 nm (Pershin et al., 2014; Reichardt et al., 2022).

Another benefit of the Raman backscatter measurements is the effective discrimination between snow and trees. Raman signals of trees are significantly weaker compared with snow. Depending on the tree types, Raman shifts of trees in general are between 100 cm⁻¹ and 2,000 cm⁻¹ (e.g., Sevetlidis and Pavlidis, 2019), depending on the tree types, which differ significantly from that of snow (3,200 cm⁻¹ to 3,500 cm⁻¹).

When the absorption coefficient of snow approaches zero, the reduction of snow bidirectional reflectance due to absorption equals the absorption coefficient multiplied by the average path-length, i.e., the distance traveled by photons during temporal intervals between their entry to and exit from the snow measured by lidar, <L>. Using a simple two-stream radiative transfer equation, this study demonstrates that, when the absorption coefficient approaches zero, the reduction of snow bidirectional reflectance due to absorption is proportional to twice the absorption optical depth of the snow layer, which equals the absorption coefficient multiplied by twice the snow depth H. This solution holds for all scattering phase functions and extinction coefficients. Thus, the two-stream solution suggests that <L> =2H, confirming the finding of Hu et al. (2022) based on Monte Carlo simulations.

Lidar also measures the extinction coefficient of snow accurately using the second moment of the distribution of the photon path-lengths within the snow. The relationship between extinction coefficient, 1 − g σ e x t = 8 < L 2 > < L > 3 , is also verified with the analytic radiative transfer solution.

Using the two-stream radiative transfer theory, we also demonstrate that snow grain size can be derived from the lidar measurements of snow bidirectional reflectance, which is equivalent to the integrated lidar backscattering of the snow layer. Snow density, which is a function of the snow extinction coefficient and snow grain size, can be measured by the multi-wavelength lidar measurements of snow. The snow grains are simplified as equivalent spheres in the two-stream radiative transfer studies. As the microstructure of snow is complicated (e.g., Ding and Tsang, 2010; Xiong and Shi, 2013), our simplification may introduce uncertainties in the snow density estimates. Alternatively, snow density can be measured effectively by adding vibrational Raman scattering channels to the lidar measurements. Raman scattering measurements also enable the discrimination between snow and trees.

In order to obtain the above analytical solutions, it is assumed in the derivations that snowpack consists of snow particles of uniform diameter. For physical reasons (e.g., Dawson et al., 2017), snow density should increase with depth in general. Future efforts are needed to address two questions. First, does the derived snow density represent the average density of snowpack? Second, could we use the same functional form (e.g., Eq. 35) but with different coefficients to represent the average density of snowpack or the density of the top snow layer? For the former, snow water equivalent (SWE) can be directly obtained from snow depth multiplied by the average snow density. For the latter, we recognize that the snow density of the top layer is most uncertain, as the density of lower layers can be estimated based on physical processes in the snowpack. For instance, if the top snow layer density can be retrieved, we could use the snow density model in Dawson et al. (2017) driven by daily temperature and precipitation to obtain the snow density for different snow layers (up to 10 layers). Then the average snow density can be obtained to compute SWE. In this way, a multi-wavelength lidar would be able to retrieve snow depth, snow density, and SWE, providing an innovative approach for one of the seven observables for the Earth System Explorers satellite mission competition, as recommended by the 2017–2027 Decadal Survey for Earth Science and Applications from Space (National Academies, 2018).