Main

Surface meltwater will play a crucial role in the response of Antarctica to climate warming, because ponded water can promote ice-shelf breakup1,2, enhance localized melting3,4 and influence grounded ice dynamics5,6,7. Although first noted by Antarctic explorers in 19098, and subsequently identified at various locations using remote sensing9,10, the recent revelation that surface water is widespread around the Antarctic margin3,11 has spawned several studies of individual ice shelves12,13,14, groups of ice catchments15,16, and at ice-sheet spatial scales11,17,18. Antarctic-wide fluxes of meltwater production have been generated using satellites19,20,21 and models22,23. However, a continent-wide surface meltwater area dataset with high temporal resolution and decadal time span is lacking, meaning that the multi-annual and seasonal evolution of surface hydrology at a continental scale remains poorly constrained24. We address this by providing a monthly, long-term (15-yr) assessment of surface meltwater across Antarctica, enabling a comprehensive assessment of trends and spatial patterns in meltwater area.

To capture seasonal and interannual variability, we mapped surface meltwater across the Antarctic continent in Google Earth Engine (GEE)14 at monthly intervals. Every available Landsat 7 and 8 image was analysed between 2006 and 2021, totalling 133,497 images. To account for meltwater probably excluded as a result of variable image coverage and cloud cover, we used image visibility assessments (Methods) to scale up our mapped area results, generating ‘estimated’ surface meltwater area totals14. This approach produced a consistent time series of surface meltwater, spanning the entire continent. All meltwater areas stated in-text are ‘estimated’ rather than ‘observed’, as these are probably more representative of the true value. However, our scaling-up method has some inherent uncertainties and true meltwater area is probably somewhere between observed and estimated values (Supplementary Discussion 1). Our analysis produces the same conclusions regardless of which of these data are used (Supplementary Discussion 2).

Distribution of Antarctic meltwater

Ice-surface meltwater systems, comprising interconnected lakes, streams and regions of slush, occur around much of the Antarctic margin (Fig. 1). Surface water covered 3,732 km2 (s.d. 1,547 km2) of Antarctica on average (mean) each melt season, with 70% of water located on ice shelves. On the flatter topography of ice shelves, lakes are often elongate and connected by streams (Fig. 1a–i), with subtle along-flow topographic corrugations encouraging lakes to extend parallel with ice flow (Fig. 1e,g,i). Fewer, more isolated lakes occur in topographically confined basins on grounded ice and often reform in the same ___location in several years (Fig. 1b,i,g). The greatest surface meltwater area (annual maximum extent) was observed during the 2016/17 melt season, covering 4,574 and 2,218 km2 of floating and grounded ice, respectively. By contrast, during the 2011/12 melt season, surface meltwater only covered 1,004 km2 of floating ice and 516 km2 of grounded ice.

Fig. 1: Surface meltwater around the Antarctic continent.
figure 1

a–i, Occurrence frequency of surface meltwater for key ice shelves (IS) and outlet glaciers (G): George VI Ice Shelf (a), Amery Ice Shelf (b), Riiser-Larsen Ice Shelf (c), Rennick Glacier (d), Shackleton Ice Shelf (e), Nivlisen Ice Shelf (f), Nansen Ice Shelf (g), Pine Island Ice Shelf (h) and Bach Ice Shelf (i). j, The ___location of each subpanel is shown in the central map: the spatial density of surface meltwater features across Antarctica, shown as mean meltwater area per 1-km2 cell using a 50-km search radius. Circles associated with ice-shelf labels are sized proportionally to the average surface meltwater area of each named ice-shelf region, split between floating ice (grey) and grounded ice within 20 km of the grounding line (black). Pie charts are only displayed for ice-shelf regions with mean meltwater area >10 km2. The thick black lines show the extent of the East and West Antarctic ice sheets and the AP. The central grey square indicates the region around the South Pole that was not mapped due to a lack of Landsat images (Methods). Ice-sheet and ice-shelf boundaries from ref. 56.

Full size image

