The Greenland ice sheet (GrIS) stores freshwater equal to more than 7 m of potential sea level rise (SLR) and strongly interacts with the Arctic, North Atlantic and global climate. Over the last few decades, the ice sheet has been losing mass at a rate that is projected to increase. Interactions between the GrIS and the climate have the potential to amplify or reduce GrIS mass balance responses to ongoing and projected warming. Here, we investigate the impact of ice sheet–climate interactions on the climate and mass balance of the GrIS using the Community Ice Sheet Model version 2 (CISM2) coupled with the Community Earth System Model version 2 (CESM2). To this end, we compare two idealized multi-century simulations with a non-evolving and evolving ice sheet topography in which we apply an annual 1 % increase in CO2 concentrations, starting from pre-industrial (PI) until stabilization at 4×PI CO2 concentrations (4×CO2). By comparing the one- and two-way coupled simulations, we find significant changes in atmospheric blocking, precipitation and cloud formation over Greenland as the GrIS topography evolves, acting as negative feedbacks on mass loss. We also attribute part of the overestimation of mass loss in the one-way coupled simulation to an overestimation of melt in the ablation area caused by the use of a uniform temperature lapse rate to reflect the elevation differences between the atmospheric and ice sheet grids. Furthermore, we investigate ice sheet–climate interactions in a simulation branched in year 350 from our two-way coupled simulation in which we annually reduce atmospheric CO2 by 5 % until PI concentrations are reached. During the 350-year 4×CO2 forcing period, the ice sheet loses a total mass of 1.1 m sea level equivalent, and part of its margins retreat landward. When the PI CO2 concentration is restored, melt decreases rapidly, leading to a small positive surface mass balance. Combined with the strongly reduced ice discharge resulting from the widespread retreat of the ice sheet margin, this halts GrIS mass loss despite a remaining global warming of 2 K. The GrIS, Arctic and North Atlantic strongly interact, causing a complex transitional phase towards a colder climate during the century following the CO2 reduction. Elevated atmospheric temperatures, larger ocean heat transport and deteriorated state of the snowpack, compared to the initial pre-industrial state, result in limited regrowth of the ice sheet under reintroduced PI CO2 conditions.

Historical aerial imagery provides valuable data from regions and periods with limited geospatial information. A common method to utilize this data is through the generation of ortho-photos and 3D models using Structure-from-Motion (SfM) techniques. However, many of these images were scanned decades after their acquisition and require geometric calibration, along with internal and external camera parameter estimation, for accurate reconstruction. Manual identification of key features, such as fiducial marks and text annotations, is labour-intensive, while existing automated methods struggle with poor-quality datasets. This paper presents an automated workflow that combines computer vision and machine learning techniques to detect and extract these key features from historical aerial images. To address challenges related to image quality, we also introduce estimation protocols that compensate for missing or unreliable detections by leveraging redundancy across multiple flight paths. The methodology was evaluated on the TMA (Trimetrogon Aerial) archive, a collection of historical images from the Antarctic Peninsula. Our test dataset comprised over 7000 images from 20 different flight paths. The workflow demonstrated high success rates in detecting and extracting fiducial marks, image subsets, and textual annotations. Approximately 70% of the images provided usable focal length data, while fiducial mark detection exhibited high accuracy except in cases of severe scanning artifacts. Altitude data extraction proved to be the most challenging, with successful results in only 15% of images due to degraded altimeter readings. Despite these limitations, the automated workflow effectively estimated missing parameters, ensuring robust image reconstruction across flight paths. The code for this workflow is open-source and publicly available on GitHub at https://github.com/fdahle/hist_meta_extraction.

The instability of the West Antarctic Ice Sheet (WAIS) is a tipping element in the climate system, and it is mainly dictated by changes in the ice flow behaviour of the outflow glaciers in the Amundsen Sea Embayment (ASE). Recent studies postulated that the vertical uplift of bedrock can delay the collapse of glaciers in this region. In West Antarctica, bedrock motion is largely caused by a fast viscoelastic response of the upper mantle to changes in ice loads over the last centuries. This glacial isostatic adjustment (GIA) effect is currently poorly understood, since Earth's rheology and the ice-loading history are both subject to large uncertainties in simulations. Moreover, results from data-driven approaches have not yet resolved GIA at a sufficient spatial resolution. We present a data-driven GIA estimate, based on data from GRACE/GRACE-FO (GRACE and GRACE-FO), CryoSat-2 altimetry, regional climate modelling, and firn modelling, which is the first to agree with independent vertical velocities in West Antarctica derived from global navigation satellite system (GNSS) data. Our data combination yields a maximum GIA bedrock motion rate of 43 ± 7 mm a−1 in the Thwaites Glacier region and agrees within uncertainties in the GNSS-derived rate. The data-driven GIA-related bedrock motion may be used in future simulation runs to quantify a potential delay of the collapse of the West Antarctic Ice Sheet due to the stabilization effects induced by GIA. Furthermore it may be used for testing rheological models with low upper-mantle viscosity in conjunction with centennial loading histories.

