Quantifying Volumetric Scattering Bias in ICESat‐2 and Operation IceBridge Altimetry Over Greenland Firn and Aged Snow

The Ice, Cloud, and Land Elevation Satellite‐2 (ICESat‐2) mission has collected surface elevation measurements for over 5 years. ICESat‐2 carries an instrument that emits laser light at 532 nm, and ice and snow absorb weakly at this wavelength. Previous modeling studies found that melting snow could induce significant bias to altimetry signals, but there is no formal assessment on ICESat‐2 acquisitions during the melting season. We performed two case studies over the Greenland Ice Sheet to quantify bias in ICESat‐2 signals over snow: one to validate Airborne Topographic Mapper (ATM) data against Next Generation Airborne Visible/Infrared Imaging Spectrometer (AVIRIS‐NG) grain sizes, and a second to estimate ICESat‐2 bias relative to ATM. We used snow optical grain sizes derived from ATM and AVIRIS‐NG to attribute altimetry bias to snowpack properties. For the first case study, the mean and standard deviation of optical grain sizes were 340 ± 65 µm (AVIRIS‐NG) and 670 ± 420 µm (ATM). A mean altimetry bias of 4.81 ± 1.76 cm was found for ATM, with larger biases linked to increases in grain size. In the second case study, we found a mean grain size of 910 ± 381 µm and biases of 6.42 ± 1.77 cm (ICESat‐2) and 9.82 ± 0.97 cm (ATM). The grain sizes and densities needed to recreate biases with a model are uncommon in nature, so we propose that additional surface attributes must be considered to characterize ICESat‐2 bias over snow. The altimetry biases are within the accuracy requirements of the ICESat‐2 mission, but we cannot rule out more significant errors over coarse‐grained snow.


Introduction
The Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) was launched in September 2018 to measure surface height over glaciers and ice sheets (Markus et al., 2017).Since then, ICESat-2 data products have been developed to estimate the surface height of land ice, vegetation canopies, and sea ice (Kwok et al., 2019;Neuenschwander & Pitts, 2019;Smith et al., 2019).The sole onboard instrument, the Advanced Topographic Laser Altimeter System (ATLAS), emits laser light at 532 nm and produces fine spatial resolution data (12 m footprint diameter and 10 kHz pulse repetition frequency) and a required accuracy of 0.4 cm yr 1 for ice sheet annual elevation change (Markus et al., 2017).Recent comparisons with ground-based data have shown that the ATLAS laser has a measured accuracy of <4 cm over ice sheet interiors (Brunt et al., 2021(Brunt et al., ). et al., 2022;;Shean et al., 2021).A complication with past snow studies is that forested and mountainous environments have significant seasonal snow, yet these regions are subject to elevation uncertainty in ground-based and airborne surveys.In recent years, digital elevation models (DEMs) from lidar have become common data sets for snow depth estimates (Deems et al., 2013), though current lidar acquisitions are limited to airborne and ground-based surveys.Provided that terrain biases are accounted for, ICESat-2 has the potential to support snow studies through inter-seasonal measurements of terrain height in forests and mountains, though existing efforts are currently limited (Besso et al., 2024;Bormann et al., 2018).
A likely complication with ICESat-2 snow depth retrievals is that a laser shot may experience multiple scattering events within a snow layer before returning to the detector (Perovich, 2007) due to weak absorption of visible light (Warren & Wiscombe, 1980).This phenomenon, which we refer to as "volumetric scattering," is greatest in clean, coarse-grained snow, where the increased path length between individual snow particles will introduce a delay time in the returned laser pulse.The optical grain size of snow, a quantity used to represent snow grains as simplified shapes, is strongly linked to photon path length.Figure 1 illustrates the volumetric scattering problem in a snowpack simulated by a Monte Carlo model (Section 3.1).The path length traveled by photons within a snowpack is similar between 532 and 1,064 nm at small optical grain sizes, but the path lengths at 532 nm increase with grain size.Near-infrared snow reflectance is low in snow with an optical grain size of 500 μm, so fewer 1,064 nm photons propagate into the snowpack.Snow impurities may attenuate the ICESat-2 signal and reduce volumetric scattering bias, though impurity content has significant variability at small spatial and temporal scales (Flanner et al., 2007;Skiles et al., 2017).
Previous studies have assessed the potential impacts of snow on lidar measurements at 532 nm.Harding et al. (2011) found that return waveforms from an airborne 532 nm lidar experienced significant pulse broadening over snow, resulting in range biases on the order of a few centimeters.A modeling study by Kerekes et al. (2012) found that centimeter-level biases occurred most frequently when the optical grain size of snow was 500 μm or more, and the amplitude of received waveforms was low relative to fine-grained snow returns.Smith et al. (2018) simulated ICESat-2 measurements over a snow-covered surface using a suite of surface height estimation techniques.The authors concluded that elevation biases may exceed 0.15 m over regions of clean, coarse-grained snow if the current ICESat-2 height estimation scheme is used for retrieved photons, though biases could decrease if other techniques are used or if snow impurities (e.g., black carbon) are present.At the time of writing, the ICESat-2 mission has collected over 3 years of altimetry measurements over highlatitude regions, yet there have been no documented efforts to quantify volumetric scattering biases over snow.As part of an extensive validation effort, Operation IceBridge (OIB) launched a series of flights over Greenland late in the 2019 melt season.The flights collected elevation measurements using the Airborne Topographic Mapper (ATM), a lidar that operated at 532 and 1,064 nm during the 2019 flights.Near-coincident flights were performed with the Next Generation Airborne Visible/Infrared Imaging Spectrometer (AVIRIS-NG) to retrieve hyperspectral reflectance and the optical grain size of snow.
Here we perform two case studies to assess bias in ICESat-2 and ATM altimetry measurements due to volumetric scattering in snow.Optical grain sizes derived from AVIRIS-NG reflectance data and ATM waveforms serve as input to a Monte Carlo ray tracing model to simulate altimetry bias over the Greenland ablation zone.In parallel, surface heights derived from ICESat-2 and the ATM 532 nm beam are compared to the ATM 1064 nm beam, which we assume also measures the surface, to estimate observed bias.The findings presented here will serve as useful insight for the eventual development of an ICESat-2 bias correction algorithm over snow-covered surfaces.

