Linking lidar multiple scattering profiles to snow depth and snow density: an analytical radiative transfer analysis and the implications for remote sensing of snow

Lidar multiple scattering measurements provide the probability distribution of the distance laser light travels inside snow. Based on an analytic two-stream radiative transfer solution, the present study demonstrates why/how these lidar measurements can be used to derive snow depth and snow density. In particular, for a laser wavelength with little snow absorption, an analytical radiative transfer solution is leveraged to prove that the physical snow depth is half of the average distance photons travel inside snow and that the relationship linking lidar measurements and the extinction coefficient of the snow is valid. Theoretical formulas that link lidar measurements to the extinction coefficient and the effective grain size of snow are provided. Snow density can also be derived from the multi-wavelength lidar measurements of the snow extinction coefficient and snow effective grain size. Alternatively, lidars can provide the most direct snow density measurements and the effective discrimination between snow and trees by adding vibrational Raman scattering channels.

1 Introduction

In recent studies (Hu et al., 2022; Lu et al., 2022), a simple relationship, $<L>=2H,$ 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.

2 $<\mathit{L}>=2\mathit{H}$: theoretical proof with two-stream radiative transfer solution

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,

$\mu \frac{dI\left(\tau ,\mu \right)}{d\tau }=I\left(\tau ,\mu \right)-\frac{\omega }{2}{\int }_{-1}^{1}P\left(\mu ,{\mu }^{\prime }\right)I\left(\tau ,{\mu }^{\prime }\right)d{\mu }^{\prime }-\frac{\omega }{4\pi }P\left(\mu ,{\mu }_0\right)F_0\exp \left(-\frac{\tau }{{\mu }_0}\right)(1)$

where I is radiance, µ the cosine of viewing zenith angle, τ the optical depth, ω the single-scattering albedo, $P\left(\mu ,{\mu }^{\prime }\right)$ 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:

$\frac{d{I}^{+}}{d\tau }=r_1{I}^{+}-r_2{I}^{-}-r_3\frac{\omega }{\pi }{\mu }_0{F}_0\exp \left(-\frac{\tau }{{\mu }_0}\right)(2)$

$\frac{d{I}^{-}}{d\tau }=-r_1{I}^{-}+r_2{I}^{+}+\left({1-r}_3\right)\frac{\omega }{\pi }{\mu }_0{F}_0\exp \left(-\frac{\tau }{{\mu }_0}\right)(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=\frac{1}{4}\left[7-\omega \left(4+3g\right)\right]$, $r_2=-\frac{1}{4}\left[1-\omega \left(4-3g\right)\right]$, $r_3=\frac{1}{4}\left[1-3g{\mu }_0\right]$, and g is the asymmetry factor of the scattering phase function $P\left(\mathit{cos}\theta \right)$, defined as $g=\frac{1}{2}{\int }_{-1}^{1}P\left(x\right)xdx.$ The parameter $x=\mathit{cos}\theta ,\theta$ 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=\frac{1}{4}P\left({\mu }_0,{\mu }^{\prime }=-1\right)=\frac{1}{4}\left(1-3g{\mu }_0\right),$ with $P\left({\mu }_0,-{\mu }^{\prime }\right)$ being the first two terms of the Legendre expansion of the Henyey-Greenstein scattering phase function, $P\left({\mu }_0,-{\mu }^{\prime }\right)={\sum }_{l=0}^1\left(2l+1\right)g^l{P}_l\left(-{\mu }^{\prime }\right)P_l\left({\mu }_0\right)$. This is different from the Eddington approximation, where the solution is for hemispheric diffuse flux measurements, thus, $r_3={\int }_0^{1}P\left({\mu }_0,-{\mu }^{\prime }\right)d\left({\mu }^{\prime }\right)=\frac{1}{2}{\int }_0^1\left(1-3g\mu {\prime \mu }_0\right)d{\mu }^{\prime }=\frac{1}{4}\left[2-3g{\mu }_0\right]$. 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).

FIGURE 1

FIGURE 1. Comparisons of absorptions in space-based lidar measurements of weak absorbing snow-like media computed from Monte Carlo simulations of lidar measurements (X-axis) and that from bidirectional reflectance calculated with the modified Eddington approximation (red line), and that from the hemispheric fluxes calculated from the standard Eddington approximation (blue line) and the quadrature method (black line).

The eigenvalues of the above two equations are