The East Antarctic Ice Sheet (EAIS) has the largest area of surface meltwater (2,812 km2 mean annual meltwater area; s.d. 1,478 km2), representing ~0.03% of total EAIS area (grounded and floating area combined). EAIS-wide totals are dominated by large lakes (up to ~90 km2) on ice shelves and their associated grounding zones, with 64% of total meltwater area situated on floating ice and 91% of water on grounded ice located within 20 km of the grounding line. Mean annual meltwater area totals are greatest for the Amery (1,110 km2; s.d. 791 km2), Roi Baudouin (516 km2; s.d. 402 km2), Shackleton (260 km2; s.d. 192 km2) and Nivlisen (138 km2; s.d. 118 km2) Ice Shelf regions, although with high variability between melt seasons (Extended Data Fig. 1). The proportion of each ice shelf covered in meltwater is nevertheless low; the highest being ~5% for Nivlisen Ice Shelf in 2007/08. The two largest EAIS ice shelves, Filchner and Ross (East), experience mean surface meltwater areas of 6 km2 (s.d. 5 km2) and 59 km2 (s.d. 32 km2), respectively.

The Antarctic Peninsula (AP) hosts the second highest meltwater area (855 km2 mean; s.d. 665 km2), although meltwater covers a greater proportion of the AP (~0.2%) than the EAIS (~0.03%). Complex meltwater networks often cover central sections of the George VI, Wilkins and Bach ice shelves (Fig. 1); the proportion of George VI Ice Shelf covered in meltwater each melt season (>2% median) is larger than any other ice shelf in Antarctica (Extended Data Fig. 1). While meltwater on the southern AP is concentrated on ice shelves, meltwater farther north on the AP is more abundant on grounded ice. For example, there is more meltwater on grounded ice surrounding the Larsen C Ice Shelf (<20 km inland of the grounding line), than on the ice shelf itself (Fig. 1j). Surface meltwater on grounded ice often includes extensive regions of slush (Extended Data Fig. 2g), which is abundant on the northeast outlet glaciers of the AP and glaciers feeding the southern AP ice shelves. Variability in the proportion of AP meltwater on grounded versus floating ice is far greater than either the West Antarctic Ice Sheet (WAIS) or EAIS, ranging from 1% on grounded ice in 2009/10 to 68% during the 2013/14 melt season (Extended Data Table 1).

Surface meltwater is less abundant on the WAIS, both in terms of absolute meltwater area (64 km2 mean; s.d. 22 km2) and proportionally (~0.002%). Peak surface meltwater area across the WAIS is typically 40–90 km2 each melt season, with 2012/13 having the highest surface meltwater area total (105 km2) of the study period (Fig. 2). Previously reported meltwater features around Pine Island, Cosgrove and Getz ice shelves17 recur throughout our time series (Fig. 1g). Surface lakes of up to ~2 km2 also commonly occur on grounded ice near nunataks surrounding Nickerson and Sulzberger ice shelves. On average, almost half (48% mean) of WAIS surface meltwater exists on grounded ice. No surface meltwater is detected around the eastern section of the Ross Ice Shelf or on any of the Siple Coast ice streams, despite the known occurrence of melt events within our study period25.

Fig. 2: The seasonal and interannual evolution of surface meltwater across the Antarctic ice sheets from 2006 to 2021.
figure 2

a–c, EAIS (a), AP (b) and WAIS (c). Meltwater area data for November–February each melt season. Each opaque-coloured bar shows the observed (O) mapped surface meltwater area, with translucent bars displaying scaled up estimated (E) meltwater area based on visibility assessments (Methods).

Full size image

Variability and trends in meltwater area

For the EAIS, WAIS and AP, meltwater area displays high interannual and seasonal variability. Around the EAIS, surface meltwater occurs throughout the austral summer (Fig. 2a) with a gradual increase in meltwater area throughout November and December, peaking in January, followed by a rapid decline during February (Fig. 2a). Annual variability in meltwater area is typically consistent between EAIS ice shelves, with notable exceptions18. For example, the greatest surface meltwater area for Nivlisen Ice Shelf occurred during the 2007/08 melt season when most of the EAIS experienced lower-than-average meltwater area (Fig. 3a). At the AP, there is typically negligible surface meltwater during November and December, before a rapid jump to a seasonal maximum in January. High meltwater area sometimes remains into February (Fig. 2b). The largest monthly meltwater area on the AP was in January 2020, consistent with record high meltwater area on George VI Ice Shelf26. Around the WAIS, surface meltwater area generally reaches a maximum in January, but peaked in December in 4 years when overall meltwater area was low (Fig. 2c).