Case Study Locations and Available Data
We performed two case studies over the Greenland ablation zone, which we refer to as CS1 and CS2 in the remainder of the paper.The locations of the study regions are shown in Figure 2, and Table 1 outlines the data sets used for each.
The first case study (CS1), performed for 6 September 2019, is located at coordinates 75.316-75.438°N, 56.528-56.778°W.This date and location correspond with a significant overlap between ATM and AVIRIS-NG flights, with ∼40 km of OIB flight data overlapping with AVIRIS-NG surveys.The ice surface features many crevasses and refreezing supraglacial lakes during this time of year, several of which were observed by ATM and AVIRIS-NG.There were no significant overlaps with ICESat-2 data over this region, so we used CS1 as a proof-of-concept to demonstrate green light penetration in snow.
The second case study (CS2) was performed for 4 September 2019 at coordinates 78.783-78.807°N,66.066-66.090°W.Across this region, ATM followed an ICESat-2 overpass for 20 min, closely intersecting two ICESat-2 beams.Local topography can have a first-order impact on our analysis (Wang et al., 2019), where ∼10 m of separation between ATM and an ICESat-2 beam may lead to significant differences in elevation estimates, particularly over rough terrain.We therefore limited our analysis to regions where ATM footprints were within the ICESat-2 footprint, or a maximum distance of 12 m (Magruder et al., 2021).This restriction minimized errors due to data separation, but it also limited the analysis to a 2.5 km region over the ice sheet interior (CS2 in Figure 2, red line).

ATM
ATM is an altimetric lidar that has been used for high-latitude elevation measurements since 1993 (Brock et al., 2002;Krabill et al., 2002).In recent years, it has been used to validate ICESat-2 surface height estimates over sea ice and the 88°S transect of Antarctica (Brunt et al., 2021;Kwok et al., 2020) as part of OIB.Originally flown using a 532 nm laser, the instrument suite flown in 2019 was composed of two laser altimeters that featured off-nadir scan angles 2.5°and 15°, which corresponds to swath widths of 40 and 245 m at typical flight altitudes.The 2.5°"narrow swath" altimeter is a dual-color laser that operates at 532 nm (green) and 1,064 nm (nearinfrared) simultaneously.The near-infrared laser has a footprint diameter of 0.91 m, or 40% larger than the 532 nm beam (0.64 m).
Here, we used two Level-1B Narrow-Swath data products: the Elevation and Return Strength with Waveforms (ILNSAW1B) and the Near-Infrared Waveforms (ILNIR1B).Both data products include information about the transmitted and received waveforms, including the amplitude and width of each waveform and the corresponding aircraft-surface range estimates (Studinger & Manizade, 2020a, 2020b).The ranges are derived using the centroid (median) time of the transmitted and received pulses, and these ranges are used with aircraft GPS information and the WGS84 ellipsoid to derive surface elevation (Brock et al., 2002).Brunt, Neumann, and Larsen (2019) found that the 532 nm laser has a mean uncertainty of ∼8.5 cm relative to ground-based measurements over the 88°S transect of Antarctica.However, the signal is sensitive to rough surface features over ice sheets, such as ice cracks or crevasses.Therefore, we applied a moving mean filter with a window size of 500 samples (30 m along-track) to mitigate noise.For the first case study, applying this filter resulted in a data reduction of 17%-21% in the ATM beams, primarily from rough surface features.Smoother ice was present for the second case study, so the data reduction was lower (4%-5%).

