^{1}

^{1}

^{2}

^{1}

^{3}

^{3}

^{4}

^{5}

^{6}

^{1}

^{1}

^{1}

^{1}

^{1}

^{1}

^{2}

^{7}

^{7}

^{7}

^{8}

^{8}

^{3}

^{1}

^{2}

^{3}

^{4}

^{5}

^{6}

^{7}

^{8}

This article was submitted to Satellite Missions, a section of the journal Frontiers in Remote Sensing

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Snow is a crucial element in the Earth’s system, but snow depth and mass are very challenging to be measured globally. Here, we provide the theoretical foundation for deriving snow depth directly from space-borne lidar (ICESat-2) snow multiple scattering measurements for the first time. First, based on the Monte Carlo lidar radiative transfer simulations of ICESat-2 measurements of 532-nm laser light propagation in snow, we find that the lidar backscattering path length follows Gamma distribution. Next, we derive three simple analytical equations to compute snow depth from the average, second-, and third-order moments of the distribution. As a preliminary application, these relations are then used to retrieve snow depth over the Antarctic ice sheet and the Arctic sea ice using the ICESat-2 lidar multiple scattering measurements. The robustness of this snow depth technique is demonstrated by the agreement of snow depth computed from the three derived relations using both modeled data and ICESat-2 observations.

Snowpack is one of the most important wintertime land-surface characteristics (

While lidar altimetry has been used in airborne measurements of snow depth as the difference between the measured snow top and snow-free ground heights (

The scientific and technological question is as follows: can space-borne lidar measure snow depth (and SWE) over land and over sea ice? If yes, the latter measurement in combination of the measurement of snow top over sea ice would also enable the global estimate of sea ice thickness, which is of substantial importance in the Earth system and national defense (e.g.,

The purpose of this study is to combine radiative transfer theory with Monte Carlo simulations to derive analytical relations of snow depth and optical properties directly from the space-borne lidar measurements of vertical backscatter profiles of snow. This would provide the theoretical foundation for space-borne lidar measurements of snow depth for the first time. We will address the robustness of our method by assessing the agreement among three derived relations and the convergence of the relations using a different approach. As a preliminary application, these relations are then used to retrieve snow depth over the Antarctic ice sheet and Arctic sea ice from ICESat-2 lidar measurements. The retrieval of SWE and comparison of SWE and snow depth with observations will be reported in separate articles.

In this section, we will introduce the simple relationships among the average path length, snow depth, and diffuse scattering coefficient.

Inspired by biological studies of the spatial distributions of ant colonies (

At 532 nm, pure snow is a weakly absorbing, diffusive scattering medium. Photons bounce inside the snow many times before exiting the snow. Thus, the high vertical resolution measurements from the ICESat-2 mission (^{2}) and the average path length <L> hold for pencil beam point source lidar measurements as it does for uniform incidence? 2) What are the relationships among the higher moments (<L^{2}>, <L^{3}>) of the path length distribution, snow depth, and diffuse scattering mean free path f (1/scattering coefficient k_{d})?

Here, we try to answer the aforementioned equations through radiative transfer simulations and introduce an algorithm to derive snow depth directly from space lidar measurements of the scattering path length distributions. For the first time, snow depths are estimated directly from lidar measurements of vertical backscatter profiles of snow. These snow depth estimates are different from the existing laser altimetry measurements (

The backscatter path length is the distance a photon has traveled inside the snow. It starts when a photon enters the snow and ends when the same photon exits the snow in the backward direction. In this subsection, we use Monte Carlo simulations to evaluate the validity of <L> = 2H, <L^{2}>/k_{sd} = H^{3}, <L^{3}>/k_{sd}
^{2} = H^{5} (where k_{sd} is the diffuse scattering coefficient) and briefly introduce the snow depth retrieval algorithms.

Specifically, we will develop 1) a simple Monte Carlo simulation model to calculate the backscatter path length distribution p(L); 2) evaluate the relationship between average path length <L> (= ^{2}> = ^{3}> = _{sd}), and snow depth H.

The Monte Carlo simulation code we developed for this study is based on a semi-analytic Monte Carlo method that is similar to the lidar radiative transfer model developed for CALIPSO cloud measurements (

1) Finding the distance D, the photon will travel before it interacts with a snow particle;

2) Estimating the new location of scattering interaction (x,y,z) from the old location and traveling direction and adding D to the photon traveling length L_{0};

3) Estimating the probability of the photon scattered directly to the lidar receiver and adding that probability to the backscatter path length distribution while registering the altitude as L′ = L_{0} + z;