Fig. 3: Interannual variability and trends in surface meltwater area from 2006 to 2021.
figure 3

a, Z-score anomalies relative to the 15-yr mean for the EAIS as a whole and for ice shelves with >0.5% meltwater coverage on average. b, Z-score anomalies as in a for the AP and individual ice shelves. c, Z-score anomalies as in a and b for the WAIS; note that no WAIS ice shelves have >0.5% meltwater coverage. Z-scores for 2018/19 are probably underestimated as a result of missing data in February 2019, highlighted by asterisks.

Full size image

For the EAIS, we observe an increase in both magnitude and variability of meltwater area across our 15-yr study period (Fig. 2a). Mean annual EAIS meltwater area between 2014 and 2021 (3,693 km2; s.d. 1,493 km2) was over double that of 2007–2013 (1,807 km2; s.d. 534 km2), despite 2015/16 and 2020/21 having low meltwater area (Fig. 2a). We find statistically significant increasing trends in EAIS meltwater area for our annual data (November–February maximum extent) and individually for the months of November and January (Extended Data Fig. 3). The greatest annual meltwater area occurred in 2016/17, when 5,945 km2 of the EAIS was covered in meltwater at some stage during the melt season. There is no overall increase in meltwater area for either the AP or the WAIS, and these regions display contrasting patterns in meltwater area variability (Fig. 2). The WAIS is characterized by consistently low meltwater area, while the AP displays large variability between melt seasons; annual meltwater area ranged from 57 km2 in 2013/14 to 2,330 km2 in 2019/20.

Links with Antarctic climate

To investigate first-order atmospheric controls on continental-scale surface meltwater area, we compared our data to three modes of atmospheric variability known to influence Antarctic surface air temperatures: (1) the Southern Annular Mode (SAM); (2) the El Niño/Southern Oscillation (ENSO); and (3) the strength and ___location of the Amundsen Sea Low (ASL). The climate impacts of the three modes vary across Antarctica27,28,29, so analysis was conducted separately for the EAIS, WAIS and AP. We also compared our data against modelled snowmelt from a statistically downscaled (2-km resolution) version30 of the regional climate model RACMO2.3p2 (ref. 23). This enabled us to identify spatial and temporal differences in the proportion of snowmelt that gets translated into ponded meltwater.

For the EAIS, we find a statistically significant negative correlation between the detrended summer SAM index and detrended maximum annual surface meltwater area (r = −0.47, P = 0.02), with greater meltwater area associated with negative SAM years (Fig. 4c). For the 6 years with the lowest SAM (<+0.7), median surface meltwater area was 20% greater than the 15-yr average (Fig. 4b). Conversely, the six lowest EAIS melt area years coincided with a strong positive summer SAM index (>+1.5) (Extended Data Table 2). EAIS surface meltwater area in the two most negative SAM years was 84% greater than the median (Extended Data Fig. 5d). The only year with a strong (<−1.5) negative SAM index (2016/17) corresponded with the highest annual meltwater area of the entire study period (Figs. 2a and 4a). No statistically significant correlations are found between EAIS surface meltwater area and ENSO or the ASL (Extended Data Fig. 4b–d).

Fig. 4: Links between Antarctic surface meltwater and climate for the EAIS and AP.
figure 4

a, Significant (P < 0.05) negative relationship (linear regression) between surface meltwater area and the summer SAM index for the EAIS (black dots, with grey bars extending to observed meltwater area), with regression (black line) shown (P = 0.02). b, Surface meltwater area anomalies across the EAIS in years with summer SAM < 0.7, showing widespread positive anomalies of >50%. c, Significant negative relationship between surface meltwater and ASL RCP on the AP (P = 0.001). The melt seasons of three distinct outliers with higher-than-expected surface meltwater are labelled. d, High percentage surface meltwater anomalies on the AP during years with ASL RCP <−8 (deep ASL years). Ice-sheet and ice-shelf boundaries in b and d from ref. 56.

Full size image

