Mapping Sea Ice Surface Topography in High Fidelity With ICESat‐2

The Advanced Topographic Laser Altimeter System on Ice, Cloud and land Elevation Satellite 2 (ICESat‐2) offers a new remote sensing capability to measure complex sea ice surface topography. We demonstrate the retrieval of six sea ice parameters from ICESat‐2/Advanced Topographic Laser Altimeter System data: surface roughness, ridge height, ridge frequency, melt pond depth, floe size distribution, and lead frequency. Our results establish that these properties can be observed in high fidelity, across broad geographic regions and ice conditions. We resolve features as narrow as 7 m and achieve a vertical height precision of 0.01 m, representing a significant advance in resolution over previous satellite altimeters. ICESat‐2 employs a year‐round observation strategy spanning all seasons, across both the Arctic and Southern Oceans. Because of its higher resolution, coupled with the spatial and temporal extent of data acquisition, ICESat‐2 observations may be used to investigate time‐varying, dynamic, and thermodynamic sea ice processes.


Introduction
Observational evidence from multiple data sources demonstrates that significant, and rapid, changes are occurring in the Arctic climate system . As air temperatures in the Arctic warm at twice the global rate, long-term declines in sea ice extent, age and volume, and the duration of the winter growth period have been observed (Perovich et al., 2019). Arctic sea ice influences global atmospheric patterns (e.g., Francis et al., 2017) and oceanic thermohaline circulation and, due to its high albedo, regulates the planetary energy balance (e.g., Curry et al., 1995). Sea ice properties (e.g., concentration, thickness, and drift velocity) and processes (e.g., growth, melt, divergence, and convergence) are, however, some of the most poorly constrained variables in global climate models; neither their magnitude nor their impact on future climate projections are well understood (e.g., Turner & Comiso, 2017). High-resolution satellite measurements offer a practical solution for achieving spatially and temporally complete monitoring of sea ice in the polar oceans (Shepherd et al., 2018).
Satellite altimeters, deployed on ICESat (2003)(2004)(2005)(2006)(2007)(2008)(2009) and CryoSat-2 (2010 up to the present), with orbital inclinations designed to observe Earth's polar regions, have delivered near-continual winter-time measurements of sea ice topography at the basin scale since 2003 (e.g., Laxon et al., 2013). These observations have revealed a decline in Arctic sea ice freeboard, thickness, and volume over the last two decades (e.g., Farrell ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2020GL090708
Special Section: The Ice, Cloud and Land Elevation Satellite-2 (ICESat-2) on-orbit Performance, Data Discoveries and Early Science  Laxon et al., 2013), during which time the ice cover has transitioned from predominantly multiyear to seasonal ice (Perovich et al., 2019).
NASA's ICESat-2 continues the polar satellite altimetry record, measuring sea ice elevation to 88°N/S. Since 14 October 2018, ICESat-2 has provided continual observations of both polar regions, with the exception of the 2019 summer melt season when an observational gap occurred between 26 June and 26 July 2019 due to a spacecraft anomaly. The data have been used to track the evolution of sea ice freeboard in winter (e.g., . Here we present a selection of high-fidelity measurements of sea ice surface topography from ICESat-2, for a variety of Arctic sites. Our goal is to demonstrate ICESat-2's unique capability to track individual floes, from which fine-scale sea ice properties may be derived. We show examples spanning the end of winter (April 2019), through summer melt (June 2019), and fall freeze-up (September 2019).
Results are validated using independent, coincident observations from airborne lidar and satellite imagery. We discuss how ICESat-2's remote sensing capabilities over sea ice will extend the utility of the data beyond fulfillment of the mission science requirement to measure freeboard (Markus et al., 2017), by enabling more detailed process studies.