ICESat-2
ICESat-2 is a polar orbiting satellite with an operational altitude of 500 km and a 91-day repeat cycle.A single ATLAS 532 nm laser pulse is split into six beams that are configured in pairs, with a 90 m separation between beams within a pair and a 3.3 km separation between pairs (Neumann et al., 2019).Each beam pair includes a strong beam and a weak beam, where the strong beam has four times the laser power of the weak beam to improve coverage of low-reflectivity targets.The beams are named according to their ground track from left to right: GT1L/R, GT2L/R, and GT3L/R.At its operational altitude, ICESat-2 has a surface footprint of 12 m, which provides significant overlap between consecutive footprints given the 0.7 m along-track sampling.
The ATLAS product used here is the ATL03 Global Geolocated Photon Data V005 (Neumann et al., 2020), which consists of latitude, longitude, and surface elevation for received photons.Each tagged photon is also classified as either signal or solar background, based on a statistical confidence algorithm (Neumann et al., 2019).Although noisier than other ICESat-2 data sets, we selected ATL03 to better capture the scattering experienced by individual photons.The number of photon returns was high over the Greenland Ice Sheet, so we only considered photons labeled by the ATL03 algorithm as high-confidence signal photons.Comparisons with ground-based measurements over the 88°S transect show a mean uncertainty of ∼4 cm for ATL03 (Brunt et al., 2021).Due to the conical scanning of the lidar, ATM acquisitions over CS2 alternated between overlapping the central strong beam and the central weak beam of ICESat-2.We selected a 2 km extent where OIB flew inline with the central strong beam (Figure 2c) to ensure that we received a high rate of photons across the study area.
We briefly consider the impacts of volumetric scattering on ATL06, the Land Ice Height Product V005 (Smith et al., 2019).The ATL06 algorithm aggregates geolocated ATL03 photons into 40 m segments, from which a mean surface height is derived.ATL06 segment values are posted every 20 m, yielding a 50% overlap between consecutive segments.Brunt, Neumann, and Smith (2019) found that ATL03 photon-based heights are generally a few centimeters higher than those from ATL06 segments due to differences in the processing algorithms.We only consider high-confidence ATL03 photons with ATM data in close proximity, so additional errors relative to ATL06 (which considers photons of low, medium, and high confidence and corrects for several instrument effects) are expected.A similar analysis could be conducted on surface heights derived from the ICESat-2 SlideRule application, with the advantage of configurable along-track resolution and photon filtering (Shean et al., 2023).We leave a detailed analysis using SlideRule to a future study.

AVIRIS-NG
The Next Generation Airborne Visible/Infrared Imaging Spectrometer is an airborne hyperspectral imager that has been used to retrieve surface radiances since 1986 (Gao et al., 1993;Green et al., 1998).Originally operating at 10 nm spectral resolution, the instrument now observes the Earth's surface between 380 and 2,510 nm at a spectral resolution of 5 nm.The spectrometer was flown at an altitude of 7.5 km over Greenland, resulting in a spatial resolution of approximately 4 m and a swath width of 3.9 km.Surface reflectances are derived from the radiances by applying an atmospheric correction and orthorectification.Reflectances from AVIRIS-NG generally have an accuracy within 2%-5% (Thompson et al., 2019).The spectrometer has been used for a suite of applications since its inception, including vegetation mapping, trace gas identification, and retrieval of snow grain size (Kokaly et al., 2003;Nolin & Dozier, 2000;Thorpe et al., 2016).
We used AVIRIS-NG reflectances for CS1 to derive the optical grain size of snow for comparison against the altimetry data.An inversion algorithm derived by Nolin and Dozier (2000) was used to relate changes in the ice absorption feature at 1,030 nm to changes in optical grain size.In short, the algorithm compares AVIRIS-NG reflectances to those derived from a radiative transfer model to show that optical grain size increases as the near-infrared snow reflectance decreases.The snow is assumed to be composed of spherical ice particles, and snow impurities are assumed to have a negligible influence on snow albedo around the 1,030 nm absorption feature.Although impurity content is assumed to have little influence on retrievals over the regions of interest, we recognize that impurities such as ice algae or cryoconite may impact retrievals over the Greenland ablation zone (Cook et al., 2020).Optical grain sizes derived from this algorithm have an uncertainty of 50 μm, as determined by Nolin and Dozier (2000) and Fair et al. (2022).

ATM
AVIRIS-NG grain sizes were unavailable over CS2, so we instead used an algorithm that infers grain size from ATM data (Smith et al., 2023).Over snow, subsurface scattering affects green ATM waveforms by reducing the maximum amplitude and increasing the width of the received pulse (Smith et al., 2018).The algorithm exploits this occurrence to compare waveforms from the ILNSAW1B product (Section 2.2.1) to an idealized waveform with no subsurface scattering.The grain size is then estimated based on differences in amplitude and pulse width.The model waveform is derived assuming that the snow has no impurities and has a density of 400 kg m 3 .While the former assumption is reasonable over this region of the Greenland Ice Sheet (Flanner et al., 2007), the assumed snow density is higher than typical values (Fausto et al., 2018;Schaller et al., 2016).The algorithm is also more sensitive to subsurface snow properties and scattering, so grain sizes derived by ATM are generally higher than those found using AVIRIS-NG, particularly over regions of coarse-grained snow (Section 4).A rough uncertainty in the recovered grain size is given by the range of grain sizes that can correspond to a range of snow densities.For density values between 200 kg m 3 (fresh snow) and 800 kg m 3 (nearly solid ice), the range of possible grain sizes is approximately a factor of two larger or smaller than that for the prescribed density of 400 kg m 3 .A more detailed description and evaluation of the ATM grain size model is given by Smith et al. (2023), which is in review at the time of publication.