At the AP, the strength and ___location of the ASL appears closely linked to the distribution and abundance of surface meltwater, as also shown for surface air temperatures and sea-ice extent31. Summer ASL relative central pressure (RCP) is negatively correlated (r = −0.38, P < 0.001) with surface meltwater area (Fig. 4c). Years with a deep ASL (<−8 RCP) result, on average, in surface meltwater area anomalies 102% greater than the median (Fig. 4d), while a weaker (>−6.5) ASL results in lower-than-average meltwater ponding on the AP (Extended Data Fig. 6a). Surface meltwater area is also typically greater when the centre of the ASL is closer to the AP (Extended Data Figs. 4h and 6c,d), consistent with the correlation to wind and temperature28. While ASL strength has a distinct broad-scale influence on AP surface meltwater, the three highest meltwater years (09/10, 12/13 and 19/20) stand out as clear outliers to this relationship (Fig. 4c). These years did not have exceptionally eastward ASL longitudes (Extended Data Table 2), suggesting that local feedbacks and variability in surface meteorological conditions, and extreme weather events, are important in controlling the most intense melt events32,33. We do not observe statistically significant relationships between surface meltwater area and either SAM or ENSO on the AP (Extended Data Fig. 4e,f), perhaps due to the known dependency of AP temperatures on the combined impact of ENSO and the SAM34, weaker ENSO links in austral summer to AP temperatures34 and/or links between the eastern AP and tropical climate variability in other regions35. No statistically significant relationships are observed for the WAIS (Extended Data Fig. 4i–l), possibly due to meltwater area being too low for discernable signals to be identified.

We find positive significant power law relationships between monthly surface meltwater area and modelled snowmelt for both the EAIS (R2 = 0.37, P < 0.01) and the WAIS (R2 = 0.60, P < 0.01) (Fig. 5a). Perhaps surprisingly, given the large area difference between the ice sheets, the magnitude of monthly modelled snowmelt integrated across the EAIS is only 43% higher on average than for the WAIS. However, the amount of snowmelt that becomes ponded surface meltwater is one-to-two orders of magnitude greater for the EAIS than the WAIS (Fig. 5a). Hence, at an ice-sheet scale, melting in East Antarctica results in a proportionally far greater increase in surface meltwater area than in West Antarctica. For the EAIS, meltwater area increases at a greater rate than modelled melt across our study period, indicating increasing susceptibility of melt to pond on the ice surface. This increase is observed at both ice-sheet-wide and ice-shelf scales across the EAIS (Fig. 6). During the first half of our study period, relatively high cumulative modelled snowmelt totals (November–February sum) in 2009/10 (74 Gt), 2010/11 (66 Gt) and 2012/13 (89 Gt) translated into peak monthly surface meltwater areas of ~1,000–1,800 km2 (Fig. 5b). By contrast, modelled snowmelt totals of a similar magnitude in 2013/14 (83 Gt), 2016/17 (79 Gt) and 2019/20 (91 Gt) coincided with total surface meltwater area of ~2,900–4,500 km2 (Fig. 5b). At a regional scale, 14 out of 17 analysed EAIS ice-shelf regions show increasing surface meltwater area relative to modelled melt, seven of which show statistically significant trends across our study period (Fig. 6a–g). These include Shackleton and Nivlisen ice shelves, the third and fourth highest contributors to total meltwater area in East Antarctica.

Fig. 5: Comparison between surface meltwater area and modelled snowmelt.
figure 5

a, The relationship between monthly surface meltwater area and monthly RACMO modelled snowmelt, integrated across the EAIS (black), WAIS (blue) and AP (red). Statistically significant (P < 0.05) results from linear regression are displayed for the EAIS (P = 2.8 × 10−5) and WAIS (P = 7.5 × 10−13). The different symbols indicate the month (November–February) of each data point. b, Time series of surface meltwater area (black) and modelled snowmelt (red) across the EAIS, for the months of November–February each melt season.

Full size image
Fig. 6: Increasing sensitivity of East Antarctic ice-shelf regions to surface meltwater ponding.
figure 6