In recent decades, satellite radar altimetry has been widely used to assess volume changes over the Greenland Ice Sheet. In particular, melt events result in drastic changes in the volume scattering of firn, which induces a pronounced change in the parameters derived from radar altimeter data. Due to the recent and increasingly frequent melt events over Greenland, the impacts of these events on the firn condition, i.e. the formation of ice lenses and reduction in firn air content, need to be better understood. This study therefore exploits the ability of long-term CryoSat-2 data to indicate changes in firn volume scattering in order to assess the spatiotemporal firn condition variations in Greenland. More specifically, this study utilises the leading edge width (LeW) parameter derived from CryoSat-2 Low Resolution Mode data, which has been proven to be a parameter strongly sensitive to changes in volume scattering, and assesses its variation between September 2010 and September 2024. With a combined analysis of remote sensing observations, in situ observations, and outputs from regional climate models, our study demonstrates that the LeW drop induced by extreme melt events in the interior of Greenland experiences a gradual recovery, which can potentially be explained by new-snow deposition. However, in many high-elevation regions of Greenland where firn layers were originally dry, the recent recurrence of extensive melt has prevented a full recovery of the firn volume scattering to pre-2012 conditions, indicating a persistent increase in firn density under a changing climate. Finally, our study also confirms the utility of radar altimeter data for long-term monitoring of the impact of melt and refreezing events on the properties of the upper firn layer

Patagonian glaciers have been rapidly losing mass in the last two decades, but the driving processes remain poorly known. Here we use two state-of-the-art regional climate models to reconstruct long-term (1940-2023) glacier surface mass balance (SMB), i.e., the difference between precipitation accumulation, surface runoff and sublimation, at about 5 km spatial resolution, further statistically downscaled to 500 m. High-resolution SMB agrees well with in-situ observations and, combined with solid ice discharge estimates, captures recent GRACE/GRACE-FO satellite mass change. Glacier mass loss coincides with a long-term SMB decline (−0.35 Gt yr⁻²), primarily driven by enhanced surface runoff (+0.47 Gt yr⁻²) and steady precipitation. We link these trends to a poleward shift of the subtropical highs favouring warm northwesterly air advections towards Patagonia (+0.14°C dec⁻¹ at 850 hPa). Since the 1940s, Patagonian glaciers have lost 1350 ± 449 Gt of ice, equivalent to 3.7 ± 1.2 mm of global mean sea-level rise.

Damage features, such as rifts and crevasses, are the first signs of a weakened ice shelf and the precursor for retreat. Yet, damage changes are not widely quantified on Antarctic ice shelves, leaving future ice shelf weakening poorly understood. Here we use satellite imagery to detect both long-term (24-year) and short-term (annual, 2015–2021) Antarctic-wide damage changes, revealing a multiyear damage development cycle strongly correlated to ice shelf area changes, and a net decline in damaged area from 1997 to 2021. We establish a data-driven link between damage and ice flow characteristics, which shows that ice flow acceleration, strain rate increases and thinning lead to more damage development, in particular under high-emission climate scenarios. This sensitivity to warming suggests that without quantification of damage impacts by detailed physical models the (timing of) ice shelf retreat and Antarctic mass loss may currently be underestimated.