Monte Carlo Modeling
We first estimated altimetry bias using a combination of optical grain size data and Monte Carlo modeling.The model fires photons into a simulated semi-infinite snowpack and records their total path length until they are absorbed or leave the medium (Schneider et al., 2019).The snowpack has user-prescribed optical grain sizes and density, and it is configured to have spherical ice particles with negligible impurity content.The model has additional inputs for solar zenith and particle surface roughness, but we assumed that (a) the snow particles were smooth and (b) the solar zenith angle was equal to the mean solar geometry observed at the time of flight for ATM.The snowpack was also assumed to have a uniform optical grain size.
The Monte Carlo model was used to benchmark lidar delay time within a snowpack.Simulations were conducted for different permutations of photon wavelength, snow density, and optical grain size.The simulations launched 10 5 photons into a snowpack at wavelengths 532 and 1,064 nm to emulate the ATM dual-colored laser interacting with a snow-covered surface.We performed these simulations for grain sizes 50-1,500 μm at 50 μm resolution.
We then applied spline curve fitting to improve the resolution to 1 μm.As a baseline, snow density was prescribed in the model to be consistent with observations by Fausto et al. (2018), that is, 315 kg m 3 .Additional densities were given to account for variability in aging snow layers, such that the full range of densities was ρ s = [100, 200, 300, 315, 400, 500].We obtained the total distance traveled by the photons that escaped from the top of the snowpack, and for each wavelength the median distance traveled by escaped photons was calculated to replicate the reference photon technique employed by ICESat-2 (Neumann et al., 2020).The median travel distance was treated as the surface height offset relative to an unbiased measurement.If we treat the 532 nm distances as analogs for biased surface height measurements (L 532 ) and the 1,064 nm distances as ideal measurements (L 1064 ), then the modeled bias estimate ΔL is simply: In this configuration, a positive ΔL implies that 532 nm photons traveled a greater distance within the snowpack, which would suggest a negative bias in the final surface height estimate.Conversely, the 1,064 nm travel distance (surface height) will be biased high (low) if there is a negative ΔL. Modeled biases were placed into lookup tables depending on the density used in the simulation.The result was six lookup tables that each had 1,451 bias estimates as a function of optical grain sizes 50-1,500 μm, given the snow densities previously defined.The biases in these lookup tables are the errors we would expect if grain size were the only factor impacting ICESat-2 and ATM observations.

Earth and Space Science
10.1029/2022EA002479 FAIR ET AL.

Observed Bias
We look for bias in the altimetry data by comparing 532 nm elevation estimates with those from the ATM 1064 nm beam.The ATM beams periodically did not record laser pulses, so we applied a co-registration algorithm to match data samples from both beams.Because the beams fire simultaneously, the algorithm co-registers shots between beams by using the time stamps recorded for each laser pulse.The co-registered shots were then filtered to match with the central strong beam of ICESat-2, which was selected over the weak beam for its greater overlap with the ATM beams.Assuming that bias is caused exclusively by volumetric scattering, we approximated observed bias between ICESat-2 and ATM elevations using Equation 2: where ΔL obs is the altimetry bias between ICESat-2 (h IS2 ) and the ATM 1064 nm beam (h ATM,1064 ).
To attribute altimetry bias to optical grain size, we co-registered ICESat-2 and corrected ATM laser pulses with AVIRIS-NG or ATM grain size estimates.We mapped each swath of grain size data with an estimate of modeled bias by matching grain sizes with the closest values found in each lookup table.In other words, each swath of coregistered data had six modeled bias estimates for each snowpack density given in Section 3.1.The observed biases (ΔL obs ) were compared to matched model biases (ΔL) at these densities.If the observations agreed with at least one of the modeled results, then we could suggest that (a) the altimetry biases are linked to the optical grain size of snow and (b) the bias is consistent with one of the given snow densities.

ATM Beam Calibration
The accuracy of the ATM beams relative to each other has not been documented, so we performed a bias assessment of the two beams in the absence of snow.Operation IceBridge was flown from Qaanaaq Air Base (formerly Thule) in September 2019, and each flight included an overpass of the aircraft ramp for calibration purposes.We selected the ramp overpass from September 6 (Track 1730, Figure 3a), and Equation 1 was used to estimate bias over a dark, flat surface.Figure 3b shows the differences between the green and NIR beams.
Comparisons over the ramp consistently feature negative bias, implying that NIR ranges surface heights are lower than those of the green beam.The bias has a nearly Gaussian distribution between 8 and 0 cm, with a slightly larger distribution toward less negative values (Figure 3c).Repeat flights over the ramp on different dates yielded similar results, with noisy transmitted pulses in the NIR beam proposed to be a factor.The median bias was 3.85 cm, and we applied this value as a correction factor to all comparisons between the ATM beams in the featured case studies.