4) Calculating the new traveling direction (

We can choose between two types of surfaces below the snow: a) a sea-ice type of surface with specular reflectance following Fresnel’s theory; b) land surfaces with Lambertian reflectance. Each Monte Carlo calculation will take about 10 million photons to generate the backscatter path length distribution p(L). We can also calculate the attenuated backscatter profile

^{−1}; 2) anisotropic scattering with a Henyey–Greenstein phase function (g = 0.7, green line) (^{−1}; 3) anisotropic scattering with a Henyey–Greenstein phase function (g = 0.88, red line) and diffuse scattering coefficient of 200 m^{−1}. Here, k_{s} is the scattering coefficient,

Monte Carlo simulation results of the lidar backscatter profile for isotropic (blue line) and anisotropic scattering with Henyey–Greenstein phase function (green: g = 0.7; red: g = 0.88) for snow depth 1 m and the same effective scattering coefficient k_{sd} = (1-g)k_{s}.

For snow, the mean free path of diffuse scattering is normally on the order of subcentimeter (_{sd} is of the order of a few hundreds per meter. For snow depth H greater than 10 cm, the diffuse scattering optical depth is far greater than 4 (

Monte Carlo simulations suggest that for snow with a variety of optical properties such as single-scattering phase functions and scattering coefficients, the average path length of the path length distribution is twice that of snow depth H (^{2}>, is (_{sd}, or equivalently the transport mean free path f/(1-g) = 1/k_{sd}, using the first and second moments <L> and <L^{2}> of the path length distribution.

^{2}> and <L^{3}> are the second and third moments of the path length distribution, respectively; k_{d} is the effective scattering coefficient; and f is the mean free path.

Monte Carlo simulations using various combinations of snow depth H and diffusive scattering coefficient k_{sd} suggest that the multiple scatter signals measured by a space-based lidar, such as ICESat-2, always follow Γ-distribution p(L) = Γ(L,α,β), which is the distribution of maximum entropy

To better understand the analytical results in

When the absorption coefficient (k_{a}) approaches zero, the total absorption of the snow layer, A, will reduce the integrated backscatter of a non-absorbing medium by <L>*k_{a},

Both the Monte Carlo simulations (_{a}*H. From Eqs,

Monte Carlo simulations _{a}*H. As _{a}. Thus, <L> = 2H.

With the theoretical development in

1) The ICESat-2 attenuated backscatter profile of snow, β(z), is improved with a correction of instrument transient response: de-convolution of the ICESat-2 snow backscatter profile using the instrument transient response (Lu at al., 2021);

2) An initial guess of the snow absorption coefficient k_{a} (e.g., k_{a} = 0.07 m^{−1}) is made;

3) The absorption-free, snow backscattering path length distribution, p(L), is estimated by removing the impact of absorption from the ICESat-2 attenuated backscatter profile of snow: p (L = 2z) = β(L)exp (k_{a}L); and snow albedo

4) Snow depth H from path length distribution is derived:

5) The diffuse scattering coefficient

6) Snow particle grain size R and diffuse flux attenuation coefficient, _{
d
}
_{
a
} (_{
d
} is a function of diffuse scattering coefficient _{
sd
} and absorption coefficient:

7) A new absorption coefficient, _{
sd
} of step 5 and diffuse flux attenuation coefficient k_{d} of step 6:

