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²>/k_sd = H³, <L³>/k_sd ² = H⁵ (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> (= ∫ 0 ∞ L p ( L ) d L ) and snow depth, and 3) evaluate the relationship among second and third moments of the path length distribution < L²> = ∫ 0 ∞ L 2 p ( L ) d L , <L³> = ∫ 0 ∞ L 3 p ( L ) d L , mean free path f (f = 1/k_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 (Hu et al., 2001). Polarization is removed from the model to reduce computational time. The Monte Carlo code uses random walk processes to mimic photon pulses (hereafter referred to as photons for short) propagating in a turbid medium. The characteristics of photons include location (x,y,z) and direction of propagation ( μ = cos ⁡ θ , ϕ ), where θ and ϕ are the zenith and azimuth angles of the direction, respectively, and a variable L denoting its travel distance. The Monte Carlo code repeats the following steps for each photon until it exits the snow layer in the backward direction and enters the receiver field-of-view:

1) Finding the distance D, the photon will travel before it interacts with a snow particle; D = f ∗ l o g ( ζ ) , where ζ is a random number with uniform distribution between 0 and 1;

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 β ( z ) of lidar measurements by adding absorption to the path length distribution: β ( z ) = p ( L = 2 z ) exp ( − k a L ) d L , where k a is the absorption coefficient of snow.

Figure 1 shows the simulated absorption-free snow backscatter profiles for 1) isotropic scattering (blue line) and for a diffuse scattering coefficient k s d = k s = 1 f   = 200 m⁻¹; 2) anisotropic scattering with a Henyey–Greenstein phase function (g = 0.7, green line) (Henyey and Greenstein, 1941), and with a diffuse scattering coefficient of k s d = ( 1 − g ) k s = 1 f   = 200 m⁻¹; 3) anisotropic scattering with a Henyey–Greenstein phase function (g = 0.88, red line) and diffuse scattering coefficient of 200 m⁻¹. Here, k_s is the scattering coefficient, P ( Θ ) is single-scattering phase function, Θ is the scattering angle between incident light and scattered light, and g is the asymmetry factor, g = ∫ − 1 1 P ( Θ ) c o s ΘdcosΘ . The lidar backscatter profiles from the aforementioned three different types of scattering phase functions align very well, which suggests that backscatter path length distributions are insensitive to the details of the single-scattering phase functions as long as the diffuse scattering coefficients k s d = ( 1 − g ) ∗ k s are the same. This is in agreement with the general invariance property of diffusive random walks, supporting that the general theory applies to ICESat-2.

For snow, the mean free path of diffuse scattering is normally on the order of subcentimeter (Kokhanovsky and Zege, 2004). Thus, the diffuse scattering coefficient k_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 ( k s d H ≫ 4 ). We can fit the lidar backscatter path length distribution using a simple Γ distribution (magenta line in Figure 1), p ( L ) = [ 2 H ( k s d H 4 − 1 ) ] − ( k s d H 4 − 1 ) − 1 Γ ( ( k s d H 4 − 1 ) − 1 ) L ( k s d H 4 − 1 ) − 1 − 1 ⁡ exp [ − L 2 H ( k s d H 4 − 1 ) ] = Γ ( L , α , β   ) = β α Γ ( α ) L α − 1 exp ( − β L )   . (1) Here, α = ( k s d H 4 − 1 ) − 1 ; β = [ 2 H ( k s d H 4 − 1 ) ] − 1 .   Both α and β are far less than 1.

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 (Figure 2A) < L >   = α / β = 2 H , (2) and the second moment of the path length distribution, <L²>, is (Figure 2B), < L 2 >   = α ( α + 1 ) β 2 = k s d H 3 . (3) Thus, scattering optical depth of the snow can be estimated from the higher moments of the path length distribution: τ s = k s d H =   < L 2 > H 2 = 4   < L 2 > < L > 2 . It is also possible to estimate the diffuse scattering coefficient k_sd, or equivalently the transport mean free path f/(1-g) = 1/k_sd, using the first and second moments <L> and <L²> of the path length distribution. τ s = k s d H =   < L 2 > H 2 = 4   < L 2 > < L > 2 ,             1 − g f = k s d = 8 < L 2 > < L > 3 . (4) The Monte Carlo simulation also suggests that the third-order moment of the path length distribution is a simple function of diffuse scattering coefficient and snow depth (Figure 2C), < L 3 >   = α ( α + 1 ) ( α + 2 ) β 3 ≈ k s d 2 H 5   . (5) Thus, snow depth can be derived from all three lower order moments, respectively, In addition, we also found that < L N >   = k s d N − 1 H 2 N − 1   for all positive integers N that is greater than 1.

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 ∫ 0 ∞ p ( L ) ln [ p ( L ) ] d L under the constraints of fixed <L> and fixed < ln(L)>: < L > = α β = 2 H , and < ln ( L ) > ≈ < ∑ N c N L N > = f ( k s d ) = d l n [ Γ ( α ) ] d α − d l n [ Γ ( β ) ] d β = ψ ( α ) − ψ ( β ) . Here, ψ ( x ) is the so-called digamma function.