Glaciers are indicators of ongoing anthropogenic climate change¹. Their melting leads to increased local geohazards², and impacts marine³ and terrestrial^4,5 ecosystems, regional freshwater resources⁶, and both global water and energy cycles^7,8. Together with the Greenland and Antarctic ice sheets, glaciers are essential drivers of present^9,10 and future^{11, 12–13} sea-level rise. Previous assessments of global glacier mass changes have been hampered by spatial and temporal limitations and the heterogeneity of existing data series^{14, 15–16}. Here we show in an intercomparison exercise that glaciers worldwide lost 273 ± 16 gigatonnes in mass annually from 2000 to 2023, with an increase of 36 ± 10% from the first (2000–2011) to the second (2012–2023) half of the period. Since 2000, glaciers have lost between 2% and 39% of their ice regionally and about 5% globally. Glacier mass loss is about 18% larger than the loss from the Greenland Ice Sheet and more than twice that from the Antarctic Ice Sheet¹⁷. Our results arise from a scientific community effort to collect, homogenize, combine and analyse glacier mass changes from in situ and remote-sensing observations. Although our estimates are in agreement with findings from previous assessments^{14, 15–16} at a global scale, we found some large regional deviations owing to systematic differences among observation methods. Our results provide a refined baseline for better understanding observational differences and for calibrating model ensembles^12,16,18, which will help to narrow projection uncertainty for the twenty-first century^11,12,18.

Most of the Greenland and Antarctic ice sheets are covered with firn — the transitional material between snow and glacial ice. Firn is vital for understanding ice-sheet mass balance and hydrology, and palaeoclimate. In this Review, we synthesize knowledge of firn, including its formation, observation, modelling and relevance to ice sheets. The refreezing of meltwater in the pore space of firn currently prevents 50% of meltwater in Greenland from running off into the ocean and protects Antarctic ice shelves from catastrophic collapse. Continued atmospheric warming could inhibit future protection against mass loss. For example, warming in Greenland has already contributed to a 5% reduction in firn pore space since 1980. All projections of future firn change suggest that surface meltwater will have an increasing impact on firn, with melt occurring tens to hundreds of kilometres further inland in Greenland, and more extensively on Antarctic ice shelves. Although progress in observation and modelling techniques has led to a well-established understanding of firn, the large uncertainties associated with meltwater percolation processes (refreezing, ice-layer formation and storage) must be reduced further. A tighter integration of modelling components (firn, atmosphere and ice-sheet models) will also be needed to better simulate ice-sheet responses to anthropogenic warming and to quantify future sea-level rise.

Despite in-situ observations of perennial firn aquifers (PFAs) at specific locations of the Antarctic ice sheet, a comprehensive continent-wide mapping of PFA distribution is currently lacking. We present an estimate of their distribution across Antarctica in the form of a probability assessment using a Monte Carlo technique. Our approach involves a novel methodology that combines observations from Sentinel-1 and Advanced SCATterometer (ASCAT) with output from a regional climate model. To evaluate our method, we conduct an extensive comparison with Operation Ice Bridge observations from the Greenland Ice Sheet. Application to Antarctica reveals high PFA probabilities in the Antarctic Peninsula (AP), particularly along its northern, northwestern, and western coastlines, as well as on the Wilkins, Müller, and George VI ice shelves. Outside the AP, PFA probability is low, except for some locations with marginally higher probabilities, such as on the Abbot, Totten, and Shackleton ice shelves.

Correction to: Nature Reviews Earth & Environment https://doi.org/10.1038/s43017-023-00507-9, published online 23 January 2024.

In the version of the article initially published, in Fig. 5, under “Radar altimeter”, “O(16–160 m)” previously read “O(16–160 km)”. This has now been corrected in the HTML and PDF versions of the article.

Recently, climate extremes have been grabbing attention as important drivers of environmental change. Here, we assemble an observational inventory of energy and mass fluxes to quantify the ice loss from glaciers on the Russian High Arctic archipelago of Novaya Zemlya. Satellite altimetry reveals that 70 ± 19% of the 149 ± 29 Gt mass loss between 2011 and 2022 occurred in just four high-melt years. We find that 71 ± 3% of the melt, including the top melt cases, are driven by extreme energy imports from atmospheric rivers. The majority of ice loss occurs on leeward slopes due to foehn winds. 45 of the 54 high-melt days (>1 Gt d−1) in 1990 to 2022 show a combination of atmospheric rivers and foehn winds. Therefore, the frequency and intensity of atmospheric rivers demand accurate representation for reliable future glacier melt projections for the Russian High Arctic.