$k=\sqrt{{r_1}^2-{r_2}^2}=\frac{1}{4}\sqrt{{\left[7-\omega \left(4+3g\right)\right]}^2-{\left[1-\omega \left(4-3g\right)\right]}^2}=\frac{1}{4}\sqrt{48\left(1-\omega -g\omega +g{\omega }^2\right)}=\sqrt{3\left(1-\omega \right)\left(1-g\omega \right)}(4)$

The bidirectional reflectance at the top of the layer (Meador and Weaver, 1980) is

$R\left({\mu }_0,\mu \right)=\frac{\omega \left[\left(1-k{\mu }_0\right)\left(a+k{r}_3\right)e^{k\tau }-\left(1+k{\mu }_0\right)\left(a-k{r}_3\right)e^{-k\tau }-2k\left(r_3-a{\mu }_0\right)e^{-\tau /{\mu }_0}\right]}{\left(1-k^2{{\mu }_0}^2\right)\left[\left(k+r_1\right)e^{k\tau }+\left(k-r_1\right)e^{-k\tau }\right].}(5)$

Here, $a=r_1{+(r}_2-r_1)\left({1-r}_3\right)$ and $\tau$ is the total optical depth of the snow layer.

For nadir-pointing lidar (${\mu }_0=1,\mu =1$) measurements of conservative scattering ($\omega =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\left({\mu }_0=1,\mu =1,\omega =1\right),$ with $r_1=r_2=a=\frac{3}{4}\left(1-g\right),k=0,$ and $\tau \gg 1,e^{-\tau /{\mu }_0}\approx 0$,

$R\left({\mu }_0=1,\mu =1,\omega =1\right)=\frac{r_1\tau +r_3-{r_1\mu }_0}{1+r_1\tau }=\frac{3\left(1-g\right)\tau -2}{3\left(1-g\right)\tau +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\left({\mu }_0=1,\mu =1\right),$ can also be expressed as a function of the vertical-integrated attenuated backscattering profile, $I\left(L\right)$, of the lidar measurement,

$R\left({\mu }_0=1,\mu =1\right)={\int }_0^{\infty }I\left(L\right)dL.(7)$

Here, L is the distance of photons traveled within the medium. $L=c\left(t_{exit}-t_{entry}\right),$ of which c is the speed of light and $t_{exit}$ and $t_{entry}$ 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=\frac{1}{2}<L>=\frac{1}{2}\frac{{\int }_0^{\infty }{I}_0\left(L,\omega =1\right)LdL}{{\int }_0^{\infty }{I}_0\left(L,\omega =1\right)dL}.(8)$

Here, I($L,\omega =1$) is the lidar backscatter profile measurements of optically thick conservative media,

$I_0\left(L,\omega =1\right)=\mathrm{I}\left(\mathrm{L},\mathrm{\omega }\ne 1\right){\mathrm{e}}^{{\sigma }_{a}L}(9)$

where $I(L,\omega \ne 1$) is the backscatter profile of snow, ${\sigma }_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,

$\begin{array}{c}R\left({\mu }_0=1,\mu =1,\left(1-\omega \right)\to 0\right)={\int }_0^{\infty }I\left(L,\omega \right)dL={\int }_0^{\infty }{I}_0\left(L,\omega =1\right)e^{-{\sigma }_{abs}L}dL\approx {\int }_0^{\infty }{I}_0\left(L,\omega =1\right)\left[1-{\sigma }_{abs}L+o\left({{{\sigma }_{abs}}^{2}L}^2\right)\right]\hfill \\ dL\approx {\int }_0^{\infty }{I}_0\left(L,\omega =1\right)dL\ast \left[1-{\sigma }_{abs}\frac{{\int }_0^{\infty }{I}_0\left(L,\omega =1\right)LdL}{{\int }_0^{\infty }{I}_0\left(L,\omega =1\right)dL}\right]={\int }_0^{\infty }{I}_0\left(L,\omega =1\right)dL\ast \left[1-{\sigma }_{abs}<L>\right]\hfill \end{array}(10)$

$R\left({\mu }_0=1,\mu =1,1-\omega \to 0\right)\approx R\left({\mu }_0=1,\mu =1,\omega =1\right)\ast \left[1-{\sigma }_{abs}<L>\right](11)$

$<L>=1-\frac{1}{{\sigma }_{abs}}\frac{R\left({\mu }_0=1,\mu =1,1-\omega \to 0\right)}{R\left({\mu }_0=1,\mu =1,\omega =1\right)}.(12)$

Similarly, we can also use the higher order Taylor expansion, $e^{-{\sigma }_{abs}L}\approx 1-{\sigma }_{abs}L+\frac{1}{2}{{{\sigma }_{abs}}^{2}L}^2,$ to derive

$R\left({\mu }_0=1,\mu =1,\left(1-\omega \right)\to 0\right)\approx R\left({\mu }_0=1,\mu =1,\omega =1\right)\ast \left[1-{\sigma }_{abs}<L>+\frac{1}{2}{{{\sigma }_{abs}}^2<L}^2>\right](13)$

$<L^2>=\frac{2}{{{\sigma }_{abs}}^2}\left[\frac{R\left({\mu }_0=1,\mu =1,1-\omega \to 0\right)}{R\left({\mu }_0=1,\mu =1,\omega =1\right)}-\left(1-{\sigma }_{abs}<L>\right)\right].(14)$

Here, ${\sigma }_{abs}$ is the absorption coefficient. From the two-stream solution, it is found here that $R\left({\mu }_0=1,\left(1-\omega \right)\to 0\right)\approx R\left({\mu }_0=1,\omega =1\right)\ast \left[1-{2\sigma }_{abs}H\right]$, and thus proves that <L>=2H. The following is the derivation:

For the limit of an extremely weakly absorbing optically thick medium where $k\tau <1,1-\omega \ll 1,e^{\pm k\tau }\approx 1\pm k\tau +\frac{{\left(k\tau \right)}^2}{2}\pm \frac{{\left(k\tau \right)}^3}{6}$, and $\tau \gg 1,e^{-\frac{\tau }{{\mu }_0}}\approx 0.$

Using Eq. 5 and the above relation, we can rewrite R as

$\begin{array}{c}R\left({\mu }_0=1,\mu =1,\left(1-\omega \right)\to 0\right)\approx \frac{\omega \left[\left(1-k{\mu }_0\right)\left(a+k{r}_3\right)\left(1+k\tau +\frac{{\left(k\tau \right)}^2}{2}+\frac{{\left(k\tau \right)}^3}{6}\right)-\left(1+k{\mu }_0\right)\left(a-k{r}_3\right)\left(1-k\tau +\frac{{\left(k\tau \right)}^2}{2}-\frac{{\left(k\tau \right)}^3}{6}\right)\right]}{\left(1-k^2{{\mu }_0}^2\right)\left[\left(k+r_1\right)e^{k\tau }+\left(k-r_1\right)e^{-k\tau }\right]}\hfill \\ \approx \frac{a\tau \left(1+\frac{k^2{\tau }^2}{6}\right)+\left(r_3-a\right)\left(1+\frac{k^2{\tau }^2}{2}\right)}{\left(1+\frac{k^2{\tau }^2}{2}\right)+r_1\tau \left(1+\frac{k^2{\tau }^2}{6}\right)}=\frac{\left[-2+3\left(1-g\right)\tau \right]\left(1+\frac{k^2{\tau }^2}{6}\right)-\frac{4k^2{\tau }^2}{6}}{\left[4+3\left(1-g\right)\tau \right]\left(1+\frac{k^2{\tau }^2}{6}\right)+\frac{{8k}^2{\tau }^2}{6}}\hfill \\ \approx \frac{\left[-2+3\left(1-g\right)\tau \right]\left[1+\frac{k^2{\tau }^2}{6}-\frac{4k^2{\tau }^2}{6\ast 3\left(1-g\right)\tau }\right]}{\left[4+3\left(1-g\right)\tau \right]\left[1+\frac{k^2{\tau }^2}{6}+\frac{{8k}^2{\tau }^2}{6\ast 3\left(1-g\right)\tau }\right]}\approx \frac{-2+3\left(1-g\right)\tau }{4+3\left(1-g\right)\tau }\left[1-2\frac{k^2{\tau }^2}{3\left(1-g\right)\tau }+O\left(k^4{\tau }^4\right)\right]\hfill \end{array}$

$=\frac{-2+3\left(1-g\right)\tau }{4+3\left(1-g\right)\tau }\left[1-2\left(1-\omega \right)\tau \right]=R\left({\mu }_0=1,\mu =1,\omega =1\right)\left(1-2{\sigma }_{abs}H\right).(15)$

Thus, when the absorption optical depth $\left(\mathbf{1}-\mathit{\omega }\right)\mathit{\tau }$ approaches zero, the reduction of snow bidirectional reflectance R due to absorption is proportional to twice the absorption optical depth of the snow layer.

$2H=1-\frac{1}{{\sigma }_{abs}}\frac{R\left({\mu }_0=1,\mu =1,\omega \to 1\right)}{R\left({\mu }_0=1,\mu =1,\omega =1\right)}(16)$

From Eq. 12, $<L>=1-\frac{1}{{\sigma }_{abs}}\frac{R\left({\mu }_0=1,\mu =1,\omega \to 1\right)}{R\left({\mu }_0=1\mu =1,,\omega =1\right)},$ thus,

$<\mathit{L}>=\mathbf{2}\mathit{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 ${\mathit{\sigma }}_{\mathit{a}\mathit{b}\mathit{s}}\mathit{\tau }<\mathbf{1}.$ Absorption reduces the backscatter signals, $\mathit{I}\left(\mathit{L},\mathit{\omega }<\mathbf{1}\right)=\mathit{I}\left(\mathit{L},\mathit{\omega }=\mathbf{1}\right){\mathit{e}}^{-{\mathit{\sigma }}_{\mathit{a}\mathit{b}\mathit{s}}\mathit{L}}$. To apply it to snow measurements, an absorption correction is needed when <L> is computed,

$<L>=\frac{{\int }_0^{\infty }I\left(L,\omega =1\right)LdL}{{\int }_0^{\infty }I\left(L,\omega =1\right)dL}=\frac{{\int }_0^{\infty }I\left(L,\omega <1\right)e^{{\sigma }_{abs}L}LdL}{{\int }_0^{\infty }I\left(L,\omega <1\right)e^{{\sigma }_{abs}L}dL}.(18)$

3 Snow density measurements from Raman lidar

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_{Raman}$, the peak Raman backscatter of snow near the top of the snow from a lidar measurement is $C_{Raman}{N}_m$. Here, $N_m$ is the number of molecules the laser light will interact with in a unit volume. $N_m=s\frac{{\rho }_{snow}}{18}\ast 6.022\ast {10}^{23}$, where ${\rho }_{snow}$ 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\approx 0.84\ast \frac{\pi d^3}{4}/\frac{\pi d^3}{6}\approx 1.26$. Thus, the peak Raman backscatter cross coefficient at the top of the snow, ${\beta }_{Raman}$, is a function of snow density,

${\beta }_{Raman}\approx 1.26\hspace{0.17em}\ast \hspace{0.17em}{C}_{Raman}\hspace{0.17em}\ast \hspace{0.17em}\frac{{\rho }_{snow}}{18}\hspace{0.17em}\ast \hspace{0.17em}6.022\hspace{0.17em}\ast \hspace{0.17em}{10}^{23}(19)$

where $C_{Raman}\approx 4.2\ast {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⁻¹).

4 Snow grain size and snow density measurements from backscatter lidar

4.1 Snow density can be derived from extinction coefficient and snow grain size

Similar to previous studies (e.g., Barkstrom, 1972; Bohren and Barkstrom, 1974), snow particles are assumed to be spherical. The diffraction truncated extinction coefficient is

${\sigma }_{ext}=\frac{N\pi d^2}{4}=\frac{3}{2d}\frac{{\rho }_{snow}}{{\rho }_{ice}}.(20)$

Here, ${\rho }_{snow}$ and ${\rho }_{ice}$ are the densities of snow and ice, respectively, N is snow particle number density (#/m³), and d is the effective diameter of snow particles. The bulk ice density, ${\rho }_{ice}$, is 916.8 kg/m³. If the average diameter of snow particles is around 0.5 mm and snow density about 300 kg/m³, the diffraction-peak-truncated snow particle extinction coefficient, ${\sigma }_{ext}$, is roughly 1,000 per meter.

Snow density, ${\rho }_{snow}$, can be derived from lidar measurements, because the diffraction-peak-truncated extinction coefficient, ${\sigma }_{ext},$ and grain size, d, can be measured from lidar measurements,

${\rho }_{snow}=\frac{2}{3}d{\sigma }_{ext}{\rho }_{ice}.(21)$

4.2 Lidar measurements of extinction coefficients of snow

The extinction coefficient, ${\sigma }_e,$ can be accurately estimated from the snow backscattering profiles of lidar measurements using a laser at a relatively weak absorbing wavelength (e.g., wavelength $\lambda <0.55\mu m).$ At such a wavelength, the distribution of the absorption-free path lengths between photons entering and exiting the snow can be derived from the lidar backscattering profile with a simple correction of snow absorption.

Previous random-walk studies (Blanco and Fournier, 2006; Hu et al., 2022) suggest that there is a simple relationship between the second moment of the snow entry-to-exit path length distribution, $<L^2>$, and the mean free path length of photons traveling between two snow particles, $<\lambda >$,

$<\lambda ><L^2>={\left[\frac{<L>}{2}\right]}^3(22)$

The diffuse extinction coefficient, $\left(1-g\right){\sigma }_{ext}=1/\lambda$, can be derived from the lowest two moments of the lidar backscattering profile,

${\left(1-g\right)\sigma }_{ext}=8\frac{<L^2>}{<L{>}^3}(23)$

Here, we prove that this equation can also be derived from the two-stream solution, with higher order Taylor expansion of the exponential terms,

$e^{\pm k\tau }=1+{\displaystyle \sum _{n=1}^5}\frac{1}{n!}{\left(\pm k\tau \right)}^n.(24)$

$R\left({\mu }_0=1,\mu =1,\left(1-\omega \right)\to 0\right)\approx \frac{\omega \left[\left(1-k{\mu }_0\right)\left(a+k{r}_3\right)e^{k\tau }-\left(1+k{\mu }_0\right)\left(a-k{r}_3\right)e^{-k\tau }\right]}{\left(1-k^2{{\mu }_0}^2\right)\left[\left(k+r_1\right)e^{k\tau }+\left(k-r_1\right)e^{-k\tau }\right]}\approx \frac{-2+3\left(1-g\right)\tau }{4+3\left(1-g\right)\tau }\left[1-\frac{2k^2{\tau }^2}{3\left(1-g\right)\tau }\left(1-\frac{k^2{\tau }^2}{12}\right)+o\left(k^6{\tau }^6\right)\right]$

$=\frac{-2+3\left(1-g\right)\tau }{4+3\left(1-g\right)\tau }\left[1-2\left(1-\omega \right)\tau +\frac{1}{2}{\left(1-\omega \right)}^2{\tau }^2\left(1-g\right)\tau \right]$

$=R\left({\mu }_0=1,\mu =1,\omega =1\right)\left[1-2{\sigma }_{abs}H+\frac{1}{2}{{\sigma }_{abs}}^2{H}^3\left(1-g\right){\sigma }_{ext}\right]$

$=R\left({\mu }_0=1,\mu =1,\omega =1\right)\times \left[1-{\sigma }_{abs}<L>+{{\sigma }_{abs}}^2\frac{{<L>}^3}{16}\left(1-g\right){\sigma }_{ext}\right].(25)$

Eq. 14 can be rewritten as

$\frac{{<L>}^3}{8}\left(1-g\right){\sigma }_{ext}=\frac{2}{{{\sigma }_{abs}}^2}\left[\frac{R\left({\mu }_0=1,\mu =1,1-\omega \to 0\right)}{R\left({\mu }_0=1,\mu =1,\omega =1\right)}-\left(1-{\sigma }_{abs}<L>\right)\right].(26)$

The right-hand side of the above equation equals $<L^2>.$ Thus,

$<L^2>=\frac{{<L>}^3}{8}\left(1-g\right){\sigma }_{ext,}(27)$

$\left(1-g\right){\sigma }_{ext}=8\frac{{<L}^2>}{{<L>}^3}=\frac{{<L}^2>}{{<H>}^3}.(28)$

4.3 Lidar measurements of snow grain size and snow density

Snow grain size, d, can be derived from lidar measurements of snow bidirectional reflectance, which is the integrated lidar backscattering of the snow layer. As discussed earlier, lidar can directly measure the snow bidirectional reflectance, $R\left({\mu }_0\right)={\int }_0^{\infty }I\left(L\right)dL$. Snow bidirectional reflectance and albedo are highly correlated with the single scattering albedo of the snow particles, which is a function of snow grain size, asymmetry factor of scattering phase function and refractive index of ice, and relatively insensitive to changes in snow density (Barkstrom, 1972; Bohren and Beschta, 1979).

Snow grain size can be estimated from multi-wavelength lidar measurements of snow bidirectional reflectance, using the relationship between snow grain size and snow bidirectional reflectance from the two-stream solution.

First, snow grain size, d, can be measured by a lidar with an infrared wavelength laser (e.g., 1,064 nm or 1,030 nm), for which the absorption is significantly higher than in the visible range. The derivation of the simple relationship between the snow bidirectional reflectance and snow grain size is as follows:

The absorption coefficient of snow particles can be estimated using a simple ray tracing technique (Barkstrom, 1972; Bohren and Barkstrom, 1974),

${\sigma }_{abs}=0.84\hspace{0.17em}\ast \hspace{0.17em}{a}_i\hspace{0.17em}\ast \hspace{0.17em}d\hspace{0.17em}\ast \hspace{0.17em}\frac{N\pi d^2}{4}=0.84\hspace{0.17em}\ast \hspace{0.17em}{a}_i\hspace{0.17em}\ast \hspace{0.17em}d\hspace{0.17em}\ast \hspace{0.17em}{\sigma }_{ext}.(29)$

Here, N is the number density of snow particles. The absorption coefficient of ice, $a_i$, at 1 μm wavelength is around 20 (1/m). For snow particles with a diameter around 0.5 mm, the co-albedo (1 –ω) of laser light at 1 μm wavelength is

$1-\omega \approx 0.84a_{i}d\approx 0.84\hspace{0.17em}\ast \hspace{0.17em}20\left(\frac{1}{m}\right)\hspace{0.17em}\ast \hspace{0.17em}0.0005\left(m\right)\approx 0.00084.(30)$

The asymmetry factor of the single scattering phase function of snow particles, after removing the diffraction peak, can be assumed to be 0.874 (Bohren and Barkstrom, 1974),

$g\approx 0.874,k=\sqrt{3\left(1-\omega \right)\left(1-g\omega \right)}\approx 0.02.(31)$

For an infrared wavelength (around 1,064 nm) lidar pointing near nadir (${\mu }_0=1$) at a thick snow layer (e.g., H > 0.3 m, $\tau >300$),

$2k\tau >6,{\mathrm{e}}^{-2\mathrm{k}\mathrm{\tau }}\approx 0,e^{-\frac{\tau }{{\mu }_0}}\approx 0.(32)$

Thus, the bidirectional reflectance at the top of the layer is

$R\left({\mu }_0=1,{\omega }_{1064nm}\right)=\frac{\omega \left[\left(1-k{\mu }_0\right)\left(a+k{r}_3\right)e^{k\tau }-\left(1+k{\mu }_0\right)\left(a-k{r}_3\right)e^{-k\tau }-2k\left(r_3-a{\mu }_0\right)e^{-\frac{\tau }{{\mu }_0}}\right]}{\left(1-k^2{{\mu }_0}^2\right)\left[\left(k+r_1\right)e^{k\tau }+\left(k-r_1\right)e^{-k\tau }\right]}$

$\approx \frac{\left[\left(1-k{\mu }_0\right)\left(a+k{r}_3\right)-\left(1+k{\mu }_0\right)\left(a-k{r}_3\right)e^{-2k\tau }-2k\left(r_3-a{\mu }_0\right)e^{-\frac{\tau }{{\mu }_0}-k\tau }\right]}{\left(1-k^2{{\mu }_0}^2\right)\left[\left(k+r_1\right)+\left(k-r_1\right)e^{-2k\tau }\right]}$

$\approx \frac{a+k{r}_3}{\left(1+k{\mu }_0\right)\left(k+r_1\right)}=\frac{\frac{3}{4}\left(1-g\right)+\frac{1}{4}\left(1-3g\right)k}{\left(1+k{\mu }_0\right)\left[k+\frac{3}{4}\left(1-g\right)\right]}\approx 1-\frac{1+\frac{3}{4}\left(1-g\right)-\frac{1-3g}{4}}{\frac{3}{4}\left(1-g\right)}k\approx 1-\frac{2k}{1-g}=1-\frac{2\sqrt{3\left(1-\omega \right)\left(1-g\omega \right)}}{1-g}\approx 1-\sqrt{\frac{12\ast 0.84a_i\mathit{d}}{1-g}}\approx 1-A\sqrt{a_i\mathit{d}}.\text{(33)}$

Here, A is dependent on the asymmetry factor, $A=\sqrt{\frac{12\ast 0.84}{1-g}}$. For spherical particles with g = 0.874 (Bohren and Barkstrom, 1974), $A\approx 9.$ $R\left({\mu }_0=1,{\omega }_{1064nm}\right)\approx 1-9\sqrt{a_i\mathit{d}}$.

Spherical particles absorb more strongly than real snow particles, which are non-spherical (Dang et al., 2016). For non-spherical particles with $g\approx 0.75$ (Dang et al., 2016), $A\approx 5.8.$ Thus, a more realistic relationship between snow bidirectional reflectance and snow grain size at relatively absorbing wavelengths (e.g., 1,064 nm wavelength) is

$R\left({\mu }_0=1,{\omega }_{1064nm}\right)\approx 1-5.8\sqrt{a_i\mathit{d}}.(34)$

Snow density (in unit kg/m³) can be derived from the lower moments of multiple scattering entry-to-exit distance distribution of a conservative medium (532 nm lidar measurements with corrections for absorption) and from the lidar measurements of 1,064 nm bidirectional reflectance (layer integrated backscatter of snow),

${\rho }_{snow}=\frac{2}{3}d{\sigma }_{ext}{\rho }_{ice}\approx \frac{16}{3}{\left[\frac{1-R\left({\mu }_0=1,{\omega }_{1064nm}\right)}{5.8}\right]}^2\frac{{<L}^2>}{{\left(1-g\right)<L>}^3}.(35)$

5 Summary and discussion

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, $\left(1-g\right){\sigma }_{ext}=8\frac{{<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).

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

YoH worked together with all the co-authors to develop the measurement concept. All authors contributed to the article and approved the submitted version.

Funding

The Funding for the lead author was provided by NASA’s ESTO and R&A programs.

Acknowledgments

The authors wish to thank the NASA ICESat-2 program, NASA Remote Sensing Theory program, and NASA ESTO’s IIP program for supporting this research.

Conflict of interest

Author CW was employed by Ball Aerospace & Technologies Corp.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Barkstrom, B. R. (1972). Some effects of multiple scattering on the distribution of solar radiation in snow and ice. J. Glaciol. 11 (63), 357–368. doi:10.1017/s0022143000022334

CrossRef Full Text | Google Scholar

Blanco, S., and Fournier, R. (2006). Short-path statistics and the diffusion approximation. Phys. Rev. Lett. 97 (23), 230604. doi:10.1103/physrevlett.97.230604

PubMed Abstract | CrossRef Full Text | Google Scholar

Bohren, C. F., and Barkstrom, B. R. (1974). Theory of the optical properties of snow. J. Geophys. Res. 79 (30), 4527–4535. doi:10.1029/jc079i030p04527

CrossRef Full Text | Google Scholar

Bohren, C. F., and Beschta, R. L. (1979). Snowpack albedo and snow density. Cold Regions Sci. Technol. 1 (1), 47–50. doi:10.1016/0165-232x(79)90018-1

CrossRef Full Text | Google Scholar

Dang, C., Fu, Q., and Warren, S. G. (2016). Effect of snow grain shape on snow albedo. J. Atmos. Sci. 73 (9), 3573–3583. doi:10.1175/jas-d-15-0276.1

CrossRef Full Text | Google Scholar

Dawson, N., Broxton, P., and Zeng, X. (2017). A new snow density parameterization for land data initialization. J. Hydrometeorol. 18, 197–207. doi:10.1175/JHM-D-16-0166.1

CrossRef Full Text | Google Scholar

Ding, K., Xu, X., and Tsang, L. (2010). Electromagnetic scattering by bicontinuous random microstructures with discrete permittivities. IEEE Trans. Geosci. Remote Sens. 48, 3139–3151. doi:10.1109/tgrs.2010.2043953

CrossRef Full Text | Google Scholar

Hu, Y., Lu, X., Zeng, X., Stamnes, S. A., Neuman, T. A., Kurtz, N. T., et al. (2022). Deriving snow depth from ICESat-2 lidar multiple scattering measurements. Front. Remote Sens. 3, 855159. doi:10.3389/frsen.2022.855159

CrossRef Full Text | Google Scholar

Lu, X., Hu, Y., Zeng, X., Stamnes, S. A., Neuman, T. A., Kurtz, N. T., et al. (2022). Deriving snow depth from ICESat-2 lidar multiple scattering measurements: uncertainty analyses. Front. Remote Sens. 3, 891481. doi:10.3389/frsen.2022.891481

CrossRef Full Text | Google Scholar

Meador, W. E., and Weaver, W. R. (1980). Two-stream approximations to radiative transfer in planetary atmospheres: A unified description of existing methods and a new improvement. J. Atmos. Sci. 37 (3), 630–643. doi:10.1175/1520-0469(1980)037<0630:tsatrt>2.0.co;2

CrossRef Full Text | Google Scholar

National Academies (2018). Thriving on our changing planet: A decadal strategy for Earth observation from space. Washington, DC, USA: National Academies Press, 716. doi:10.17226/24938

CrossRef Full Text | Google Scholar

Pershin, S. M., Lednev, V. N., Klinkov, V. K., Yulmetov, R. N., and Bunkin, A. F. (2014). Ice thickness measurements by Raman scattering. Opt. Lett. 39 (9), 2573–2575. doi:10.1364/ol.39.002573

PubMed Abstract | CrossRef Full Text | Google Scholar

Reichardt, J., Knist, C., Kouremeti, N., Kitchin, W., and Plakhotnik, T. (2022). Accurate absolute measurements of liquid water content (LWC) and ice water content (IWC) of clouds and precipitation with spectrometric water Raman lidar. J. Atmos. Ocean. Technol. 39 (2), 163–180. doi:10.1175/jtech-d-21-0077.1

CrossRef Full Text | Google Scholar

Sevetlidis, V., and Pavlidis, G. (2019). Effective Raman spectra identification with tree-based methods. J. Cult. Herit. 37, 121–128. doi:10.1016/j.culher.2018.10.016

CrossRef Full Text | Google Scholar

Stamnes, K., Thomas, G., and Stamnes, J. (2017). Radiative transfer in the atmosphere and ocean. Cambridge: Cambridge University Press. doi:10.1017/9781316148549

CrossRef Full Text | Google Scholar

Wang, Z., Whiteman, D. N., Demoz, B. B., and Veselovskii, I. (2004). A new way to measure cirrus cloud ice water content by using ice Raman scatter with Raman lidar. Geophys. Res. Lett. 31 (15), L15101. doi:10.1029/2004gl020004

CrossRef Full Text | Google Scholar

Xiong, C., and Shi, J. (2013). Simulating polarized light scattering in terrestrial snow based on bicontinuous random medium and Monte Carlo ray tracing. J. Quantitative Spectrosc. Radiat. Transf. 133, 177–189. doi:10.1016/j.jqsrt.2013.07.026

CrossRef Full Text | Google Scholar

Appendix A

The relationship $<L>=2H$ can also be derived as the following:

$I\left(L,\omega <1\right)=I\left(L,\omega =1\right)e^{-{\sigma }_{a}L},(\mathrm{A1})$

$-\frac{d{\int }_0^{\infty }I\left(L,\omega \to 1\right)dL}{d{\sigma }_a}={\int }_0^{\infty }I\left(L,\omega =1\right)L{e}^{-{\sigma }_{a}L}dL={\int }_0^{\infty }I\left(L,\omega \to 1\right)LdL,(\mathrm{A2})$

thus,

$-\frac{d\left\{\log \left[{\int }_0^{\infty }I\left(L,\omega \to 1\right)dL\right]\right\}}{d{\sigma }_a}=\frac{{\int }_0^{\infty }I\left(L,\omega \to 1\right)LdL}{{\int }_0^{\infty }I\left(L,\omega \to 1\right)dL}=<L>(\mathrm{A3})$

From Eqs 5, 7,

${\int }_0^{\infty }I\left(L,\omega \to 1\right)dL=\frac{\omega \left[\left(1-k{\mu }_0\right)\left(a+k{r}_3\right)e^{k\tau }-\left(1+k{\mu }_0\right)\left(a-k{r}_3\right)e^{-k\tau }-2k\left(r_3-a{\mu }_0\right)e^{-\tau /{\mu }_0}\right]}{\left(1-k^2{{\mu }_0}^2\right)\left[\left(k+r_1\right)e^{k\tau }+\left(k-r_1\right)e^{-k\tau }\right]}(\mathrm{A4})$

As $k^2=3\left(1-\omega \right)\left(1-g\omega \right)\ll 1$, ${\sigma }_a=\frac{\left(1-\omega \right)\tau }{H}=\frac{k^2\tau }{3\left(1-g\omega \right)H}\approx \frac{k^2\tau }{3\left(1-g\right)H},$ and $d{\sigma }_a=\frac{2k\tau }{3\left(1-g\right)H}dk$, $e^{-\frac{\tau }{{\mu }_0}}\approx 0,$ it is not hard to derive that

$-\frac{d\left\{\log \left[{\int }_0^{\infty }I\left(L,\omega \to 1\right)dL\right]\right\}}{d{\sigma }_a}=2H.(\mathrm{A5})$

From Eqs A3, A5, $<L>=2H.$