Case Study 1
Model-derived results of altimetry bias have a strong dependence on snow optical grain size and density, as seen in Figure 4.At smaller grain sizes, bias is less sensitive to changes in snow density, particularly at grain sizes below 400 μm.Larger grain sizes exhibit greater dependence on snow density, especially when ρ s ≤ 200 kg m 3 .The largest modeled biases occur for ρ s = 100 kg m 3 , up to a maximum of 37.84 cm at the largest grain sizes.However, a combination of large snow grains and low density was not observed over either study location.Densities of ρ s = 315-500 kg m 3 are more representative of Greenland snow and firn layers, which generally corresponded to biases of 7.55-11.88cm.The bias asymptotically approaches zero with decreasing grain size at all densities, implying that little altimetry bias should be expected over fine-grained snow.
The AVIRIS-NG optical grain sizes co-registered with ATM are shown in Figure 5.The southern reaches of CS1 are characterized by grain sizes of ≤200 μm that typically increase near crevassed terrain or near melt ponds.In the northern portions of the flight track, grain sizes increase to 300-400 μm.This increase corresponds with a general decrease in surface elevation (Figure 6), with lower elevations implying warmer temperatures, greater melt, and faster snow metamorphism.Subsurface scattering on the order of 1-10 m is evident throughout the study area (green dots in Figure 6), indicating the presence of heavy crevassing.Although ATM waveform-fitted grain sizes exhibit similar trends to those from AVIRIS-NG, the derived values are much larger, with a mean and standard deviation of 653 ± 422 μm for ATM and 338 ± 65 μm for AVIRIS-NG.Figure 7 shows that the small grain sizes at the start of CS1 correspond to negligible model bias and uncertainty due to snow density.In regions with larger grain sizes, the bias increases, with the full extent dependent on the optical grain sizes used as model input.The lower grain sizes of AVIRIS-NG correspond with a maximum bias of 1.95 ± 0.39 cm, whereas bias peaks at 5.58 ± 1.12 cm with ATM waveform-fitted grain sizes.The green-NIR path length differences generally show agreement with the model when the ATM grain size algorithm is used, with the ATM-derived model estimates accounting for 71% of the observed bias.The best agreement between the model and the observations is in regions of larger optical grain size.The model underestimates bias relative to the observations at smaller grain sizes, suggesting that (a) the observations may agree better with snow densities of ρ s = 200-250 kg m 3 , or (b) other factors are influencing the bias in this region.

Case Study 2
The optical grain sizes and along-track surface heights for CS2 are given in Figures 8 and 9a.The region features gently sloped terrain that decreases in surface height over the track.At the large scale, co-registered ATM 532 nm and ICESat-2 data show general agreement in surface height trends.The ICESat-2 data has slightly larger variability among individual photons which   (Schneider et al., 2019).Bias is given as a function of snow optical grain size and snowpack density.The snow density ρ s = 315 kg m 3 is used to represent the average snow density reported over Greenland by Fausto et al. (2018) (ρ s = 315 kg m 3 ).The simulated bias at ρ s = 300 kg m 3 is nearly identical to that of ρ s = 315 kg m 3 and was thus omitted.

Earth and Space Science
10.1029/2022EA002479 FAIR ET AL. may be attributed to the inherent noisiness of the ATL03 product.The mean bias between the ATM 532 nm heights and ICESat-2 heights is ∼1 cm, with a precision of 10 cm.The grain sizes are comparable to those derived over CS1, with a mean value of 910 ± 381 μm.A surface feature at 1.3 km along-track corresponds with a localized increase in grain size and greater disagreement between the altimeters.Utilizing coincident imagery Figure 6.Surface heights for CS1, as given by ILNSAWL1B ("Raw," green dots).A moving mean filter with a sampling window of 30 m (black line) was applied to remove features significantly below the surface.from the IceBridge Continuous Airborne Mapping By Optical Translator (CAMBOT) system (Figure 9b), we found no obvious signs of a crevasse.Because the discrepancy occurs at a relatively small scale (˜100 m), the ICESat-2 photons may have interacted with a rough surface feature at the meter scale differently than ATM 532 nm.Despite the variability in optical grain size, the modeled bias over CS2 peaks at 6.7 cm at the start of the track before decreasing to ∼5.5 cm as grain size decreases.The mean modeled bias across the region is 4.93 ± 1.89 cm.Between the two altimeters, ICESat-2 bias trends show the closest agreement to modeled estimates.The ICESat-2 biases peak at 17.64 cm at the start of the track before reducing to 6.42 ± 1.77 cm.The ATM green-NIR range differences show weaker agreement with the model than in CS1, with a mean bias of 9.82 ± 0.97 cm.However, ATM trends resemble those seen in the model, implying that optical grain size still has an influence on ATM signals.Overall, the model accounts for 66% (ATM) and 95% (ICESat-2) of the observed bias.The surface feature at 1.3 km produces the greatest disagreement between the altimeters and the model, with ICESat-2 showing the greatest agreement with the 1,064 nm beam and ATM 532 nm having the weakest.It is possible that ICESat-2 and both ATM beams experienced significant scattering within the region of rough terrain, thereby causing large discrepancies in bias.