8) If abs (k_{a}′-k_{a}) > 0.001, k_{a} is replaced by k_{a}′ and procedures 3–7 are repeated until

To derive the backscatter path length distribution from the nighttime measurements of ICESat-2 photon heights (ATL03, version 5,

After de-convolution, the attenuated lidar backscatter profile β(z = L/2), scaled by the snow layer–integrated attenuated backscatter, equals path length distribution p(L) reduced by absorption,_{a}) from the lidar measurements. For introduction, we first assume the 532-nm absorption coefficient around 0.07 m^{−1} as a first guess for relatively weakly contaminated snow in the Arctic in springtime (^{2}>/k_{sd})^{1/3}, and (<L^{3}>/k^{2}
_{sd})^{1/5} (

^{2}> and diffuse scattering coefficient method: H = (<L^{2}>/k_{sd})^{1/3}; Black: <L^{3}> and diffuse scattering coefficient method: H = (<L^{3}>/k_{sd}
^{2})^{1/5}.

^{2}>/k_{sd})^{1/3}; Method 3: [<L^{3}>*k_{sd}
^{2}]^{1/5};) also show reasonable agreement of ICESat-2 freeboard measurements (altitude difference between sea water and snow top). The freeboard is defined as the total height of the snow cover and sea ice above the ocean. If proven, this technique will enable accurate lidar measurements of both freeboard and snow depth, and thus estimates of sea-ice depth directly without having to rely on radars to provide altitude of the snow/sea-ice interface. The purple data points of

Snow grain size (blue line in _{d} (blue line in

_{d} (blue line) estimated from the 532-nm lidar backscatter measurements.

For snow thicker than 1 m, model simulations suggest that multiple scattering path lengths can last a hundred meters or more. Due to limited satellite downlink bandwidth and the finite number of detected photons, ICESat-2 does not send down all the snow multiple scattering signals. As a result, the long tail parts of the multiple scattering path lengths are truncated, and the amount of photons reaching the deeper part of the snow is under-estimated. We are assessing its impact on snow depth retrievals and testing different extrapolation methods using simulated data.

This technique retrieves the depths of snow layers comprising diffusively scattering snow particulates with large scattering coefficients and small single-scattering co-albedos (absorption coefficient/scattering coefficient) (

Previous studies (

To answer this question, a Monte Carlo lidar radiative transfer model is developed for simulating ICESat-2 measurements of laser light propagation in snow. For normal incidence, the model results confirmed the theory predicting that 1) the average snow backscattering path length <L> = 4V/S derived from the ICESat-2 measurements equaled two times the snow depth H; 2) the diffuse scattering mean free path (1/k_{sd}) times the second moment of the path length distribution; <L^{2}> equals the third power of snow depth (H^{3}); (3)

The path length distribution follows a simple Γ (L, α, β) distribution function and _{sd}), both <L> and <ln(L)> are fixed as well. The multiple scattering photon path length distribution p(L) maximizes entropy

For snow depth greater than a few centimeters,

Initial ICESat-2 data analysis suggests that snow depth can be derived from the average diffuse scattering path length derived from the ICESat-2 lidar measurements. More studies in retrieval algorithm improvements and validations are needed to mature the new retrieval concept. For example, we are working on improving solar background removal algorithms to reduce uncertainties in daytime snow depth retrievals. We are also working on robust along-track averaging algorithms to reduce snow depth retrieval uncertainties associated with snow surface roughness and various sampling errors.

This snow depth retrieval technique may open a new path for future remote sensing of snow globally.

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

YH initiated the theoretical concept and conducted Monte Carlo modeling studies. XL provided ICESAT-2 data analysis under TN and NK guidance. XZ is the science advisor of the research. SS and KS provided DISORT modeling analysis. CW and all co-authors contributed various ideas during discussions.

Authors CW and JL were employed by Ball Aerospace 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.

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.

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. The ICESat-2 ATL03 data used in this study are available at the National Snow and Ice Data Center: