Iceberg imaging and sampling
During expedition PS126 of RV Polarstern, the opportunity arose to visit a large iceberg in the vicinity of HAUSGARTEN carrying numerous dark stones. The iceberg was accessed by helicopter from the ship on 14 June 2021 at the coordinates 78° 35.66′ N, 3° 32.92′ W, allowing us to collect samples of the transported stone material.
A series of overlapping, downward-facing images of the iceberg and its stone load were recorded from the helicopter at an altitude of about 100 m using a Canon EOS digital camera with a lens operated at about 50 mm. The images were stitched together using Agisoft Photoscan Pro. Stone piles placed approximately 5 m apart by field personnel served as a size reference for both individual images and the final mosaic (Extended Data Fig. 2a). Additional close-range images of stone clusters were acquired on the iceberg using a Nikon D850 with a 50 mm lens, with a ruler included in each frame for scale. The planar area of 20 randomly selected stones per image (n = 500 in total) was measured using ImageJ. Furthermore, stones were sampled to determine their mineralogical composition. Representative specimens spanning the observed morphological range were collected for classification (Extended Data Fig. 2c).
An ice core (9 cm in diameter) was collected from the upper 2 m of the iceberg (after snow removal) using a hand auger with an electric drill (Kovacs Enterprise). The coring location was about 50 m away from the nearest stone pile to avoid stone interference during drilling. Three replicate measurements of electrical conductivity were made from sub-cores 27–57 cm long (n = 6) to determine the salinity of the ice. Identical methods were used to assess the salinity of two sea ice floes nearby for reference (78° 54.11′ N, 3° 9.39′ W and 79° 01.44′ N, 5° 42.49′ W).
Seafloor observations
Images were recorded from the seafloor at HAUSGARTEN station EG-IV (2,500 m depth) using either the OFOS or the Ocean Floor Observation and Bathymetry System between 2015 and 2021 (ref. 49). The exact transect location varied by year due to ice conditions (Extended Data Fig. 1b). In 2015 and 2017, the same transect was surveyed, allowing for a comparison between years in an area of seafloor about 5,400 m2. OFOS included a downward-facing camera (Canon EOS), strobes that flashed to illuminate each image and downward-facing laser points that served as a size scale for the image and stone plan area. Telemetry data showed distance to the seafloor, and all images were recorded at a target altitude of 1.5 m. OFOS images were recorded automatically every 30 s while the ship transited at 0.5 knots. Images that were too bright, too dark or at an anomalous altitude (that is, outside the range of 1.3– 1.6 m) were considered ineligible for analysis.
We observed glacial dropstones on the seafloor in OFOS images, ranging from 0.6 cm to 73 cm. The vast majority of stones on the deep seafloor north of 45° N are glacial in origin50. Near-bottom currents in the bathyal Fram Strait are <10 cm s−1, not fast enough to uncover buried stones or bedrock14,21.
An exception occurs at a steep reef in the eastern Fram Strait, where stronger near-bottom currents expose bedrock and mobilize locally derived stones51. However, this feature lies several hundred kilometres from station EG-IV, where our imagery was collected. The stones observed on both the seafloor and the visited iceberg consist of shale and quartz (Extended Data Fig. 2). Taken together, the geological, bathymetric and oceanographic context indicates that stones observed on the seafloor at station EG-IV represent glacial dropstones.
The densities of stones in 30 randomly selected images from the beginning, middle and end of the transect were calculated each year (total of 90 images per year). Dropstones and dropstone-associated fauna were annotated in each image using the online tool BIIGLE (Bio-Image Indexing and Graphical Labelling Environment; biigle.de; ref. 52). OFOS images recorded at 3 m altitude reliably show seafloor features and fauna >1 cm across49. Images in this study were recorded at about 1.5 m altitude, allowing for visualization of fauna and features to 0.6 cm across. Fauna visible in our OFOS images included sponges, bryozoans, tunicates, anemones, soft corals, serpulid polychaetes, barnacles, sea stars and crinoids. Some taxa could be identified based on taxonomic voucher samples collected in previous studies throughout Fram Strait53. Taxa for which no voucher had been identified were given morphotype descriptors. The densities of stones and of each dropstone-associated species in 2015 and 2017 were calculated by dividing the number in each image by the planar area of that image. Plan areas of 250 stones in randomly selected 2015 and 2017 images (n = 250 per year) were measured using the rectangle select tool in BIIGLE. Plan areas of 500 stones from randomly selected 2021 images were measured using the straight line tool for measurement of length and width in ImageJ.
Dropstones are typically deposited on the seafloor in clusters, and slight variations in the swing of the OFOS camera could cause these clusters to appear in one year but not another. To avoid biasing the results by a small number of images containing dense stone clusters, we removed outlier images with anomalously high stone densities (>5 stones per m2) before analysis (n = 1 image, 2015; n = 3 images, 2017). Differences in univariate metrics (that is, stone density, fauna density) were evaluated using parametric analysis of variance (ANOVA) tests or non-parametric Kruskal–Wallis tests, when the assumptions of ANOVA were violated. Chi-square tests were used to test for differences in the size distributions of stones on the iceberg and on the seafloor and between years (2015 and 2017). Statistical analyses were conducted in the R environment using the packages vegan, pairwiseAdonis and ggplot2.
Ship-based glacial ice observations
Small icebergs, growlers and bergy bits generally evade detection by satellite imagery, necessitating visual observations for monitoring their presence. On board the research vessel Polarstern, iceberg observations have been routinely carried out every 3 h by the German Weather Service (DWD). These observations are part of the Surface Synoptic Observations, a globally standardized meteorological data collection program maintained according to the guidelines set by the World Meteorological Organization (WMO). This continuous observational practice has been operational since the commissioning of Polarstern in 1984, providing essential information for onboard activities such as navigation, route planning and aircraft operations, while simultaneously contributing to global meteorological networks.
Synoptic observations from ships include key atmospheric parameters, current weather conditions and basic oceanic parameters, such as air and water temperature, air pressure, wind speed and direction, humidity, precipitation, cloud cover and cloud types, visibility, wave height and periodicity, all encoded in a globally standardized reporting format ‘FM 13 SHIP’54. In ice-covered waters, observations include details on the ice-going ability of the vessel, the type and developmental stage of encountered sea ice, and sightings of glacial ice classified according to WMO54 alphanumeric table 0439/BUFR table 020035. For glacial ice, this code specifically describes the type (icebergs, bergy bits and growlers) and quantity of observed glacial ice, reported in increments of five (see legend in Extended Data Fig. 5). The absence of glacial ice and the inability to perform observations due to restricted visibility conditions are also explicitly noted.
The processing and analysis of glacial ice observations were conducted as follows. Initially, all synoptic records collected between 1981 and 2024 were compiled separately by campaign (data for individual campaigns are accessible through PANGAEA55). Subsequently, all valid glacial ice observations within the Fram Strait region (70°–85° N, 30° W–30° E) were filtered. Observations were defined as valid if visibility exceeded 1 km. This resulted in approximately 8,200 valid observations, with glacial ice reported around the vessel in about 32% of these cases. The spatial distribution of observations is shown in Fig. 3a. Note that repeated sightings at the same position or within a radius smaller than the visibility at the time of observation were excluded to avoid double-counting.
Following this filtering, campaign-specific and annual mean frequencies of iceberg sightings were calculated (Fig. 3b), as well as frequency distributions of individual glacial ice types. Extended Data Fig. 5 shows the frequencies of glacial ice observations, along with the corresponding classification codes. Most frequently observed were 1–5 icebergs accompanied by growlers and bergy bits (class 6), followed by occurrences of only icebergs (1–5, class 1) or only growlers and bergy bits (class 4).
Visual glacial ice observations inherently involve uncertainties affecting their representativeness. First, sightings may vary among individual observers, particularly under challenging visibility conditions or low sun angles. Moreover, observer performance might evolve over time, for example, as observers gain experience. Although this subjective error is difficult to quantify precisely, it is considered small because of the consistent use of trained meteorological experts. Second, weather conditions such as precipitation or fog can limit the number of visual observations carried out during expeditions. On average, about 20.4% of planned observations are hindered by visibility conditions below the 1 km threshold, although this percentage can be considerably higher for individual expeditions. However, analysis of the 40-year observational record shows no significant trend in visibility conditions, indicating that no systematic long-term biases are present in the observational time series (Extended Data Fig. 6). Third, larger and more prominent icebergs are proportionally detected more frequently than smaller growlers and bergy bits. This uncertainty also remains difficult to precisely quantify.
Given these limitations, definitive quantification of absolute iceberg frequencies, interannual variability or spatial distribution patterns remains uncertain and should be interpreted cautiously. Nevertheless, Fram Strait, and specifically the HAUSGARTEN region, benefits from frequent observation efforts, on average traversed by the AWI 2.5 times per year, yielding a continuous and extensive dataset. The density and duration of this observational time series enable robust detection of relative changes in glacial ice occurrence despite the discussed observational constraints.
Backtracking of glacial ice
To identify the source regions (calving sites) of glacial ice and, consequently, the origin of stones deposited on the deep seafloor, we backtracked ship-based glacial ice observations using satellite-derived ice drift data. For this purpose, we used IceTrack56, a tool originally developed for backtracking sea ice to the sites where it was formed.
For the backtracking approach, ice is traced backwards on a daily basis from observed locations using three independent low-resolution sea ice drift products: (1) OSI-405-c ice-motion data from the Ocean and Sea Ice (OSI) Satellite Application Facility57 (SAF); (2) MERGED ice-motion vectors provided by the Center for Satellite Exploitation and Research58 (CERSAT); and (3) Polar Pathfinder daily motion vectors (v.4.1) available from the National Snow and Ice Data Center59 (NSIDC). Backtracking is terminated whenever trajectories intersect coastlines, encounter fast-ice edges, or when the satellite-derived ice concentration (from OSI SAF products OSI-430-a and OSI-450-a) falls below 50%.
This Lagrangian tracking method is applicable to glacial ice as long as icebergs, growlers and bergy bits are firmly embedded within the pack ice and carried along by it. Under looser ice-pack conditions or for deep-drafted icebergs, the validity of this assumption diminishes. Here, ocean currents begin to strongly influence the drift, causing icebergs to diverge from the general motion of the surrounding, mainly wind-driven, sea ice60. To account for this limitation, we restricted backtracking to sightings of glacial ice firmly entrained within the pack ice at the time of observation and throughout their drift. Specifically, we imposed a criterion that satellite-derived ice concentration along each trajectory should not fall below 80% for more than 30 cumulative days. A total of 106 out of the 1,362 glacial ice sightings fulfilled this condition, with their derived drift tracks shown in Fig. 3a.
The applied concentration criterion is grounded in an evaluation of the tracking approach against independent drift observations. First, trajectories of 50 GPS-equipped sea ice buoys deployed in the Arctic Ocean between 2012 and 2016 (Meereisportal.de data repository61) were reconstructed using IceTrack and directly compared with their observed GPS positions. The mean great-circle deviation between reconstructed and observed locations amounts to approximately 47 ± 48 km after 100 days and 109 ± 68 km after 300 days (Extended Data Fig. 10a), indicating reliable reproduction of large-scale drift pathways under compact pack-ice conditions. To assess tracking performance under looser ice conditions, iceberg drift was additionally evaluated using a small set of GPS-equipped icebergs that were partially exposed to open-water conditions. The applied dataset comprises 11 icebergs tracked between 2012 and 2025 in the Beaufort Gyre and north of the Canadian Arctic Archipelago (data sources: International Arctic Buoy Programme Ice Island Archive, https://iabp.apl.uw.edu/Ice_Islands_2025.html; Carleton University Ice Island Drift Tracking Database, https://wirl.carleton.ca/research/ice/ice-islands/ibtd/). The comparison between reconstructed and observed trajectories indicates larger deviations (55 ± 70 km after 100 days and 286 ± 123 km after 300 days), particularly during periods when sea ice concentration falls below 80% (Extended Data Fig. 10b). Despite the higher deviations under loose ice conditions, positional deviations remain sufficiently small to distinguish iceberg source regions at basin and sector scale, although attribution to individual marine-terminating glaciers is not supported, particularly for trajectories involving extended drift durations.
To assess whether the increasing number of icebergs observed in Fram Strait could, in part, be explained by an enhanced input of glacial ice from the Canadian Arctic, we applied the same Lagrangian tracking framework. Sea ice crossing a virtual gate at 83° N (25° W–25° E) was traced backwards in time using an identical Lagrangian backtracking configuration to that described in ref. 45. The Beaufort fraction was defined as the proportion of tracked trajectories that intersected the Beaufort Gyre region prior to export through Fram Strait (Extended Data Fig. 9). Temporal changes in this fraction were analysed to evaluate whether the contribution of Beaufort-origin ice increased after 2000.
Glacier changes
In this study, we relate the observed increase in iceberg frequency to changes at the contributing marine-terminating glaciers. These changes can be quantified through several complementary approaches that differ in physical meaning, temporal resolution and data requirements. Although the main outlet glaciers in northeast Greenland are comparatively well observed, many marine-terminating glaciers in the Russian High Arctic lack the consistent, long-term datasets required for detailed dynamic assessments62.
Frontal ablation
Calving is most directly reflected in changes in frontal ablation. Estimates of frontal ablation rely on high-resolution satellite imagery used to track terminus positions over time, typically through manual or semi-automated delineation. For northeast Greenland, these observations document the extensive retreat and fragmentation of the ZI and NG glaciers since the early 2000s (Extended Data Fig. 4, derived from Landsat 5–8 optical imagery), providing an indication of enhanced calving over the past two decades. By contrast, for many Russian High Arctic glaciers, terminus positions have been mapped only episodically or over relatively short time windows35,39,63.
Ice discharge
A second approach quantifies ice discharge, defined as the flux of grounded ice across a gate near the grounding line. Discharge isolates the dynamic component of glacier change and, over multi-decadal timescales, typically correlates with changes in calving activity. Although discharge is not identical to frontal ablation, increased discharge generally reflects accelerating flow and dynamic thinning that favour enhanced calving. This method requires ice velocity, ice thickness and bedrock topography and is therefore available only on longer timescales for the major outlet glaciers in northeast Greenland. For ZI and NG, we use the discharge estimates in ref. 64, which extend back to 1987 and allow quantification of long-term changes in dynamic mass loss. In Fig. 3c, cumulative discharge estimates relative to the stable period from 1985 to 2000 are shown, providing a robust baseline against which the post-2000 increase can be evaluated.
Mass-change estimates
Where neither frontal-ablation nor discharge records exist, glacier change can be assessed through mass-change estimates derived from remote-sensing products. These measurements reflect the combined effects of surface mass balance and dynamic thinning. For northeast Greenland, we calculate mass changes using NASA Operation IceBridge ATM data and satellite altimetry (ICESat, ICESat-2, CryoSat-2 and ENVISAT), referenced to a 1978 DEM32,65. Extended Data Fig. 7 presents the time series of mass loss for the lower sectors of ZI and NG32. The two mass-change time series indicate that there was no significant ice loss between 1978 and 2000, an assumption that is also applied when calculating the cumulative discharge shown in Fig. 3c. For SZ and FJL, we rely on published mass-change assessments based on DEM differencing, radar and optical altimetry36,38,39,62. Although these studies differ in methodological foundations, resolution and temporal coverage, they consistently report sustained mass loss of marine-terminating glaciers since 2000.
Surface velocities
Where direct estimates of frontal ablation or ice discharge are unavailable, glacier dynamics can be assessed using satellite-derived surface velocity time series as indicators of changes in ice flow. We use glacier surface velocities for selected marine-terminating glaciers on FJL and SZ from the NASA MEaSUREs ITS_LIVE project, generated using the auto-RIFT feature-tracking algorithm48,66. Each velocity estimate shown in Fig. 4 and Extended Data Fig. 8 represents a temporally averaged surface speed over the image-pair separation interval. Seasonal mean velocities were computed using a moving temporal window centred on 1 July of each year (±30 days). Within each window, velocities were averaged using weights proportional to the image-pair temporal separation to account for unequal temporal coverage.
FESOM2-iceberg drift model
To quantify potential changes in the large-scale transport of icebergs towards Fram Strait, we conducted a controlled sensitivity experiment using the FESOM267 equipped with an interactive iceberg module68,69. The experiment prescribes temporally constant calving fluxes from predefined glacier sectors in NEGIS, FJL, and SZ, allowing us to examine changes in iceberg drift dynamics and associated processes along their trajectories over the simulation period.
FESOM2 is a global, multi-resolution ocean-sea ice model that supports regional mesh refinement in areas of interest. In this study, we use a model configuration with a horizontal resolution of approximately 4.5 km in the Arctic Ocean and Nordic Seas. The model is forced with the JRA-55 atmospheric reanalysis70 and has previously been evaluated against observational datasets for Arctic Ocean circulation and sea ice dynamics71,72,73.
The iceberg module represents icebergs as Lagrangian particles with prescribed geometric properties, including height, length and width. These dimensions evolve dynamically during the simulation through thermodynamic processes such as basal and lateral melting, wave-induced erosion and buoyant convection. The momentum balance of model icebergs includes Coriolis and gravitational forces as well as drag forces exerted by the atmosphere, ocean and sea ice. Sea ice drag follows the step-function parameterization described in ref. 74. For sea ice concentrations below 15%, sea ice drag is set to zero. When sea ice concentration exceeds 90%, and sea ice strength is greater than 10,000 N m−3, icebergs are assumed to be frozen into the pack ice and drift with the surrounding sea ice. A detailed description of the iceberg module is provided in refs. 68,69.
In the Southern Ocean, simulated iceberg drift, both in open water and while embedded in pack ice, has been evaluated against observed iceberg trajectories derived from satellite75,76 and buoy datasets77 in refs. 68,69. However, in the Arctic, direct validation is difficult because of the limited availability of observational data. Nevertheless, icebergs that become embedded in dense pack ice are expected to drift largely with the surrounding sea ice cover. Because FESOM2 realistically reproduces observed Arctic sea ice drift and variability71, we assume that the drift of icebergs entrained in pack ice is captured realistically by the model.
Whenever the simulated iceberg draft exceeds the local bathymetry, the iceberg is assumed to be grounded. Basal melting continues during grounding and can eventually reduce the draft sufficiently for the iceberg to refloat. Discrete icebergs are generated from an integrated calving flux. The prescribed ice discharge is converted into an equivalent calving area by dividing the flux by the assumed iceberg thickness (250 m; refs. 78,79). Individual iceberg sizes are sampled from a −2.2 power-law size distribution following ref. 80. Icebergs are assumed to have a quadratic surface and cuboidal geometry. For icebergs smaller than 250 m in horizontal extent, height is set equal to the horizontal dimensions, whereas larger icebergs are assigned a maximum height of 250 m to reduce the risk of instantaneous grounding after release. In an iterative adjustment procedure, iceberg dimensions are adjusted so that the total iceberg volume matches the prescribed integrated discharge.
To isolate transport dynamics from variations in iceberg production, iceberg release was prescribed with temporally constant annual fluxes throughout the simulation period (1970–2023), thereby preventing changes in prescribed calving rates from influencing the simulated iceberg distribution. Icebergs were released continuously throughout the year rather than at a fixed time of release to avoid introducing artificial seasonal biases. Release locations were defined in sufficiently deep offshore waters rather than directly at glacier grounding lines, thereby minimizing artificial grounding effects and reducing the influence of poorly constrained near-coastal processes such as fast-ice storage or mélange-controlled release.
Icebergs were released in three source sectors identified by the backtracking analysis (Fig. 3a): NEG, FJL and SZ. These regions represent the primary sectors from which icebergs may enter Fram Strait. The simulations, therefore, provide a physically consistent large-scale representation of the dominant drift corridors linking these source regions Fram Strait. For the three source regions, the following integrated calving fluxes, averaged over the time period 2000–2010, were prescribed:
-
SZ: 3.15 gigatonnes per annum (Gt a−1) (ref. 63)
-
FJL: 8.78 Gt a−1(ref. 63)
-
NEG: 11.15 Gt a−1 (ZI64); 8.41 Gt a−1 (NG81); 0.33 Gt a−1 (NEG peripheral glaciers63).
The first 16 years (1970–1985) of the simulation are considered a spin-up period to reach an equilibrated number of active icebergs, as suggested in ref. 79. To separate the effect of changing ocean and sea ice conditions on iceberg drift from the increase in active icebergs during the spin-up phase, the analysis of iceberg dynamics focuses on the period from 1986 to 2023. To quantify spatial patterns of iceberg transport (Fig. 5), we calculate iceberg transit frequency as the relative frequency of iceberg passages through each grid cell over the simulation period.