Modeled and Observed Bias
The relationship between optical grain size and ICESat-2 bias is a function of the path length of signal photons incident upon a snow surface.When photons interact with snow, there are two potential outcomes: reflection after subsurface scattering or absorption by the snowpack.The first outcome is more frequent for ICESat-2 and the ATM 532 nm beam over coarse-grained snow, and it is responsible for the largest biases over our study area.The second outcome occurs with the ATM 1064 nm beam when grain size is large, due to the low reflectance of aged snow in this part of the electromagnetic spectrum.Snow absorption reduces the occurrence of multiple scattering in NIR signals, so bias will be low over snow surfaces.The presence of snow impurities, such as black carbon and ice algae, may increase the probability of absorption for 532 nm signals and reduce bias (Smith et al., 2018), but further research is needed to confirm this hypothesis.
The waveform-based grain size estimates are generated based on a nominal surface density of 400 kg m 3 .This density is used to calculate (a) the velocity of light in the snow and (b) the distance that a photon must travel between scattering events for a given snow grain size, both of which are used to estimate the mean time required for the photon to travel between snow grains.As described in Section 2.3.2, a range of densities is given to the algorithm to generate an uncertainty that is a factor of two smaller or larger than the recovered grain size.If independent information were available about the surface density, this uncertainty could be reduced, but no such measurements were carried out during the IceBridge or AVIRIS-NG campaigns.
Model-derived results of altimetry bias have a strong dependence on optical grain size and density.As seen in Figure 4, the dependence on snow density is minor at small grain sizes, particularly at snow densities expected over Greenland.The modeled bias asymptotically approaches zero when grain size is small, which implies that the bias between two altimeters should be negligible when excluding other factors.Although discrepancies in bias become more evident between snow densities at larger grain sizes, we note that lower densities (i.e., ρ s = 100-200 kg m 3 ) only occur for fresh snow.The case studies presented here take place at the end of the Greenland melt season, where snow densities of ρ s = 300-500 kg m 3 and larger grain sizes are more common, as supported by the observations of Schaller et al. (2016) and Fausto et al. (2018).However, the model results also suggest that the optical grain sizes and densities needed to fully reproduce observed bias are implausibly large, so other factors, such as surface roughness, must be examined to characterize ICESat-2 biases.

Sources of Uncertainty
In CS1, we found that optical grain size retrievals differed substantially in magnitude between AVIRIS-NG and ATM.The exact cause may be related to the respective retrieval methods.Grain sizes from AVIRIS-NG are derived from the near-infrared surface reflectance of a location, so they are more sensitive to changes in surface snow properties.The ATM algorithm estimates grain size through received waveform pulses, which contain backscatter from below the snow surface if volumetric scattering is significant.The differences between the two algorithms are evident in Figure 5, which can be separated into fine-grained and coarse-grained regions, which correspond to sections 0-17 km and 17-42 km along the CS1 transect, respectively.In the fine-grained region, the pore spacing between snow grains is small, so ATM beam penetration beyond the surface layer will be minimal, and the derived optical grain size will be smaller as a consequence.Volumetric scattering becomes more significant at larger optical grain sizes, as is reflected in the coarse-grained region of Figure 5.If the ATM beam penetrates into the subsurface, grain sizes may be larger than those observed by AVIRIS-NG, assuming that grain size increases with depth.Although ICESat-2 was not considered in CS1, the strong agreement between ICESat-2 and modeled bias in Figure 10 indicates that retrieved ATL03 photon path lengths are more sensitive to subsurface grain sizes over aged or melting snow.
The data sets used here each have different approaches to estimating surface height, which may influence the biases given in Figures 7 and 10.Both ATM beams estimate surface height from the centroid of received waveforms, and a signal strength threshold is applied to filter noise.Although most background noise is removed with the threshold, sufficiently coarse snow or complex topography (i.e., crevasses or large slopes) may broaden waveforms and shift the centroid by nanoseconds, or centimeters in height change.ICESat-2 and the Monte Carlo model use similar approaches by estimating the median surface height (ICESat-2) or travel time (model) for photon aggregates.Subsurface scattering increases the distribution of photon delay times, therefore also increasing uncertainty and bias in ATL03 and the model.Thus, the differences between ATM and the model may be partly explained by these different approaches in signal processing.The two bias estimates show better agreement in the coarse-grained region of Figure 7, which may indicate that the model is neglecting snowpack features that impact the ATM signal, such as rough snow surfaces.A sufficiently rough snow surface may trap photons at either 532 nm or 1,064 nm, and photons returning to the detector would experience time delays comparable to coarse-grained snowpacks (Bair et al., 2022;Larue et al., 2020;Warren et al., 1998).A similar situation could occur over ice cracks or rough terrain, as suggested by the notable discrepancies in bias between ICESat-2, ATM, and the model in Figure 10.A model that allows for more complex scenarios, such as a rough surface layer or layer-dependent optical grain sizes, could help to answer these questions, though we leave the development of such a model to a future study.
The results of this study assume that the ATM NIR beam has negligible bias over the study location.Although we anticipate negligible bias caused by volumetric scattering, we also note that the NIR beam required calibration prior to use (Section 3.2).We determined that the NIR beam was biased relative to the green beam by 3.85 cm over a dark, flat surface, and Figure 3b shows that there is an uncertainty of ˜2 cm in the median bias.The spread in the NIR beam is likely caused by noise in the transmitted NIR waveforms, with the exact reason for the noise currently under investigation by the ATM team.
Additional uncertainties in ICESat-2 acquisitions may be categorized into time-variant and time-invariant uncertainty (Neumann et al., 2020).Both have been determined to be nearly constant over ICESat-2 tracks of 100 km or less.The time-variant biases are about 12 cm, and are caused by random errors from noise photons.These biases have been accounted for in current releases of ICESat-2 data products.The time-invariant biases are smaller at 1-2 cm, and they are likely caused by beam pointing errors inherent to all ICESat-2 data products.These pointing errors are expected to be corrected in a future data release.
As noted above, the ATL06 algorithm aggregates ATL03 photons, applies several corrections, and produces 40 m segment heights that are 3 cm lower than metrics based solely on ATL03 photon heights.Consequently, ATL06 biases in the presence of volumetric scattering should be 3 cm larger than the ATL03-based biases reported here.The impact of these biases should be greatest when comparing times of year with relatively little and relatively significant volumetric scattering, for example, summer versus winter surface heights.However, given the magnitude of seasonal elevation change in the ablation zone of Greenland, it may be difficult to isolate the magnitude of volumetric-scattering-based biases from the height change due to seasonal melt and accumulation.