a, Map indicating the correlation strength between time (years since start of study period) and surface meltwater ponding susceptibility (quantified as annual surface meltwater area per unit Gt of modelled melt) for selected ice shelves around East Antarctica. The correlation strengths of ice-shelf regions with statistically significant relationships are displayed on a red–blue divergent colour-scale, with non-significant ice shelves shaded in light (positive correlation) or dark (negative) grey. b–h, Scatter plots for ice shelves labelled in a, which display a statistically significant increasing trend (linear regression) in meltwater ponding susceptibility throughout our study period: Brunt/Stancomb (b), Riiser-Larsen (c), Jelbart (d), Nivlisen (e), West (f), Shackleton (g) and Totten (h). Grey bars indicate the difference in result if using observed meltwater area. Note the large differences in y axis scale between plots. Ice-sheet and ice-shelf boundaries in a from ref. 56.

Full size image

At the AP, we observe hysteresis rather than a power law relationship between modelled melt and meltwater area, with surface meltwater area lagging behind modelled snowmelt (Fig. 5a). Relatively high snowmelt in December typically coincides with total meltwater areas of 5–100 km2, whereas surface meltwater areas of 300–800 km2 in February occur with less snowmelt. While a similar lag effect is observed for the EAIS and WAIS, these signals are less pronounced (Fig. 5a). There is far greater variability in AP surface meltwater area values relative to the same magnitude of melt than is observed for either the EAIS or WAIS. No clear patterns in ice-shelf susceptibility to melt were identified for the AP or WAIS.

Discussion

The likelihood of surface meltwater being able to pond on an ice surface is primarily controlled by the balance between snow accumulation and snowmelt36. Thicker snowpacks have a greater capacity for holding meltwater, meaning more melt is required to saturate the pore space and trigger surface ponding. This could explain the lagged effect on the AP (Fig. 5a), where high accumulation rates30 result in thick snowpacks, which take time to saturate with meltwater before ponding can occur. Despite having only marginally higher ice-sheet-wide integrated modelled melt rates, the EAIS hosts substantially more surface meltwater (Fig. 5a) than the WAIS, suggesting that, in general, WAIS snowpacks are currently less susceptible to meltwater saturation and surface ponding.

Our results show an overall increase in the magnitude and variability of surface meltwater area across the EAIS throughout our study period, despite ice-sheet-wide modelled snowmelt rates remaining relatively consistent (Fig. 5b). Regional analysis suggests that this increased sensitivity to meltwater ponding has occurred around the entire EAIS margin, but has been most pronounced at ice shelves in Dronning Maud Land and Wilkes Land (Fig. 6a). While there are several climatic and glaciological controls that could explain this increased sensitivity, we hypothesize that the cause was an increase in the abundance of low permeability ice surfaces in East Antarctica. Reduced ice-sheet surface permeability can occur for several reasons: enhanced wind scouring exposing ice surfaces14,37; multi-annual melting and freezing of water within the snowpack forming impermeable ice lenses38,39; and reductions in firn air content due to reduced accumulation or the snowpack becoming saturated with melt40,41. These processes have already been observed at regional scales in Antarctica18,38, and have resulted in widespread reductions in permeability across the Greenland ice sheet42. Our results suggest that a similar effect may be occurring in East Antarctica, although detailed analysis of in situ, modelled and satellite data would be required to assess the specific surface conditions for any given ___location. High variability in meltwater area (sometimes despite similar magnitudes of melt; Fig. 6b–h) between different EAIS ice shelves suggests strong regional variability in susceptibility to meltwater ponding, probably due to regional differences in firn thickness and saturation.

The associations between surface meltwater area and the SAM in East Antarctica, and the ASL on the AP, demonstrate the broad-scale influence that these modes of climate variability can have on melt43. At the AP, a deeper (lower pressure) ASL typically brings relatively warm, northerly air flows down the western side of the AP34, enhancing surface air temperature and the likelihood of melting44. At the EAIS, the negative correlation between the SAM and surface meltwater area is consistent with the known influence of the SAM on EAIS surface air temperatures29,43 and observed relationships between surface air temperature and meltwater ponding18. Extreme melt years that cannot be explained by these drivers may be related to localized factors, such as solar radiation44,45, foehn winds45,46 and low cloud cover47, or to extreme weather events, often driven by atmospheric rivers48. Although related to large-scale drivers, atmospheric rivers often have quite specific configurations32, which would not be detected by our correlation or composite analysis.