Data
ICESat-2 operates in a 91-day exact repeat orbit, with 1,387 orbits per cycle. Over the Arctic Ocean, orbit subcycles of 4 and 29 days offer complete, basin-scale coverage. ICESat-2 carries one primary instrument, the Advanced Topographic Laser Altimeter System (ATLAS). Six ATLAS beams, arranged in three pairs, span approximately 6.6 km in the across-track direction. Beam locations are defined relative to spacecraft Reference Ground Tracks (RGTs). We use the convention ttttccss_gtxy to identify specific RGTs, where t is the RGT number, c is the orbit cycle, s is segment number, gt indicates "ground track," x is beam number, and y indicates either left (l) or right (r) beam. Controlled pointing to the RGTs began in March 2019.
Here we use the ATLAS Release 003 Level 2 ATL03 product that contains geolocated photon heights relative to the WGS84 reference ellipsoid (Neumann et al., 2020). Geolocation of individual photons results in a vertical range accuracy of 0.05 m and a precision better than 0.13 m . We also use Release 002 Level 3 ATL07 sea ice heights (Kwok, Cunningham, et al., 2019), derived from ATL03.
ICESat-2 retrievals are validated using two independent data sets. Dedicated Operation IceBridge (OIB) under-flights of ICESat-2 were conducted in April 2019 to obtain high-resolution (2 m) Airborne Topographic Mapper (ATM) lidar data along ICESat-2 RGTs in the Canada Basin .
Here we present results from 22 April 2019, when near-coincident (less than~38 min) ATM data were acquired. We also utilize high-resolution (10 m) visible imagery from the Sentinel-2 MultiSpectral Instrument (MSI) for validation.

Fine-Scale Sea Ice Properties
ICESat-2's predecessor ICESat carried an analog laser altimeter that had a large footprint (~50 m) with 172 m spacing between shots, which limited the resolution of sea ice observations (Farrell et al., 2011). Overlapping~12-m diameter ICESat-2/ATLAS footprints (L. Magruder, personal communication, 6/8/ 2020), sampled every~0.7 m along-track, offer a unique opportunity for adaptive sampling of the surface, at length scales suitable for discriminating discrete sea ice features. Following previous work using an airborne simulator for ATLAS (Farrell et al., 2015), we exploit the innovation of photon-counting laser altimetry to map the rough, topographically complex, sea ice surface at high resolution. Prior to applying basic sea ice retracking algorithms to ATL03 photon heights, we preprocess the data to remove background noise. This is achieved by retaining only photons between the 15th and 85th percentile of the per-shot height distribution. ATL03 atmospheric, tidal, and geoid corrections are applied to obtain corrected elevation. Following Duncan et al. (2018) and Farrell et al. (2015), elevation is the height anomaly relative to the local level ice/water surface.
In the following subsections, we explore ICESat-2's capabilities to observe signatures of both sea ice dynamics (ice pack convergence and divergence) and thermodynamics (sea ice melt/freeze). We include a brief description of the retrieval of six sea ice parameters from ATL03 data: surface roughness, ridge height, ridge frequency, melt pond depth, floe size distribution, and lead frequency. We then discuss (in section 4) their utility in sea ice process studies.