Implications for Snow Studies
This study was performed to assess ICESat-2 measurements over snow.There has been increasing interest in using ICESat-2 to derive spaceborne measurements of snow depth (Bormann et al., 2018).Currently, DEMs from lidar are commonly used to estimate snow depth (Deems et al., 2013), though current lidar applications are restricted to airborne and ground-based surveys.There is a critical need to measure deep snow in forests and mountains using spaceborne instrumentation (Bormann et al., 2018;National Academies of Sciences, Engineering, and Medicine, 2018), and progress is this area is an objective of the NASA-sponsored Snow Experiment (SnowEx).ICESat-2 has the potential to support SnowEx objectives through interseasonal measurements of terrain height.
Field campaigns conducted for the SnowEx mission require measurements from the mid-latitude melting season, when the optical grain size of snow is largest.The results in Figure 10 indicate that ICESat-2 measurements over melting snow should be accurate to within ∼10 cm, assuming that the snow has compacted prior to melt.In contrast, ATM shows low bias for CS1, where the grain size is smaller.Although we were unable to consider ICESat-2 for CS1, the reasonable agreement between ICESat-2 and ATM (Figure 10) implies that ICESat-2 would experience minimal bias over locations with fresh, fine-grained snow.Hence, utilizing ICESat-2 for accurate measurements of snow depth is possible, though melting or aged snow may introduce bias and uncertainty in uncorrected ATL03 measurements.Higher-level products, such as ATL06 or ATL08, may reduce noise from the ATL03 photon cloud, but they will retain biased snow surface heights if subsurface scattering is unaccounted for, particularly over shallow snow.ICESat-2 tracks also rarely repeat in the mid-latitudes, so an accurate snowfree DEM would be required (Deschamps-Berger et al., 2023;Enderlin et al., 2022).

