Gaussian Process Interpolation Turns Sparse ICESat-2 Data Into Detailed Glacier Change Maps

A new study demonstrates that a statistical technique borrowed from machine learning can transform the sparse, discrete measurements of NASA’s ICESat-2 satellite laser altimeter into spatially continuous maps of glacier surface height change, even across some of the most topographically demanding ice terrain on Earth. Thorsten Seehaus of Friedrich-Alexander-Universität Erlangen-Nürnberg, together with Alex S. Gardner and Johan Nilsson of the Jet Propulsion Laboratory and Uppsala University, applied Gaussian Process Regression to ICESat-2 elevation data over the Larsen-B embayment on the Antarctic Peninsula and central Southern Svalbard. Their results, published in The Cryosphere, show that the method reproduces independent radar-based measurements of ice loss at Larsen-B within roughly 13 percent, while revealing surge activity on Svalbard glaciers earlier than established detection approaches.

The challenge the researchers set out to address is fundamental to satellite altimetry. Laser altimeters such as ICESat, which operated from 2003 to 2009, and its successor ICESat-2, launched in September 2018, deliver highly accurate height measurements but only at discrete points along narrow ground tracks. ICESat-2 carries the Advanced Topographic Laser Altimeter System, which fires 532-nanometer laser pulses at a repetition rate of 10 kilohertz, producing roughly 0.7 meters of along-track spacing between pulses in a six-beam configuration arranged as three pairs. Vertical precision reaches about 2 centimeters over flat, hard surfaces, and the satellite repeats its orbit every 91 days. Yet between the ground tracks lie vast unsampled areas, and cloud cover further thins the data, because laser altimetry requires clear skies. Radar altimeters such as ESA’s CryoSat-2 can measure through clouds but lose accuracy over steep slopes and narrow glacier valleys, where waveform interpretation and surface slope effects degrade performance. Converting the sparse laser sampling into continuous height-change fields therefore requires interpolation, which underpins glacier mass balance assessments and data assimilation.

Gaussian Process Regression is a non-parametric, Bayesian framework closely related to classical geostatistical kriging, but with distinct practical advantages. Both approaches represent the interpolated field as a weighted linear combination of observations derived from an assumed spatial covariance structure. However, classical kriging generally operates in two-dimensional geographic space, requires manual variogram fitting, and needs separate regression steps to incorporate explanatory variables such as elevation or slope. GPR instead treats spatial coordinates and auxiliary terrain or glacier-specific predictors jointly within a single covariance kernel, allowing complex terrain dependencies, anisotropy, and non-stationary behavior to be encoded directly in the kernel design. Kernel hyperparameters are optimized automatically by maximizing the marginal likelihood, avoiding subjective variogram fitting. Crucially, GPR yields a full posterior predictive distribution rather than a single point estimate, providing spatially distributed uncertainty estimates for every predicted location. Its main drawback is computational cost, which scales cubically with the number of training points, though sparse approximations mitigate this limitation.

The two study regions were chosen precisely because they stress-test any interpolation scheme. The Larsen-B embayment on the eastern Antarctic Peninsula has undergone dramatic change since the disintegration of the Larsen-B Ice Shelf in 2002, which triggered acceleration, dynamic thinning, and retreat of its tributary glaciers. Although flow speeds later declined, persistent sea ice appears to have temporarily stabilized the glacier fronts until its breakup in early 2022, after which several outlet glaciers renewed their retreat and acceleration. Rugged topography and steep, narrow glacier valleys challenge both radar and laser altimetry there. Central Southern Svalbard, by contrast, hosts a diverse mix of land-terminating and tidewater glaciers, including several surge-type glaciers that have been active in recent years. Three decades of atmospheric warming, increased surface melt, and widespread frontal retreat have driven increasingly negative mass balances across the archipelago, and moderately rugged terrain with narrow valleys adds further interpolation difficulty.

The team processed ICESat-2’s ATL11 product, which provides surface height time series at fixed reference points along repeat ground tracks, using all available data from March 2019 through December 2024. Height change rates were computed for multi-annual, annual, and seasonal intervals, with quality filters excluding points whose combined uncertainty exceeded 2 meters and removing extreme rates above plus or minus 30 meters per year for long-term datasets. An additional hypsometric filter, binning measurements in 50-meter elevation intervals and flagging values beyond three times the normalized median absolute deviation, dismissed only about 0.4 to 0.7 percent of measurements. To select the best predictors, the authors employed a permutation feature importance analysis, randomly shuffling candidate variables and measuring the resulting degradation of the output. Shuffling surface elevation and glacier identity produced strong noise, whereas shuffling ice velocity and slope had little effect, so elevation and glacier IDs were retained. The final model used a Matérn 5/2 kernel with Automatic Relevance Determination across four input dimensions, with correlation lengths of 10,000 meters for the horizontal coordinates, 100 meters for elevation, and 1 for glacier ID, all optimized by maximum likelihood.