Surface Roughness and Pressure Ridges
Sea ice roughness (σ h ) provides an indication of both the mechanical deformation history of the ice cover and snow distribution across the surface. It is also a proxy for ice thickness. Knowledge of σ h is required to understand the exchange of turbulent energy between the ice and atmosphere and drag-induced ice dynamics (Zwally et al., 2002). Here, σ h is the standard deviation of ATL07 sea ice height within 25-km-long segments. To illustrate ICESat-2's capability for obtaining σ h , we consider the Arctic Ocean as a whole but also highlight two regions of the ice cover (A, spanning 79.5-83°N, 100-120°W, and B, spanning 76.25-81°N, 140-160°W) with distinct roughness characteristics ( Figure 1a). In April 2019, Arctic-wide σ h averaged 0.18 m (Figure 1b) but showed a spatial pattern ( Figure 1a) consistent with the known geographic locations of the seasonal and multiyear ice zones. Perovich et al. (2019) show that Region A, north of Borden Island, contained multiyear ice ≥3 years old, while Region B was an area of seasonal ice within the Beaufort Gyre. ICESat-2 shows the seasonal ice in Region B was half as rough (σ h = 0.15 m, Figure 1b) as the multiyear ice in Region A (σ h = 0.30 m; Figure 1b).
Focusing on representative, 1-km-long segments within each region, we apply the University of Maryland-Ridge Detection Algorithm (UMD-RDA) to the preprocessed (section 3) ATL03 photon heights. The UMD-RDA retains the 99th percentile of the photon height distribution for a 5-shot aggregate, applied on a per-shot basis so as to retain full along-track resolution (0.7 m). When applied to ICESat-2 retrievals over sea ice, we obtain an elevation profile from which individual pressure ridges may be detected. Following previous studies (e.g., Duncan et al., 2018), we define a pressure ridge sail as any local maxima occurring 0.6 m above the local level ice surface. This threshold distinguishes ridges from lower-amplitude surface features (e.g., snow dunes or sastrugi). Local minima lying above the threshold height are also flagged. Ridge width is the along-track distance between minima or the point(/s) at which elevation drops below the threshold height, whichever is closer to the local maxima. Maxima separated by ≤10 m are not considered unique and are instead counted as a single ridge (e.g., Ridge 8; Figure 1c).
Results from the UMD-RDA reveal an average sail height (h s ) of 1.5 m for nine ridges spanning 10.7-51.8 m in width for a 1-km-long segment in Region A (Figure 1c). Coincident OIB ATM elevations show h s averaged 1.6 m, verifying the ICESat-2 UMD-RDA result. The altimeter height comparison, as well as insight from OIB imagery (Figure 1d), confirms both the location and number (n r ) of distinct ridges. The ATL07 algorithm accurately detects individual deformation features in this region of multiyear ice; however, h s is underestimated by 0.4 m (0.3 m) when compared with ATM (UMD-RDA) (Figure 1c). The results from Region A contrast with the UMD-RDA statistics obtained over the smoother surface of Region B. Assessing a 1-kmlong segment in B, h s averaged 0.8 m for n r = 5, ranging 7.1-35.7 m in width (Figure 1e). In this example the ATL07 data set performs poorly, enabling the detection of only one ridge, with h s = 0.61 m, suggesting surface roughness on seasonal ice may be underestimated by the ATL07 algorithm.
Regions A and B contain approximately the same number of 1-km segments (n seg ; Figure 1f ), derived from between 64 and 78 RGTs in April 2019. Aggregating these measurements illustrates distinctions in the number of ridges (n r ) and their frequency (f r ; Figure 1f ) and in h s (Figure 1g), as a function of ice type. The rougher, older ice in Region A was more heavily deformed than the ice in Region B, with~2.5 times more ridges that were, on average, 0.28 m (0.25 m) taller in mean (modal) h s . We found that the 99th percentile of sail height (h s 99 ) was 0.67 m larger in Region A than in Region B. These results are consistent with an earlier study ) that found h s 99 is a strong indicator of the predominant ice type in which a pressure ridge forms.