Conclusions
In this study, we used altimetry data from ICESat-2 and ATM to estimate volumetric scattering bias in snow.A combination of airborne optical grain size retrievals and Monte Carlo modeling was used to predict altimetry bias over the western Greenland ablation zone at the end of the melt season.ICESat-2 and the green ATM beam were compared to the near-infrared ATM beam to estimate observed bias.Our results suggest a positive relationship between the optical grain size of snow and altimetry bias over two case studies.The modeled results show that snowpack density is an important driver for volumetric scattering, though actual biases in the study locations remained consistent with densities of ∼315-500 kg m 3 .Although bias in both altimeters was generally within 10 cm, we cannot rule out more significant biases near the peak of the Northern Hemisphere melting season, when snow grain coarsening will enhance volumetric scattering at all snow densities.
The results in CS1 indicate that retrieved snow grain size is dependent on the instrument used.Grain sizes from AVIRIS-NG are based on the 1,030 nm absorption feature, thereby weighting snow grains near the surface, whereas ATM retrieves larger grain sizes within the snow subsurface.When combined with the Monte Carlo model, both data sets adequately reproduce trends in volumetric scattering bias, though the magnitude of the observed bias is not fully captured.Although used here, both instruments have infrequent coverage over midlatitude field sites, and ATM is not expected to collect data in the near future.Other sources of effective grain size, such as MODIS or Sentinel-3 (Mei et al., 2021;Painter et al., 2009), will therefore be needed for future volumetric scattering assessments over snow.
Further research is needed to identify altimetry biases in the presence of snowpack impurities or rough topography.The full impact of dust and black carbon on altimetry signals is not known, so there is a need for accurate airborne and satellite retrievals of surface impurity content.Similarly, a correction for rough or sloped terrain is needed, given that both factors increase height uncertainties for both ICESat-2 and ATM.To address these problems among others, a follow-up study is in preparation that validates the results in this paper over midlatitude snow.The SnowEx mission is conducting airborne lidar surveys for its 2023 Alaska campaign, several of which are expected to have significant overlap with ICESat-2 tracks.The flights will overpass coastal and forested regions of Alaska, so we anticipate a more rigorous analysis of ICESat-2 over multiple terrain types.The expected result is a bias correction algorithm that ideally will be applicable to all snow surfaces.

Figure 1 .
Figure1.A demonstration of volumetric scattering for two lidar wavelengths (532/1,064 nm) in a semi-infinite snowpack simulated using a Monte Carlo model.The snow is assumed to be clean (i.e., no impurities) with ρ s = 400 kg m 3 and r eff = 50 μm (top), r eff = 500 μm (bottom).The histograms were generated using 10 5 photon path simulations.

Figure 2 .
Figure 2. (a) ESRI imagery showing the location of the two case studies over the Greenland Ice Sheet, with a reference relative to the entire ice sheet given in the upper right.(b) Zoomed perspective of the first case study (CS1), using Landsat-8 imagery from 31 August 2019.The red line is the path flown by both Airborne Topographic Mapper (ATM) and AVIRIS-NG, and the false color overlay is the snow grain sizes derived by AVIRIS-NG.(c) Same as (b) but for the second case study (CS2), using Landsat-8 imagery from 2 September 2019.The red line here represents the region where ATM and the central strong ground track of ICESat-2 (GT2L) intersect.

Figure 3 .
Figure 3. (a) Location of the aircraft ramp used to assess the bias between the Airborne Topographic Mapper (ATM) 532 and 1,064 nm beams.The black box highlights the ramp overpass.(b) Along-track scatter plot of the 532 nm bias relative to 1,064 nm measurements.Negative values indicate lower surface heights measured by the 1,064 nm beam.The red dashed line depicts the median bias observed across the overpass.(c) Bias distribution between the ATM beams, as derived from ramp overpass measurements.

Figure 4 .
Figure 4. Modeled altimetry bias derived using median path lengths estimated from a Monte Carlo model(Schneider et al., 2019).Bias is given as a function of snow optical grain size and snowpack density.The snow density ρ s = 315 kg m 3 is used to represent the average snow density reported over Greenland byFausto et al. (2018) (ρ s = 315 kg m 3 ).The simulated bias at ρ s = 300 kg m 3 is nearly identical to that of ρ s = 315 kg m 3 and was thus omitted.

Figure 5 .
Figure 5. Snow optical grain sizes derived along-track from AVIRIS-NG (light purple) and Airborne Topographic Mapper (ATM) waveforms (orange) for CS1.The shading for both curves represents uncertainty in grain size.The AVIRIS-NG uncertainty is an assumed constant of 50 μm, and the ATM uncertainty is derived from the retrieval algorithm assuming ρ s = 200 kg m 3 (lower) and 800 kg m 3 (upper).The dashed black line is the difference in grain size between the two data sets (AVIRIS-ATM).

Figure 7 .
Figure7.Observed Airborne Topographic Mapper (ATM) green-NIR range differences (green) as compared to modeled estimates using optical grain sizes from AVIRIS-NG reflectances (purple) and ATM waveform fitting (orange).The solid lines for the modeled estimates represents ρ s = 400 kg m 3 , whereas the shading is the uncertainty due to realistic changes in snow density, given ρ s = 315-500 kg m 3 .

Figure 8 .
Figure 8. Snow optical grain sizes for CS2, as derived from 532 nm Airborne Topographic Mapper waveforms.As with Figure 5, the shading represents the uncertainty due to snow density, given ρ s,min = 200 kg m 3 and ρ s,max = 800 kg m 3 .

Figure 9 .
Figure 9. (a) Along-track surface heights from ATL03 (light blue), the Airborne Topographic Mapper (ATM) 532 nm beam (green), and ATL06 (black).The plot only includes spot measurements where ATM was within 12 m of ICESat-2 footprints.The black box highlights a surface feature where ATM and ICESat-2 show greater disagreement.(b) IceBridge Continuous Airborne Mapping By Optical Translator imagery centered on the boxed feature in panel (a).The colored line is the ICESat-2 biases relative to ATM 1064 nm.

Earth and Space Science
FAIR ET AL.