At Larsen-B, the GPR-derived average height change rate for 2021 to 2024 was minus 0.61 plus or minus 0.02 meters per year, corresponding to a volume loss of 4.55 plus or minus 0.17 cubic kilometers per year, excluding sub-aqueous losses at potentially floating termini. The strongest thinning occurred along the marine-terminating glaciers that once fed the collapsed ice shelf, particularly Hektoria, Green, and Crane glaciers. Validation against an independent TanDEM-X radar dataset covering 95 percent of the study area, corrected for elevation-dependent signal penetration bias, yielded a mean rate of minus 0.72 plus or minus 0.01 meters per year and a volume change of minus 5.26 plus or minus 0.07 cubic kilometers per year, with a mean offset of just 0.09 meters per year between the two products. The annual time series revealed moderate lowering of roughly 3 to 5 meters per year at the Hektoria and Green glacier tongues before 2022, followed by rapid surface lowering of 15 to 36 meters per year at their termini in 2022 to 2024, consistent with renewed dynamic acceleration after the sea ice breakup. Crane glacier showed a delayed response, with intensified lowering of about 17 meters per year appearing only in 2023 to 2024; the authors hypothesize that its narrower fjord and small calving front, factors tied to bedrock geometry, explain the lag.

Over Svalbard, the method detected widespread thinning averaging minus 1.62 plus or minus 0.03 meters per year, equivalent to an ice volume change of minus 4.85 plus or minus 0.10 cubic kilometers per year. That rate is roughly twice the minus 0.80 and minus 0.84 meters per year reported by earlier studies of Southern Svalbard covering 2003 to 2019 and 2019 to 2021 respectively, and the associated uncertainty is about six times lower. The authors attribute the elevated loss partly to the fact that all glacier surges in Southern Svalbard identified in a recent survey fall within their study region, since surges relocate ice from higher to lower elevations where ablation rates are greater. The height-change maps clearly captured both down-glacier propagating surges, marked by surface lowering upstream and thickening toward the glacier tongue, and up-glacier propagating surges, initiated by drastic lowering at the terminus. Notably, the onset of surging at Scheelebreen and Vallåkrabreen could be identified in the 2019 to 2020 annual fields, earlier than the surge start dates derived from Sentinel-1 radar coherence analysis, and precursors of mass downwasting at Nordsysselbreen, Doktorbreen, and Paulabreen appeared years before their reported surge onsets after June 2024.

The validation exercises revealed both strengths and limits of the approach. Compared with the ICESat-2 ATL15 gridded product at 1-kilometer resolution, the GPR fields resolved small-scale features, including narrow dynamic glacier tongues and surge-related structures, that ATL15 missed entirely, with local offsets exceeding 10 meters per year in places. Region-wide, however, mean offsets were modest, 0.13 meters per year for Svalbard and 0.09 meters per year for Larsen-B, and ATL15 captured 94 percent of the GPR-derived volume loss at Svalbard but only 73 percent at Larsen-B, where spatial regularization dampened the strong thinning at the retreating glacier fronts. In artificial gap experiments, circular voids of increasing radius up to 7,000 meters were carved into the data and the withheld measurements used as truth. At Larsen-B, mean offsets stayed between minus 0.12 and plus 0.04 meters per year with standard deviations below 0.7 meters per year, while Svalbard showed growing error spread beyond 3,000-meter gap radii, reflecting its heterogeneous glacier dynamics. Interpolation errors increased with gap size, surface slope, the magnitude of the height-change signal, and distance to the nearest observation. Classical benchmarks, inverse distance weighting and nearest-neighbour interpolation, produced striping artifacts aligned with the satellite ground tracks and smeared signals in dynamically active zones; although their regional mean offsets were as small as 0.01 to 0.02 meters per year, local standard deviations of 0.87 to 1.16 meters per year exceeded those of GPR, confirming its superior local fidelity.

The study’s authors conclude that Gaussian Process Regression offers a robust, scalable, and statistically rigorous framework for future glacier monitoring, one that enhances the spatial and temporal resolution of altimetry-based observation while quantifying its own uncertainty. Because regular large-scale surface height change information has lagged behind the near-global availability of glacier velocity fields, the ability to generate area-wide, at least yearly elevation-change products could improve surge detection, process understanding, and early warning systems, with relevance for natural hazard management in regions such as High Mountain Asia, where glacier surges have caused flooding events. The authors note that meaningful results at seasonal resolution, roughly one 91-day ICESat-2 cycle, were hampered by persistent cloud cover at both test sites but may prove feasible in less cloudy regions, and that glaciers of different sizes, from small mountain glaciers to large ice caps, may require adjusted kernel parameters. As ICESat-2 data accumulate and additional altimetry missions launch, the approach promises refined mass balance assessments and sharper eyes on the localized, transient dynamics that reshape the world’s glaciers.

Subject of Research: Gaussian Process interpolation of ICESat-2 altimetry data to reconstruct glacier height change and surge dynamics in complex terrain

Article Title: Reconstructing Glacier Dynamics in Complex Terrain with ICESat-2 and Gaussian Process Interpolation

Article References: Seehaus, T., Gardner, A. S., & Nilsson, J. (2026). Reconstructing Glacier Dynamics in Complex Terrain with ICESat-2 and Gaussian Process Interpolation. The Cryosphere, 20(10), 6107-6130. https://doi.org/10.5194/tc-20-6107-2026

Image Credits: AI Generated

DOI: 10.5194/tc-20-6107-2026

Keywords: ICESat-2, glacier height change, Gaussian Process Regression, Larsen-B, Svalbard, glacier surges, satellite altimetry, ice shelf collapse, mass balance, TanDEM-X, The Cryosphere, glacier monitoring