Melt Ponds
Following the end of winter, as air temperatures warm, both thermodynamic and dynamic processes introduce meltwater to the system. The presence of low-albedo ponds on the sea ice surface increases the ice melting rate (Perovich & Polashenski, 2012) and enhances the ice-albedo feedback by increasing absorption of shortwave radiation, altering Earth's energy budget (Curry et al., 1995). The detection of melt ponds with spaceborne sensors has proved challenging since ponds are radiometrically similar to open water/leads and cover small areas (~5-100 m 2 ; Perovich et al., 2002). Early ICESat-2 observations demonstrated an unexpected capability to penetrate shallow, low turbidity water to measure coastal bathymetry and identify glacial melt ponds on Antarctic ice shelves Parrish et al., 2019). These early results, coupled with ICESat-2's high resolution, suggest the possibility of measuring sea ice melt pond depth and motivate the following analysis.
We examine ten 1-km-long ATL03 segments (Figure 2a Following Buckley et al. (2020), we classified open water, ponded surfaces, and ice floes in the Sentinel-2 scene ( Figure 2b) and used this to validate the presence of ponds in the ICESat-2 data. By tracking the movement of 10 floes between two overlapping Sentinel-2 images acquired 50 min apart (not shown), we estimated an average ice drift rate of 9.3 cm s −1 . To account for the time elapsed between the Sentinel-2 and ICESat-2 acquisitions, we applied a drift correction of 206 m to the imagery. This provided the exact geolocation of melt features in the Sentinel-2 scene at the time of the ICESat-2 overpass. Assessment of the ICESat-2 segments (Figure 2a) reveals strong returns from the approximately level sea ice surface and classic concave pond features (Perovich et al., 2003). The latter are a result of ICESat-2 returns from melt pond bottoms (MP1-10; Figure 2a). Four ponds (MP7-10) can be identified in both the Sentinel-2 (Figure 2b, insets) and ICESat-2 data (Figure 2a).
The small-scale pond features, ranging~60-280 m wide (Figure 2a), are not captured by higher-level ICESat-2 products, such as ATL07 (Figure 2a, cyan). Hence, we developed the University of Maryland-Melt Pond Algorithm to identify pond surfaces (Figure 2a, black) and bathymetry (Figure 2a, magenta) in the ATL03 data. To determine pond depth (h mp ), the algorithm utilizes a two-dimensional histogram with 10-m alongtrack and 0.1-m vertical resolution. Ponds occur where a bimodal distribution exists, with a primary mode close to the mean segment elevation indicating the ponded water surface and a secondary mode below the water surface (e.g., see MP2; Figure 2a). The leading edge of the secondary mode defines pond bathymetry since this represents the first photon returns from the pond bottom. Its selection mitigates the impact of photons with delayed arrival times at the detector. Initially, h mp is derived by subtracting bottom elevations from pond surface elevation. We note, however, that photon heights for photons returned from within ponds are not inherently corrected for the refraction of light at the air-water interface. We therefore apply a refraction correction (following Parrish et al., 2019), wherein estimated pond depth is scaled by 0.749, the ratio of the refractive index of air (1.00029) to water (1.33567). After correcting for refraction and linearly interpolating at 5-m along-track resolution, the h mp distribution for MP1-10 indicates that h mp ranged 0.04-2.4 m, with a modal (mean) depth of 0.35 m (0.80 m) (Figure 2c). Seventy-five percent of the h mp retrievals were ≤1.1 m. MP9 (inset, Figure 2b), an approximately circular pond, is likely younger than the other geometrically more complex ponds (Perovich et al., 2002). Combining maximum h mp (1.73 m; Figure 2a) with Sentinel-2 pond area (37,000 m 2 ; Figure 2b) and assuming pond volume is approximated by the volume of a spherical cap, we estimate that MP9 contains~32,000 m 3 of melt water.
The ICESat-2-derived estimates of maximum h mp are deeper than those typically observed in the field (e.g., Perovich et al., 2003). Pond depths reported here can, however, be explained by their geographical setting on rough multiyear ice (Skyllingstad et al., 2009) and the atmospheric conditions under which the ponds formed. Regional temperatures in May 2019 were~4-6°C above average (Vose et al., 2014) allowing for enhanced snow melt and mature pond evolution. Sophisticated simulations of pond evolution (Scott & Feltham, 2010) have suggested that rapid pond deepening, of over 0.5 m in 10 days, can occur on thick, rough multiyear ice, with mean pond depth reaching 0.85 m, in line with the observations shown here.

Floe Size Distribution and Lead Frequency
As the melt season progresses, mechanical breakup continues, and the unconsolidated ice pack comprises discrete floes in free drift. Open water fraction increases rapidly, amplified by lateral melt, and floe size decreases. Solar heat input to the upper ocean increases, further enhancing melt . Lead and floe size statistics, and their temporal and regional variability, are needed to understand ice-ocean-atmosphere heat fluxes in summer. Previously, satellite altimeters faced challenges observing summer ice processes. ICESat operated in campaign mode and did not obtain summer data (Farrell et al., 2009), while CryoSat-2 has limited along-track resolution (~300 m) and cannot distinguish between summer melt features (ponds and leads) on the basis of radar altimeter return power, since the radar backscatter coefficient of sea ice is sensitive to meltwater (Wingham et al., 2006).
Here, we revisit the ice cover at the onset of fall freeze-up, focusing on ICESat-2 observations in the Canada Basin (Region D; Figure 1a) in early September. Floes in the region range from tens to ten thousands of meters wide (Figure 3a) and are surrounded by thin nilas (WMO, 1970). Applying the UMD-RDA to a 200-km-long ICESat-2 transect, we obtain elevation profiles for the three strong beams. Aggregating heights from 10 short (~1 km-long) lead segments crossing smooth nilas (cyan diamonds, Figure 3a), we obtain a zero-mean elevation distribution with a standard deviation of 0.009 m (blue curve, Figure 3b), illustrating the vertical height precision of ATLAS over leads. This is a 50% improvement in capability compared with ICESat, which had a demonstrated precision of~0.02 m over leads (Kwok et al., 2004). We note that while the major mode of the ATL07 height distribution for the same leads is consistent with the UMD-RDA results (gray curve, Figure 3b), 23% of the data fall into a secondary mode, with a mean elevation of 0.05 m. The reason for this secondary mode is currently unknown, but its impact is a positive bias in ATL07 sea surface heights.
To demonstrate ICESat-2's ability to discriminate individual floes, we examine a shorter (~7.5 km) representative area. In this region, we compare ICESat-2 elevations with a coincident Sentinel-2 image acquired just 11 min after the ICESat-2 pass (Figure 3c). The imagery reveals nilas between floes of varying sizes, with some evidence of finger rafting. The ICESat-2 lead/floe locations, and deviations in their elevation, accurately correspond with the local ice conditions revealed in the Sentinel-2 MSI data. The ICESat-2 retrievals demonstrate floes ranging 20 m to 3.034 km wide. ICESat-2 modal freeboard, computed at the floe scale, averaged 0.44 m and ranged 0.05 to 1.35 m, for the floes shown in Figure 3c. If we suppose that the floes are in hydrostatic equilibrium with an ice density of 880 kg m −3 and that a thin (~0.05 m), low density (~220 kgm −3 ) dusting of snow has accumulated on these floes, we can estimate an average ice thickness of 2.85 m ± 2.19 m, which is reasonable when compared with ice mass balance estimates .
We extend the analysis to~600 km by combining ICESat-2 retrievals from the three strong beams and compute lead and floe statistics. Because of the orientation of the track relative to the floes (Figure 3a), we do not strictly measure lead width or floe diameter. But due to a >600-km sample size, the statistics are regionally and seasonally representative. Classifying open water and ice floes in the Sentinel-2 scene (following Buckley et al., 2020), we tagged lead and floe pixels along the ICESat-2 track and used these for validation. Based on the results in Figure 3b, leads are identified in the ICESat-2 data as level ice surfaces with ≥15 contiguous retrievals (~10-m along-track width) within 0.1 m of local reference surface with a standard deviation of ≤0.01 m. Lead retrievals accounted for 27.6% of the ICESat-2 data, which is consistent with a regional open water fraction of 25.1% derived from Sentinel-2. ICESat-2 retrievals indicate zero to two distinct leads per kilometer, with an average lead frequency of 1.3 km −1 , in close agreement with Sentinel-2 ( Figure 3d). While leads ranged 10 m to >3 km, average (median) lead width was 235 m (71 m), and 75% of leads were <200 m wide ( Figure 3e). Floes, on the other hand, averaged 479 m, and 75% were <600 m wide (Figure 3f ). Floe and lead widths differed by 0-14 m between the two independent estimates (Figures 3e and  3f), demonstrating the quality of the altimeter-derived metrics. The ICESat-2 statistics are consistent with recent studies wherein high-resolution optical and synthetic aperture radar (SAR) imagery revealed summer ice floe diameters ranging tens to thousands of meters, averaging <200 m wide (Arntsen et al., 2015;Hwang et al., 2017).

Discussion
The evidence provided here establishes that the small footprint and high-pulse repetition frequency (10 kHz) of ATLAS on ICESat-2 are capable of resolving individual floes, sails, and ponds on the sea ice cover. ICESat-2 retrievals deliver unprecedented quality in the measurement of sea ice properties including floe size distribution, lead frequency, and sail height and frequency. They also offer the new capability to measure sea ice melt pond depth, derived from pond bathymetry, the first such measurements to be retrieved using spaceborne altimetry. We have shown that sea ice features as narrow as 7.1 m may be detected. ICESat-2's capability to retrieve fine-scale sea ice properties year-round will be transformational in understanding time-varying sea ice processes and will advance interpretation of lower-resolution remote sensing data. Evaluating the skill of sea ice process models has been heretofore hindered by a lack of high-resolution observations covering large spatial and temporal scales (Roach et al., 2018). The ICESat-2 observation strategy will also address this need.
By applying customized surface retracking algorithms to ATL03 photon heights, we captured signatures of dynamic ice convergence (pressure ridges), divergence (leads), and thermodynamic melt. Our examples show that over level surfaces (recently refrozen leads), an elevation precision of 0.01 m can be achieved, representing a considerable advance over ICESat. Discrimination of leads, and accurate measurement of their height, is critical for retrieval of sea ice freeboard and thickness using altimeter techniques, because lead elevation approximates local sea surface height and thus provides a reference level from which freeboard can be derived (Farrell et al., 2009). Here we have demonstrated that floe-scale freeboard may be retrieved with ICESat-2. Statistical analysis of individual floes, and their freeboard, will inform algorithm development for current and future satellite altimeter missions. The data may be used to study the joint sea ice floe size thickness distribution, which is required to understand the impact of the geometrical sampling error in lower-resolution radar altimeter retrievals (Envisat, CryoSat-2) over sea ice (Wingham et al., 2006). Quantifying this is critical for combined sea ice thickness estimates from multiple altimeters with varying resolutions.
High-resolution observations from ICESat-2, such as those shown here, will fill a gap in knowledge required to advance sea ice modeling (Horvat & Tziperman, 2017). For example, ICESat-2 measurements of sea ice growth in winter compliment those obtained by CryoSat-2 and will enable investigations of both dynamic and thermodynamic thickening. Routine retrieval of surface roughness, lead frequency, and the floe size thickness distribution will help advance drag parameterization in sea ice models. Observing sea ice evolution during summer melt and fall freeze-up will also improve our understanding of thermodynamically driven mass loss . Extending the present analysis of melt ponds to determine regional variability in pond depth and volume, and their temporal evolution, will be useful for assessing current parameterizations in melt pond models (Hunke et al., 2013). Though Sentinel-2 coverage across the Arctic is limited (by orbit inclination, cloud interference, and crossover timing with ICESat-2) and restricts the extension of the present approach to estimate pond volume at basin scale, the potential exists to estimate surface topography, melt pond depth, and ice thickness simultaneously with ICESat-2 data alone. This may be helpful in further constraining our understanding of ice albedo evolution during summer (Eicken et al., 2004). Tracking both dynamic and thermodynamic processes with ICESat-2, across a more comprehensive range of sea ice conditions, and seasons, than was heretofore possible with remote sensing techniques, will also permit examination of ice-ocean-atmosphere exchanges of energy, mass, and momentum, supporting improved understanding of the connections between sea ice variability and climate forcings.