Assessing the Surface Mass Balance (SMB) of the Antarctic Ice Sheet is crucial for understanding its response to climate change. Synthetic Aperture Radar observations from Sentinel-1 provide the potential to monitor the variability of SMB processes through changes in the scattering response of near-surface and internal snow layers. However, the interplay between several factors, such as accumulation, wind erosion, deposition, and melt, complicates the interpretation of scattering changes of the microwave signal. Additionally, lack of reliable ground truth measurements of the snow surface limits our capability to associate the SMB processes with dominant scattering mechanism. In this study, we aim to quantify the dominant scattering in Sentinel-1 signal and evaluate the scattering changes in drifting snow-dominated regions of East Antarctica. We introduce a scattering indicator, alpha -{text{scat},varepsilon }, derived from scattering-type and entropy descriptors, providing a measure between volume and pure scattering. By relating the field measurements to alpha -{text{scat},varepsilon }, we establish that the evolution of dominant scattering in the presence of snowdrift is complex. First, alpha -{text{scat},varepsilon } strongly correlates with surface roughness (R^{2}=0.92, RMSE =2^circ). Spatially variable erosion patterns significantly increase the roughness and result in a strong affinity towards pure scattering despite net accumulation. Second, high surface densities also tend to influence pure scattering; however, the effect is dependent on the accumulation rate. With more accumulation, we observe an increasing dominance of volume scattering from internal snow layers. Long-term trends in alpha -{text{scat},varepsilon } (2017/2023) further suggest that it is challenging to address the causes behind the scattering source based on a single snow surface process. We thus demonstrate the potential and limitations of alpha -{text{scat},varepsilon } to infer the variability in dominant scattering from changes in surface processes.

The utility of historical image repositories is often limited due to the lack of geo-referencing. A good example is the TriMetrogon Aerial (TMA) archive, a collection of historical aerial images of Antarctica between 1940 and 2000. These images are difficult to use, as their geolocation is only approximately, with location errors in the order of tens of km. This study addresses this challenge by developing an automated geo-referencing workflow that leverages recent advancements in machine-learning-based tie-point matching. We use the algorithm LightGlue, to establish tie-points between geo-referenced Sentinel-2 satellite imagery and historical non-geo-referenced aerial images. To aid the process, we use already known approximate positions of the historical images. Due to the sub-optimal and inhomogeneous quality of the aerial images, only a portion of the images can be geo-referenced directly by matching. For the remaining images, we employed alternative means of geo-referencing, again based on tie-point matching. Out of a subset of 4,459 images located inside the research area, 3,393 images could be geo-referenced, a percentage of 76%. Reasons for the geo-referencing failing are insufficient contrast or an approximate position too far away from the real position. The workflow can easily be applied to historical images from other archives, to enhance the usability of historical image repositories for scientific research.

The TriMetrogon Aerial (TMA) archive is an archive of historical images of Antarctica taken by the US Navy between 1940 and 2000 with analogue cameras. The analysis of such historic data can give a view of Antarctica's glaciers predating modern satellite imagery and provide unique insights into the long-term impact of changing climate conditions with essential validation data for climate modelling. However, the lack of semantic information for these images presents a challenge for large-scale computer-driven analysis. Such information can be added to the data using semantic segmentation, but traditional algorithms fail on these scanned historical grayscale images, due to varying image quality, lack of colour information and artefacts in the images. To address this, we present a deep-learning-based U-net workflow. Our approach includes creating training data by pre-processing and labelling the raw images. Furthermore, different versions of the U-net are trained to optimize its hyperparameters and augmentation methods. With the optimal hyper-parameters and augmentation methods, a final model has been trained for a use-case to segment 118 images covering Adelaide Island. We tested our approach by segmenting challenging historical images using a U-net model with just 80 training images, achieving an accuracy of 73% for 20 validation images. While no test data is available for our use case, a visual examination of the segmented images shows that our method performs effectively. The comparison of the hyper-parameters and augmentation methods provides directions for training other U-net-based models so that the presented workflow can be used to segment other archives with historical imagery. Additionally, the labelled training data and the segmented images of the test are publicly available at https://github.com/fdahle/antarctic_segmentation.