To better understand the analytical results in Section 2.1, here, we use a different theoretical explanation of the average path length <L> and snow depth H relationship <L> = 2H using the weak absorption limit simulations from both Monte Carlo and DISORT (DIScrete Ordinate Radiative Transfer) (Stamnes et al., 1988) calculations, a general and versatile plane-parallel radiative transfer program.

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, A ( k a → 0 ) = 1 − R ( k a ) R ( k a = 0 ) = ∫ 0 ∞ p ( L ) ( 1 − e − k a L ) d L   ≈ ∫ 0 ∞ p ( L ) k a L d L =   < L > k a , (6) where R (ka = 0) is the reflectance for a non-absorbing medium.

Both the Monte Carlo simulations (Figure 3A) and DISORT (DIScrete Ordinate Radiative Transfer) calculations (Figure 3B) suggest that A ( k a → 0 ) = 1 − R ( k a ) R ( k a = 0 ) ≈ 2 ∗ ( 1 − ω ) τ = 2 ∗ k a ∗ H . (7) Here, ω is single-scattering albedo and τ is total optical depth of snow. 1- ω is the fraction of absorption and (1-ω)*τ is the absorption optical depth, which equals k_a*H. From Eqs, 6,7, we can conclude that < L > k a =   2 ∗ k a ∗ H . Thus, < L > =   2 H .

With the theoretical development in Section 2.1 and Section 2.2, we will use the following iterative procedure to derive snow depth and snow albedo, snow grain size, absorption coefficient, and the diffuse scattering coefficient estimates:

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);

To derive the backscatter path length distribution from the nighttime measurements of ICESat-2 photon heights (ATL03, version 5, Neumann et al., 2019), we need to correct for lidar receiver transient response and absorption by snow. Lidar receiver optic components and laser pulse shape can affect the ICESat-2 snow backscatter profile measurements. The results here are from the strong beam measurements of ICESat-2. To derive snow backscatter path length distribution, the lidar backscatter profile is corrected through a deconvolution using the lidar receiver transient response derived from ICESat-2 measurements of flat, hard targets (Lu et al., 2021).

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, p ( L = 2 z ) e − k a L = β ( z ) / ∫ 0 ∞ β ( z ) d z . (8) p(L) is insensitive to radiometric calibration and atmospheric attenuation. We are testing a couple of algorithms for estimating the absorption coefficient (k_a) from the lidar measurements. For introduction, we first assume the 532-nm absorption coefficient around 0.07 m⁻¹ as a first guess for relatively weakly contaminated snow in the Arctic in springtime (Warren et al., 2006). Thus, snow backscatter path length distribution can be derived from the deconvoluted ICESat-2 snow measurements, p ( L ) = β ( z = L 2 ) e k a L / ∫ 0 ∞ β ( z ) d z . (9) < L n > = ∫ 0 ∞ ( 2 z ) n β ( z ) e 2 k a z d z / ∫ 0 ∞ β ( z ) d z ( n = 1,2,3 , … ) . (10) Figure 4A shows the ICESat-2 snow measurements near the South Pole after correcting for the detector transient response with deconvolution. Snow depths over Antarctica derived from the three methods, <L>/2, (<L²>/k_sd)^1/3, and (<L³>/k² _sd)^1/5 (Figure 4B) agree reasonably well.

Figure 5A is similar to Figure 4A, but from ICESat-2 measurements of snow above first-year sea-ice in the Arctic. The snow depths derived from the multiple scattering photon path length distribution [Figure 5B. Method 1: <L>/2; Method 2 (<L²>/k_sd)^1/3; Method 3: [<L³>*k_sd ²]^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 Figure 5 are snow depths from collocated AMSR-2 snow depth data product (Rostosky et al., 2018). The spatial resolution of the AMSR-2 snow depth data product is 25 km.

Snow grain size (blue line in Figure 6A), diffuse flux attenuation coefficient k_d (blue line in Figure 6B), and snow albedo (red lines in Figure 6) are estimated from the multiple scattering signal of the Arctic ICESat-2 flight as in Figure 5 and theoretical radiative transfer, as discussed in the previous section. The estimated diffuse flux attenuation coefficients are consistent with diffuse attenuation obtained from in situ measurements (e.g., Schwerdtfeger and Weller, 1977).

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) (Räisänen et al., 2015) at the 532-nm laser wavelength. More studies are needed to examine the sensitivity of the technique to lower layers of multiyear snow with reduced scattering coefficients due to melting, infiltration of rainwater, and retention of snowmelt (Boone and Etchevers, 2001). More theoretical studies and field measurements concerning ice sheets are needed to examine if the snow depth we derived here is the depths of first-year snow only.