The presence of surface meltwater ponds can provide sufficient meltwater supply to initiate hydrofracture, either on floating49,50 or grounded ice7,51. While most ice-shelf regions which currently host surface meltwater are thought to be resilient to hydrofracture50, overall increases in surface meltwater area could result in expansion into more vulnerable areas, as has already been observed at Shackleton Ice Shelf18. The proportion of Antarctic ice shelves that are currently covered in surface meltwater is low (Extended Data Fig. 1), and future anticipated increases in melt52 may partly be mitigated by increased evacuation of meltwater off ice shelves53. Increases in surface meltwater on grounded ice would increase the likelihood of surface-to-bed connections developing, which could limit the amount of meltwater that drains from the ice sheet onto ice shelves. However, the delivery of surface meltwater to the bed could have important implications for grounded ice dynamics, as it can induce transient7 and seasonal54,55 accelerations in ice motion.

Methodological advances and computing capabilities have enabled us to produce a monthly, continent-wide dataset of Antarctic surface meltwater area for 2006–2021. Our results show that there has been an increase in the magnitude and variability of surface meltwater area in East Antarctica in recent years, which, in the absence of a coincident increase in modelled snowmelt, suggests that the ice-sheet surface is becoming increasingly prone to meltwater ponding. These results create an imperative for our surface meltwater dataset to be combined with modelled, satellite and in situ data, to better understand the controls and potential future impacts of meltwater ponding in regions most at risk to increases in surface meltwater area.

Methods

Surface meltwater mapping

Surface meltwater, ice, rock and cloud were automatically detected from Landsat 7 and 8 optical satellite imagery using a band-thresholding technique15. Cloud and rock areas were masked from image tiles, and an ice-specific version of the normalized difference water index (NDWIice) was subsequently used to delineate areas of surface meltwater. We applied the same thresholds as those developed by ref. 15, except for one minor change; we altered the threshold value for ‘green (B3)–red (B4)’ from >0.07 to >0.10 (Supplementary Table 3). This threshold is applied following NDWI classification to exclude areas of cloud shadow and shaded snow from surface meltwater classification results15. Our testing found that a threshold of >0.07 was insufficient at removing all rock and cloud shadow, resulting in widespread instances of misclassification error (Supplementary Fig. 4). By increasing this threshold to >0.10, we were able to successfully eradicate these misclassifications, while still retaining correctly identified surface meltwater within our dataset (Supplementary Discussion 3).

Mapping was undertaken in GEE, using Landsat level 1 tier 2 top-of-atmosphere images stored in the GEE data catalogue. By using a cloud-based computational platform, we were able to conduct rapid processing over spatial and temporal scales that would have otherwise been unachievable14. Surface meltwater was mapped monthly between 2006 and 2021. Before 2006, image coverage was insufficient to generate comprehensive continent-wide data. Only images that had a sun elevation angle >20° were processed15, hence limiting surface meltwater data to austral summer months. This filtering step was taken to avoid misclassification errors, since surface water is difficult to differentiate spectrally in low-light conditions57. Surface meltwater is rare outside austral summer months regardless58. Our mapping procedure incorporated a robust method for assessing image visibility, enabling us to quantify variability in image coverage and cloud cover14. Using this approach, we estimated the maximum area of surface hydrology that would be expected under cloud-free conditions14, in addition to the observed mapped values. A discussion of the uncertainties associated with this method are provided in Supplementary Discussion 1. We mapped both monthly and ‘annual’ (November–February, to cover the melt season) maximum area of surface meltwater, for ease of comparison with climate datasets.

A total area of ~12.32 million km2 was covered during surface meltwater mapping, including both the grounded ice sheet and surrounding ice shelves. Every Landsat image covering this study area between 2006 and 2021 was used during analysis, totalling 133,497 images. Our mapping procedure involved automatically looping through an Antarctic-wide grid of 1,151 region-of-interest (ROI) shapefile tiles, ordinarily 108 × 108 km2 in area. This tile size maximized spatial coverage for mapping, while remaining within the memory capacity of an individual task within GEE. Landsat image coverage does not extend to latitudes greater than ~85° south; hence ROI tiles that overlapped this area (~1.28 million km2) were not mapped (Fig. 1). However, given the high elevation and very low temperature of the region around the South Pole, little meltwater ever exists at this ___location. The grid was clipped to the Antarctic coastline56,59, meaning coastal ROI tiles varied in shape and area. We did not account for changes in ice-shelf area throughout our study period. Our results show that surface meltwater rarely ponds adjacent to ice shelf calving fronts (it typically exists further inland near the grounding line); hence we deemed that minor changes in ice-shelf area will have had negligible influence on mapped meltwater totals. Tiles were given a unique ID based on longitude and latitude for identification purposes. In instances where the coastline clipping process split a tile into multiple portions, tile segments were merged to adjacent tiles to ensure no two tiles had the same ID.