Because Antarctic surface melt is mostly driven by local processes, its simulation necessitates high-resolution regional climate models (RCMs). However, the current horizontal resolution of RCMs (≈25–30 km) is inadequate for capturing small-scale melt processes. To address this limitation, we present SUPREME (SUPer-REsolution-based Melt Estimation over Antarctica), a deep learning method to downscale surface melt to 5.5 km resolution using a physically-informed super-resolution model. The physical information integrated into the model originates from observations tied to surface melt, specifically remote sensing-derived albedo and elevation. These remote sensing data, in addition to a Regional Atmospheric Climate Model (RACMO) run at 27 km resolution, account for the diverse drivers of surface melt across Antarctica, facilitating effective generalization beyond the training region of the Antarctic Peninsula. A comparison of SUPREME with a dynamically downscaled RACMO run at 5.5 km over the Antarctic Peninsula shows high accuracy, with average yearly RMSE and bias of 5.5 mm w.e. yr⁻¹ and 4.5 mm w.e. yr⁻¹, respectively. Validation at five automatic weather stations reveals SUPREME's marked improvement with substantially lower average RMSE (81 mm w.e.) compared to RACMO 27 km (129 mm w.e.). Beyond the training region, SUPREME aligns more closely with remote sensing products associated with surface melt than super-resolution models lacking physical constraints. While further validation of SUPREME is needed, our study highlights the potential of super-resolution techniques with physical constraints for high-resolution surface melt monitoring in Antarctica, providing insights into the impacts of localized melting on processes affecting ice shelf integrity such as hydrofracturing.

The joint ESA/NASA Mass-change And Geosciences International Constellation (MAGIC) has the objective to extend time-series from previous gravity missions, including an improvement of accuracy and spatio-temporal resolution. The long-term monitoring of Earth’s gravity field carries information on mass change induced by water cycle, climate change and mass transport processes between atmosphere, cryosphere, oceans and solid Earth. MAGIC will be composed of two satellite pairs flying in different orbit planes. The NASA/DLR-led first pair (P1) is expected to be in a near-polar orbit around 500 km of altitude; while the second ESA-led pair (P2) is expected to be in an inclined orbit of 65^◦–70^◦ at approximately 400 km altitude. The ESA-led pair P2 Next Generation Gravity Mission shall be launched after P1 in a staggered manner to form the MAGIC constellation. The addition of an inclined pair shall lead to reduction of temporal aliasing effects and consequently of reliance on de-aliasing models and post-processing. The main novelty of the MAGIC constellation is the delivery of mass-change products at higher spatial resolution, temporal (i.e. subweekly) resolution, shorter latency and higher accuracy than the Gravity Recovery and Climate Experiment (GRACE) and Gravity Recovery and Climate Experiment Follow-On (GRACE-FO). This will pave the way to new science applications and operational services. In this paper, an overview of various fields of science and service applications for hydrology, cryosphere, oceanography, solid Earth, climate change and geodesy is provided. These thematic fields and newly enabled applications and services were analysed in the frame of the initial ESA Science Support activities for MAGIC. The analyses of MAGIC scenarios for different application areas in the field of geosciences confirmed that the double-pair configuration will significantly enlarge the number of observable mass-change phenomena by resolving smaller spatial scales with an uncertainty that satisfies evolved user requirements expressed by international bodies such as IUGG. The required uncertainty levels of dedicated thematic fields met by MAGIC unfiltered Level-2 products will benefit hydrological applications by recovering more than 90 per cent of the major river basins worldwide at 260 km spatial resolution, cryosphere applications by enabling mass change signal separation in the interior of Greenland from those in the coastal zones and by resolving small-scale mass variability in challenging regions such as the Antarctic Peninsula, oceanography applications by monitoring meridional overturning circulation changes on timescales of years and decades, climate applications by detecting amplitude and phase changes of Terrestrial Water Storage after 30 yr in 64 and 56 per cent of the global land areas and solid Earth applications by lowering the Earthquake detection threshold from magnitude 8.8 to magnitude 7.4 with spatial resolution increased to 333 km.

Ice shelves play a pivotal role in stabilizing the Antarctic ice sheet by providing crucial buttressing support. However, their vulnerability to basal melting poses significant concerns for ice sheet and shelf stability. Our study focuses on assessing basal melt rates at a 50 m posting of 12 ice shelves where earlier studies have identified high melt rates. We make use of the Reference Elevation Model of Antarctica (REMA) strips to generate surface elevation- and melt rates using the Basal melt rates Using Rema and Google Earth Engine (BURGEE) methodology. BURGEE reveals higher melt rates in areas with thinner ice than existing remote sensing basal melt products. This is for instance the case for basal channels on both Dotson, Totten and Pine Island ice shelves. Modelling studies have already shown that remote sensing inferred basal melt rates are underestimated at the thinnest part of basal channels, and that this underestimation scales with resolution coarsening. Since the thinner parts of an ice shelf also represent its weakest part, it is crucial that we capture its melting well to fully grasp the vulnerability of the ice shelf. Our work, therefore, represents a crucial step in uncovering the vulnerability of Antarctic ice shelves. By exposing detailed melting patterns, particularly in areas like basal channels, we highlight not just extensive melting but also potential weak points, significantly contributing to our understanding of ice shelf stability. These findings bear substantial importance in comprehending the broader implications of ongoing climate changes on Antarctica's ice sheet integrity and, consequently, global sea levels.