Processing was performed on a yearly basis for up to ~350 ROI tiles at a time. Memory capacity and timeout limits within GEE were exceeded when attempting to process larger regions than this. Following the method of ref. 14, vector outputs were exported from GEE as geoJSON files. Post-processing steps were then undertaken in MATLAB to clean the raw data and to produce final shapefile outputs. For full details of the lake detection, image visibility assessment and post-processing methods applied here, see ref. 14.

Surface meltwater quality control

Surface meltwater quality control was conducted to identify any artifacts within the continent-wide dataset. Annual maximum extent surface meltwater shapefiles were merged to produce a single continent-wide map. Quality control was then performed via manual inspection of this 15-yr maximum extent raster layer (Extended Data Fig. 7), as this offered the quickest method for flagging invalid meltwater polygons. Polygons that (1) existed in highly improbable locations (for example, at high elevations in the cold interior of the ice sheet), (2) had straight or cornered edges or (3) formed an ordered or repeated pattern were searched for during manual inspection. Polygons that satisfied some or all of the above criteria were highlighted for verification and checked against the optical images from which they originated. In all cases highlighted by the above process, flagged lake polygons were the result of artifacts in Landsat 7 imagery and were removed from the final dataset. In total, 699 polygons were removed during this process, totalling 9.87 km2 of incorrectly identified meltwater area.

Owing to the large spatial and temporal scale of our study, it was not realistic to conduct systematic manual verification of each individually mapped instance of surface meltwater against its corresponding optical imagery. The thresholds we applied to delineate surface meltwater, modified from ref. 15, were established on the basis of a range of spectral conditions around Antarctica. Their verification against optical imagery showed an accuracy of >95% (ref. 15) and mapped outputs showed high similarity with meltwater area data produced by alternative methods57. Our minor adjustment to their method (Supplementary Discussion 3) reduces the likelihood of cloud and rock shadow misclassifications, providing further confidence in the accuracy of our meltwater classifications. Furthermore, ref. 14 showed that the thresholds were suitable for use on both Landsat 7 and 8 imagery (Supplementary Discussion 4). It should be noted that our dataset only captures surface meltwater identifiable in the visible spectrum and cannot be used to make conclusions about melt rates or subsurface meltwater storage. The 30-m resolution of Landsat additionally means that we will not have captured the smallest-scale surface meltwater features, such as streams or small lakes. Areas of slush (which are particularly common on ice shelves and sections of grounded ice on the northern AP, Extended Data Fig. 2g,h) identified by the thresholds have been retained within the dataset. However, while our dataset contains some areas of slush, the thresholds we applied were designed specifically for mapping surface meltwater (not slush); hence, alternative methods60,61 would be required to comprehensibly map slush across Antarctica.

While systematic method performance was not assessed, sample mapped results were visually compared against optical imagery when trialling our method in different regions of Antarctica. Overall, the thresholds were extremely effective at identifying surface meltwater across all regions of Antarctica, with <1% of mapped meltwater pixels identified as misclassification error based on 100 randomly sampled features across Antarctica. Errors were highest in regions with dirty ice (such as McMurdo Ice Shelf), high crevasse shadow or variable slush presence (Extended Data Fig. 2). However, misclassified meltwater features were typically small (<5 pixels in size) and limited in spatial extent, and were therefore deemed to have minimal influence on regional or ice-sheet scale results (Supplementary Discussion 5). On the basis of all the stated reasons, we hence considered the method to be sufficiently accurate for continent-wide application, without needing to adjust thresholds for different regions.

Statistical analyses

We generated regional surface meltwater area totals by aggregating mapped meltwater areas across the EAIS, WAIS and AP, and for individual ice shelves and their associated grounding zones. Ice-sheet scale (EAIS, WAIS and AP) and ice-shelf regions were defined using basins from ref. 59. Grounded catchments and ice shelves were merged in ArcMap on the basis of their ice-sheet classification, to produce polygons for the EAIS, WAIS and AP. Ice-shelf regions included the ice shelf as defined by ref. 59 plus a 20-km buffer onto grounded ice inland from the grounding line. We took this additional step to ensure that surface meltwater in the grounding zone of an ice shelf, where it is typically most abundant7, was included in regional analysis.

We assessed changes in meltwater area over time using robust linear regression in MATLAB (https://uk.mathworks.com/help/stats/robustfit.html). Our data are unusual in that we only have data during 4 months of the year; hence we have repeated seasonal gaps in our time series. Regression analyses, using time as the independent variable, were therefore used instead of traditional trend analysis techniques (such as the Mann–Kendall test), which require data to have no serial correlation or seasonality. Data were tested to check that they met the assumptions of robust linear regression. Robust regression is less sensitive to outliers than standard linear regression, allowing us to test for the presence of broad-scale trends across our study period, without being disproportionally skewed by one or two outlier years. Regression analyses were performed on individual monthly surface meltwater area data (November–February) and on annual maximum extent data (Extended Data Fig. 3). Analyses were run on both our mapped and estimated meltwater area data, to investigate how the results varied (Supplementary Discussion 2). Standardized Z-scores of annual surface meltwater area totals were additionally calculated in MATLAB, enabling quantification of surface meltwater area for individual years relative to the study period mean (Fig. 3).

Comparison with climate data

Climatic indices for three modes of Antarctic climate variability were compared against our surface meltwater area data: (1) the SAM index, from the observation-based index following ref. 62; (2) the Oceanic Niño index, one measure for the ENSO, from the National Oceanic and Atmospheric Administration (NOAA); and (3) ASL indices of RCP and longitude, following ref. 63. RCP data were used in preference to ‘actual’ central pressure data to isolate changes in the ASL from the SAM, which strongly modulates the actual pressure across all seasons29. Each climate index provides a simple diagnostic quantity, which is used to characterize the state of the geophysical system in question. For all the datasets, monthly index values for December, January and February were averaged to generate single, austral summer values. These ‘annual’ values were then compared against the annual maximum extent of surface meltwater. Climatic and surface meltwater data were detrended before conducting robust regression analysis, which enabled us to better assess the statistical relationship for each individual mode of climatic variability, without results being skewed by anonymously high- or low-melt years, which could have been strongly influenced via teleconnections with other climatic modes64.

We performed composite analysis29,65, a commonly used statistical technique to identify potential impacts of climatic phenomena, to explore the large-scale impact of modes of climatic variability on surface meltwater area. For each climatic mode, we determined thresholds to categorize subsets of ‘low’ and ‘high’ yearly values. Thresholds were chosen on the basis of visual inspection of plotted time series for each climatic index, with thresholds selected to separate clear ‘peaks’ and ‘troughs’ in the data. We produced composite surface meltwater fields for each category by averaging the surface meltwater totals for each subset of years. Mean surface meltwater area totals were calculated on a tile-by-tile basis, using the Antarctic-wide grid used during mapping.

RACMO comparison

We compared our meltwater area data to a statistically downscaled version30 (2-km resolution) of the regional climate model RACMO2.3.p2 (ref. 23). Reference 30 produced this statistically downscaled product by refining the spatial distribution of surface mass balance components using high-resolution (2 km) datasets, including the reference elevation model of Antarctica (REMA) and an albedo map from the moderate resolution imaging spectroradiometer (MODIS). These data therefore offer higher spatial resolution estimates of meltwater production than previous modelled datasets and are shown to have closer agreement with in situ and satellite records30. Monthly meltwater totals were integrated across each respective ice-sheet area to produce monthly totals of meltwater production for the EAIS, WAIS and AP (Fig. 5). We additionally compared RACMO snowmelt against surface meltwater area for each individual ice-shelf region (Fig. 6) and each ROI tile used during mapping, to assess spatial variability in agreement between the two datasets. For any tile where meltwater was mapped in at least 5 months throughout the study period, we correlated monthly meltwater area against RACMO snowmelt (Extended Data Fig. 8).