The intrusion of Circumpolar Deep Water in the Amundsen and Bellingshausen Sea embayments of Antarctica causes ice shelves in the region to melt from below, potentially putting their stability at risk. Earlier studies have shown how digital elevation models can be used to obtain ice shelf basal melt rates at a high spatial resolution. However, there has been limited availability of high-resolution elevation data, a gap the Reference Elevation Model of Antarctica (REMA) has filled. In this study we use a novel combination of REMA and CryoSat-2 elevation data to obtain high-resolution basal melt rates of the Dotson Ice Shelf in a Lagrangian framework, at a 50 m spatial posting on a 3-yearly temporal resolution. We present a novel method: Basal melt rates Using REMA and Google Earth Engine (BURGEE). The high resolution of BURGEE is supported through a sensitivity study of the Lagrangian displacement. The high-resolution basal melt rates show a good agreement with an earlier basal melt product based on CryoSat-2. Both products show a wide melt channel extending from the grounding line to the ice front, but our high-resolution product indicates that the pathway and spatial variability of this channel is influenced by a pinning point on the ice shelf. This result emphasizes the importance of high-resolution basal melt rates to expand our understanding of channel formation and melt patterns. BURGEE can be expanded to a pan-Antarctic study of high-resolution basal melt rates. This will provide a better picture of the (in)stability of Antarctic ice shelves.

Our awareness of ice caps' and mountain glaciers' sensitivity to climate change has driven major advances in the application of remote sensing techniques during the past decade. Regarding ESA's SARIn altimeter CryoSat-2, processing the full waveform to generate swaths of elevation estimates has become standard practice in regions of complex topographies. This technique provides information on areas where we would be blind otherwise. In this article, we discuss systematic errors and analyze their impact on surface elevation measurements and change rates of two test areas. In particular, we focus on periodically occurring errors in elevation swaths, caused by the superposition of coherent signals from range-ambiguous surfaces. They can lead to measurement errors in excess of 10 m, affect most measurements in mountainous regions, are difficult to exclude with established post-processing techniques, and occur repeatedly for satellite revisits introducing a 369-day periodicity-difficult to distinguish from the annual cycle. We show a correlation between derived elevation swaths and the sensor view angle and explore the influence of common data exclusion choices on higher level products. Our results indicate that these systematic errors hold a substantial share of the error budget and that the choice of thresholds impacts higher level products. We conclude that error correlations need to be considered to characterize the data accuracy. With the established data editing strategies, systematic errors prevent resolving seasonal mass changes of single mountain glacier basins and impact aggregates over larger areas or longer periods.

Geoscientific models are facing increasing challenges to exploit growing datasets coming from remote sensing. Universal differential equations (UDEs), aided by differentiable programming, provide a new scientific modelling paradigm enabling both complex functional inversions to potentially discover new physical laws and data assimilation from heterogeneous and sparse observations. We demonstrate an application of UDEs as a proof of concept to learn the creep component of ice flow, i.e. a nonlinear diffusivity differential equation, of a glacier evolution model. By combining a mechanistic model based on a two-dimensional shallow-ice approximation partial differential equation with an embedded neural network, i.e. a UDE, we can learn parts of an equation as nonlinear functions that then can be translated into mathematical expressions. We implemented this modelling framework as ODINN.jl, a package in the Julia programming language, providing high performance, source-to-source automatic differentiation (AD) and seamless integration with tools and global datasets from the Open Global Glacier Model in Python. We demonstrate this concept for 17 different glaciers around the world, for which we successfully recover a prescribed artificial law describing ice creep variability by solving ∼ 500 000 ordinary differential equations in parallel. Furthermore, we investigate which are the best tools in the scientific machine learning ecosystem in Julia to differentiate and optimize large nonlinear diffusivity UDEs. This study represents a proof of concept for a new modelling framework aiming at discovering empirical laws for large-scale glacier processes, such as the variability in ice creep and basal sliding for ice flow, and new hybrid surface mass balance models.