Large subglacial source of mercury from the southwestern margin of the Greenland Ice Sheet

The Greenland Ice Sheet is currently not accounted for in Arctic mercury budgets, despite large and increasing annual runoff to the ocean and the socio-economic concerns of high mercury levels in Arctic organisms. Here we present concentrations of mercury in meltwaters from three glacial catchments on the southwestern margin of the Greenland Ice Sheet and evaluate the export of mercury to downstream fjords based on samples collected during summer ablation seasons. We show that concentrations of dissolved mercury are among the highest recorded in natural waters and mercury yields from these glacial catchments (521–3,300 mmol km−2 year−1) are two orders of magnitude higher than from Arctic rivers (4–20 mmol km−2 year−1). Fluxes of dissolved mercury from the southwestern region of Greenland are estimated to be globally significant (15.4–212 kmol year−1), accounting for about 10% of the estimated global riverine flux, and include export of bioaccumulating methylmercury (0.31–1.97 kmol year−1). High dissolved mercury concentrations (~20 pM inorganic mercury and ~2 pM methylmercury) were found to persist across salinity gradients of fjords. Mean particulate mercury concentrations were among the highest recorded in the literature (~51,000 pM), and dissolved mercury concentrations in runoff exceed reported surface snow and ice values. These results suggest a geological source of mercury at the ice sheet bed. The high concentrations of mercury and its large export to the downstream fjords have important implications for Arctic ecosystems, highlighting an urgent need to better understand mercury dynamics in ice sheet runoff under global warming. Meltwaters from the southwestern margin of the Greenland Ice Sheet contain exceptionally high concentrations of mercury, exporting up to more than 200 kmol of dissolved mercury every year, suggest mercury measurements from three glacial catchments.

M ercury (Hg) is a toxic element of global concern with limited biological function. Hg bioaccumulates in organisms and biomagnifies in aquatic food webs largely as the neurotoxin methylmercury (MeHg, CH 3 Hg). Overexposure to Hg, primarily due to consumption of seafood, has major environmental and human health implications with a socio-economic cost estimated to exceed US $5 billion per year 1 . The global Hg problem has worsened significantly due to anthropogenic pollution 2,3 , although there are natural sources of Hg to the environment such as volcanic emissions and weathering of Hg-bearing minerals in rock (for example, cinnabar) 2 . Hg is of particular concern in the Arctic 1 , where Hg concentrations in marine biota have increased by an order of magnitude over the past 150 years (ref. 4 ). This has led to health risks in local communities, where diets rely heavily on marine animals 1 . The Arctic is particularly vulnerable to anthropogenic Hg perturbations because it may be a global Hg sink as prevailing atmospheric circulation carries Hg to northern latitudes 5 . Although anthropogenic pollution is a major contributor to the elevated Hg concentrations in Arctic biota, inputs from climatically vulnerable naturally occurring pools have also received more recognition during the last decade [6][7][8][9][10] . For example, Arctic rivers are a significant Hg source to the Arctic Ocean and are climatically sensitive due to their intensifying hydrological cycles (variability and magnitude), widespread catchment permafrost cover and disproportionate warming in the Arctic 9,11 . As yet, little consideration has been given to the Greenland Ice Sheet (GrIS) in the Arctic Hg cycle.
The GrIS is the second largest body of ice on Earth, covering nearly 25% of the land surface area in the Arctic region. Observed rapid and accelerating GrIS mass loss is producing increased annual freshwater runoff to coastal waters, which is predicted to nearly double in magnitude by 2100 (ref. 12 ). The biogeochemical implications of this increasing runoff are uncertain but likely important 13,14 .
Only a handful of studies of Hg concentrations from Arctic proglacial rivers exist from the Canadian High Arctic Archipelago and disconnected glaciers (that is, those not part of the ice sheet) in east Greenland [15][16][17] , all showing moderate and climatically sensitive Hg concentrations in glacial runoff. Glaciers in the Himalayas and Alaska have also been highlighted as Hg stores and sources to downstream ecosystems [18][19][20][21] , with all studies implying an atmospheric or ice-marginal Hg source. There are no data from runoff draining the large, polythermal-based catchments of the GrIS. Determining the role of the GrIS in the Arctic Hg cycle is therefore crucial in assessing Arctic ecosystem health both now and into the future.
Here, we evaluate Hg and MeHg concentrations in meltwaters originating from the southwestern margin of the GrIS and export to downstream fjords (Fig. 1). The study focuses on a >4,000 km 2 region of the GrIS that covers three glacial catchments ranging from ~85 to 3,200 km 2 in size (Russell Glacier (RG), Leverett Glacier (LG) and Isunnguata Sermia (IS)), and three fjord systems (Nuup Kangerlua (NK), Ameralik Fjord (AF) and Søndre Strømfjord (SS)), which receive substantial GrIS meltwater inputs ( Fig. 1 and Methods). All GrIS catchments sampled are warm based from the ice sheet margin and have well-developed subglacial hydrological systems, likely with regions of anoxia 22,23 . Data are presented for total Hg (THg), particulate Hg (pHg, >0.45 µm), total dissolved Hg (dHg, <0.45 µm), size-fractionated dissolved Hg (cHg = 0.02-0.45 µm, sHg = <0.02 µm), dissolved inorganic Hg 2+ (dIHg, <0.45 µm) and dissolved methylmercury (dMeHg, <0.45 µm) over four field seasons. Complementary metagenomic sequencing of DNA from the subglacial microbial community sampled from subglacial meltwaters exiting LG is used to hypothesize Hg biogeochemical pathways through the identification of hgcAB, merA and merB genes. Our results provide evidence for the importance of ice sheet runoff in regional Hg cycling and indicate that coastal Hg dynamics are likely to change under climatic warming scenarios.

Mercury concentrations and speciation in meltwater rivers
Measurements of Hg in GrIS subglacial runoff over all sampling periods and sites were similar to contaminated and urban watersheds but much higher than those reported for glacierized and non-glacierized environments elsewhere (Table 1 and Supplementary Table 1 Table 1). Discharge-weighted mean concentrations of dIHg (Hg 0 + Hg 2+ ) were 603-1,260 pM (LOADEST modelled mean concentrations were 607-1,310 pM; Supplementary Table 2). There appear to be systematic differences in the dIHg concentrations observed between glacial catchments that may be explained by underlying bedrock type (Extended Data Fig. 1). Two of the glaciers (IS and RG) located predominantly on the Nagssugtoqidian mobile belt have significantly lower (P < 0.05; n = 28) dIHg concentrations (mean 632 and 603 pM, respectively) than LG, which predominantly overlies an Archaean block geological belt to the south (1,260 pM; Extended Data Fig. 1). The larger glacial catchments (LG and IS) show a gradually increasing trend in dIHg over time (R 2 = 0.20 (P < 0.15; n = 12) and 0.54 (P < 0.05; n = 14), respectively), concurrent with increasing discharge and seasonal subglacial hydrological system expansion (Fig. 2). The smallest catchment (RG) displays relatively invariant dIHg concentrations and discharge during the observation period (Fig. 2)  concentrations reported here are one to two orders of magnitude greater than dHg concentrations in Arctic rivers (1.5-14 pM; Supplementary Table 1), and more than an order of magnitude higher than Arctic river THg concentrations (28-74 pM; Table 1). The differences between GrIS meltwaters and Arctic rivers are also reflected in extremely high dIHg catchment yields (1,440-2,310 mmol km −2 year −1 ; Table 1), approximately two orders of magnitude higher than large Arctic rivers (4-20 mmol km −2 year −1 ). Only polluted rivers and estuaries in Asia have been shown to have similar or higher concentrations of dIHg (up to 50,000 pM; Supplementary Table 1) 24,25 . Methylmercury concentrations from the three catchments ( Fig. 2d) indicate substantial export of bioavailable Hg and potential Hg cycling pathways in subglacial ecosystems. Concentrations of dMeHg are relatively high from all three glacial catchments (>5 pM), exceeding values from most pristine freshwaters, and similar to those found in wetlands (below detection in rivers to 2.5 pM in the Everglades) 2,16,26 . Dissolved dMeHg accounts for 0.74 ± 0.17% (LG), 1.5 ± 0.97% (IS) and 2.1 ± 0.43% (RG) of the total dHg concentration, which is in the range of other freshwater systems (~1-10%) 2 , including Arctic rivers (~0.2-5 %) 10 , but lower than marine environments (4-22%) 27 . The higher concentrations of dMeHg in RG outflow coincide with higher mean dissolved organic carbon (DOC) concentrations compared with LG and IS (16 ± 2 µM versus 12 ± 4 µM, respectively; Extended Data Fig. 2), and higher methane (CH 4 ) concentrations compared with LG 22,28 . These differences suggest that variations in subglacial hydrology and/or organic matter content and supply beneath each catchment are important for dMeHg cycling. In addition, a greater fraction of meltwaters may originate from hypoxic/anoxic regions of the glacier bed at RG (hypothesized due to the high CH 4 concentrations 22,29 ) where potential for methylation is higher 1,30 . There may also be an influence of river-marginal inputs in the RG catchment (for example, proglacial lakes and tundra) given that the meltwater river travels along the ice margin for ~9 km before arriving at the sampling site. High concentrations of dMeHg in the glacial catchments sampled here therefore point to spatially heterogeneous subglacial water sources and hydrological flow pathways, a psychrophilic mercury cycling microbial community capable of mercury cycling and/or abiotic dIHg methylation 2,31 . The dMeHg concentrations also have detrimental implications for bioaccumulation of ice sheet-derived Hg in Arctic ecosystems. Collectively, the Hg concentrations here lead to first-order estimates of dHg fluxes from the southwestern region of the GrIS (Extended Data Fig. 3 (Table 1), which are within the same order of magnitude as Arctic rivers and with much higher yields (382 mmol km −2 year −1 dIHg and 4.3 mmol km −2 year −1 dMeHg; Table 1), highlighting the importance of previously unquantified ice sheet runoff in the Arctic Hg cycle.
A 2015 time series dataset from May to the end of July from LG ( Fig. 3) indicates temporal variability and distinctive size fractionation in dHg. Concentrations followed the seasonal hydrological drainage evolution of LG and are divided into four time periods 22 22 corresponding to sHg increasing from ~1,800 pM to ~4,200 pM, (3) the 'outburst' period, where the subglacial hydrological system was perturbed by water inputs from supraglacial meltwater drainage events 22,33 (day 170-190; 1,000-5,000 pM sHg and 2,000-18,000 pM cHg) and (4) the post-outburst period (day 190 onward), where discharge was largely controlled by the intensity of surface melt and Hg concentrations fell to a seasonal minimum before increasing, with the decline in discharge observed around day 208 (Fig. 3). There appeared to be a weak relationship between outburst events and Hg concentration. For example, sHg was elevated during outburst events 1, 2 and 4 (blue regions in Fig. 3) compared with immediately before, similar to trends observed in CH 4 and DOC (refs. 22,32 ), indicating periodic flushing of isolated subglacial waters high in Hg. A rise in Hg concentrations during a period of decreasing discharge in period 4 reinforces a previous hypothesis that lower subglacial water pressure in subglacial channels as discharge falls results in drainage of concentrated channel-marginal waters from more isolated regions of the subglacial environment into main drainage channels along a hydraulic gradient 23,33 . Total dHg was dominated by cHg species (0.02-0.45 μm; mean 82% of dHg), similar to other trace elements in ice sheet meltwaters 34 . High colloidal Fe and Al concentrations (>1 µM) 34 indicate cHg species are likely adsorbed to inorganic nanoparticle minerals (for example, Fe and Al (oxy)hydroxides).

Mercury cycling genes in the subglacial microbiome
DNA metagenomic analysis from LG meltwater suggests that subglacial microbial communities are genetically equipped to deal with Hg, consistent with other Arctic environments 35 . Microbes have been shown to cope with Hg 2+ toxicity by transforming it into volatile Hg 0 using mercuric reductase encoded by the merA gene and respond to MeHg stress by using the organomercurial lyase encoded by the merB gene to de-methylate MeHg into Hg 0 (and CH 4 ) 36,37 . Analysis of the microbial metagenomes at LG during 2015 showed the presence of merA and merB genes (Extended Data Fig. 4 displays taxa bearing the merA gene, a subset of which also carry the merB gene), while genes encoding for known Hg methylation (hgcAB) 38 were not detected. The presence of dMeHg and absence of hgcAB has also been reported for marine environments 37 , indicating that other poorly described microorganisms and metabolic pathways may be important for Hg methylation, or that abiotic methylation predominates 31 . Clearly, there is a need to better constrain the activity of Hg cycling microbes in subglacial environments and their impact on the speciation of exported Hg.

Environmentally high mercury concentrations in fjords
Hg concentrations in Greenlandic fjords help elucidate whether inputs from glacial meltwaters have the potential to impact downstream ecosystem health. We find that environmentally high surface dIHg and dMeHg concentrations are sustained downstream of glacial inputs in NK and AF in a region with bedrock geology similar to the Kangerlussuaq field sites (Figs. 1 and 4 and Supplementary Table 1). Surface concentrations of dIHg ranged from 12.5 to 30.6 pM ( Fig. 4a). Lower values were measured at higher salinities, consistent with seawater dilution of a glacial meltwater Hg source. At the salinities sampled in these fjord systems (17.0-32.4), dIHg appears to mix linearly ( Fig. 4a; R 2 = 0.70 (P < 0.05; n = 9)) with lower concentrations at higher salinities implying seawater dilution, but this does not follow a two end-member mixing model using the data from the northern ice sheet catchments sampled. We hypothesize non-conservative behaviour at lower salinities due to removal of Hg by scavenging on flocculated particles, alongside biological and/ or photochemical reduction of Hg 2+ to Hg 0 with subsequent degassing 26 , which is implied by a more limited dataset from SS ( Fig. 1) downstream of RG and LG inputs (see the dHg concentrations from fjord waters with salinity 0.4-8.3 in Extended Data Fig. 5). Dissolved methylmercury (dMeHg) in fjord surface waters accounted for a larger proportion of dHg than meltwater runoff (6-14% versus 0.6-4.7%), and the percentage composition was positively correlated with salinity (R 2 = 0.49; P < 0.05; n = 9), consistent with trends observed in marine environments elsewhere 27 . The dMeHg concentrations were slightly lower than those found in meltwaters (1.56-2.59 pM versus 4.47-16.2 pM), but much higher than open ocean concentrations (0.02-0.13 pM) 27 , indicating non-conservative mixing and dMeHg enhancement in the fjord (Fig. 4b). Additional inputs may include microbial methylation in the water column or fjord sediments, which are likely to be anoxic 39 . The dIHg and MeHg concentrations were more than an  40,41 . These data therefore highlight a role for Greenlandic fjords as Hg hotspots and demonstrate the potential importance of local Hg sources to coastal Arctic ecosystems.

Likely source of mercury to ice sheet meltwaters
The precise origin of Hg in GrIS meltwaters is unknown. The majority of Hg delivered by GrIS meltwater runoff is unlikely to come from anthropogenically derived atmospheric Hg deposited on the ice sheet surface and transported through surface melt because dHg concentrations in meltwater rivers exceed all reported surface snow and ice values (<0.3-150 pM) 18 , including samples from the GrIS (<0.3-10 pM) 42 . Additionally, a large proportion of supraglacial Hg is likely to be photochemically reduced to Hg 0 and volatilized under the high-irradiance conditions 18   a poorly understood mechanism of Hg enrichment, the existence of relatively large quantities of Hg-bearing minerals such as cinnabar in subglacial sediments and/or the influence of overridden Hg-rich permafrost that is being gradually eroded. Overridden permafrost is an unlikely source because pHg concentrations would have to vastly exceed those recorded in recent studies (0.06-0.4 nmol g −1 ) 11 . Instead the unique physicochemical properties of glacial environments are likely to enhance mobilization of Hg in Hg-rich bedrock because of intense bedrock comminution by physical erosion that produces an abundance of freshly weathered microparticles 45 . There is also evidence of geothermal activity under large portions of the GrIS (refs. [46][47][48], potentially manifested locally at LG by an iron-rich groundwater spring located 400 m in front of the glacial terminus 46 with a sHg concentration >20,000 pM. This geothermal activity may be associated with the presence of the geological divide between the Nagssugtoqidian mobile belt to the north and the Archaean block to the south and/or geothermal hotspots under southwestern and central Greenland 46,47 (Extended Data Fig. 1). Data collected from meltwater rivers draining a ~4000 km 2 region of the GrIS and two fjords receiving meltwater runoff indicate a large export of Hg from southwestern Greenland. Extremely high Hg yields, coupled with dMeHg concentrations of environmental concern, highlight the need to better understand ice sheet processes in the Hg cycle, especially the elevated pHg concentrations, which likely fuel the ice sheet Hg cycle. Hypothesized microbial cycling of Hg in the subglacial environment, including unidentified methylation pathways and as a mechanism of protection from the high levels of dissolved inorganic Hg present, provides additional support for the importance of Hg in these systems, but additional evidence is still needed to demonstrate activity of the merA and merB metabolic pathways. Our data indicate that glacial meltwater-sourced Hg impacts downstream fjords, where much lower but environmentally high concentrations persist, with high potential for Hg bioaccumulation in coastal food webs. It is highly uncertain how projected increases in ice sheet melt rates due to climatic warming in the Arctic will impact Hg export; however, we postulate that rising meltwater runoff will increase Hg yields and therefore downstream flux. This large, unaccounted for and climatically sensitive Hg source has not been considered in current global Hg budgets and Hg management strategies, but it should be assessed urgently given the human and economic implications of elevated Hg exposure.

online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41561-021-00753-w.

Site description. Kangerlussuaq (LG, RG and IS).
Glacial meltwater samples were collected from three land-terminating glacial catchments on the southwestern margin of the GrIS during the 2012 (LG), 2015 (LG) and 2018 (LG, RG and IS sampled concurrently) summer ablation seasons. The three glacial catchments are east and northeast of Kangerlussuaq and drain west into the Isortoq River and Watson River (Fig. 1, ED1 and ED6). The hydrology of this sector of the GrIS is relatively well documented 49,51-55 , especially that of LG. All samples were collected from a fast-flowing location at the edge of the meltwater rivers draining the three glacial catchments.
LG extends ~80 km into the ice sheet interior, has a hydrologically active catchment area of ~600-900 km 252,50 and drains through a single portal on the northern side of the glacier snout (Extended Data Fig. 6a). Detailed documentation of the sampling site for LG can be found elsewhere 22 RG is located to the north of LG (Fig. 1) and drains into Akuliarusiarsuup Kuua (Extended Data Fig. 6c). It was sampled in 2018 at 67.104° N, 50.217° W at the same location where hydrochemical gauging took place. The RG catchment is relatively poorly defined and comprises several sub-catchments with multiple outlets. More recent modelling estimates put the active catchment area at ~85 km 2 (including the 'Point 660' sub-catchment defined in ref. 53 ) and the distance it extends inland at ~10-15 km (refs. 50,57 ), IS is located to the north of both the LG and RG catchments. The glacier drains into Isortoq River to the west (Extended Data Fig. 6b NK/AF. NK (Godthåbsfjord) is one of the best-studied fjords in Greenland ( Fig. 1) 14,58-60 . It is ~190 km long, up to 8 km wide, has a mean depth of 260 m and covers an area >2,000 km 2 . Glacial meltwater enters the fjord through a combination of land-and marine-terminating glaciers, of which three are major marine-terminating glaciers (Narssap Sermia, Akullersuup Sermia and Kangiata Nunaata Sermia) and three are land-terminating glacial inputs (Kangilinnguata Sermia, Qamanaarsuup Sermia and Saqqap Sermersua). Freshwater inputs during the summer ablation season drive an estuarine circulation with a low-salinity surface plume flowing toward the fjord mouth to the west 61 . Meltwater input into the fjord is estimated to be ~20 km 3 year −1 , with roughly a 50% contribution from marine-and 50% from land-terminating glaciers, and minor inputs from non-glacial tundra runoff 60 . All glaciers have experienced thinning and accelerating retreat over the past two decades, echoing other outlet glaciers of the GrIS 60 .
AF is located to the south of NK (Fig. 1). Glacial meltwater enters the fjord from a single land-terminating glacier, Kangaasarsuup Sermia. The freshwater input from the river at the head of AF is estimated to be 2.4 km 3 year −1 (1.8-3.7 km 3 year −1 ), with meltwater runoff contributing >75 % of total discharge (mean of Modèle Atmosphérique Régional (MAR)-generated liquid water runoff from 2008 to 2017 (ref. 50

)).
Hydrological monitoring. Hydrological monitoring of LG has been described in detail multiple times elsewhere 49,51,52 , and the same approach was taken for all catchments in this study. In brief, LG was hydrologically gauged at a stable bedrock section throughout the 2012 and 2015 ablation seasons (late April to end of August) and LG, IS and RG during a 3-week period in 2018. A package of hydrochemical sensors connected to a datalogger (Campbell Scientific CR1000 or CR800) were deployed to record pH (Honeywell Durafet), water temperature (Campbell Scientific), electrical conductivity (Campbell Scientific 547 or Keller-Druck DCX-22-CTD, calibrated with an 84 µS cm −1 conductivity standard before deployment), turbidity (Partech C or Turner Designs Cyclops-7 T turbidity sensor) and stage (HOBO pressure transducer) at each sampling site. Stage was converted to discharge using a rating curve determined from rhodamine dye-dilution experiments (n = 41 during 2012, n = 25 during 2015 at LG; and n = 8, 6 and 9 at LG, IS and RG, respectively, during 2018). Uncertainties on discharge measurements were calculated using the root-mean-square error and are estimated as ±10.7-12.1%. Turbidity was converted to suspended sediment concentration (SSC) using a calibration of manual sediment samples collected against the turbidity recorded at the time of sampling 52 . Uncertainties in SSC are estimated to be ±6% (ref. 52 ).
Sample collection, processing and storage. Trace metal clean sampling. All samples were collected according to strict trace metal clean protocols. Sampling bottles (15/250-mL low-density polyethylene (LDPE) bottles and 250/500/1,000-mL fluorinated high-density polyethylene (Fl-HDPE) bottles) and syringes were cleaned sequentially in 1% DECON (~24 h), trace metal grade 3 M HCl (~48 h), followed by trade metal grade 3 M HNO 3 (~48 h) and thoroughly rinsed with ultra-pure water (UPW; Milli-Q, 18.2 MΩ cm −1 ) in between cleaning stages before drying under a laminar flow hood (ISO 5). Whatman GD/XP polyethersulphone (PES) 0.45 µm and Whatman Anotop 25 0.02 µm syringe filters were cleaned with ultra-trace metal grade HCl (Optima T ). The 0.45 µm syringe filters were cleaned by passing through 20 mL of 1.2 M HCl, with the final ~1 mL acidic cleaning solution allowed to sit in the filter for ~2 h before rinsing with 40 mL of UPW and flushing with laminar flow filtered air to dry. The 0.02 µm filters were cleaned by passing 20 mL of 0.02 M HCl, followed immediately by 20 mL of UPW and clean laminar flow filtered air to dry. A portable peristaltic pump or Teflon diaphragm pump and filter capsules were used for sample collection in 2018 (see 'LG/RG/IS 2018' and 'NK/AF 2018'). Peristaltic pump tubing was soaked in 1.2 M trace-metal grade HCl for 24 h and rinsed by pumping several litres of ultra-pure water though it before use. AcroPak 500 capsule filters (0.45 µm with Supor PES membrane) were thoroughly flushed with sample water (several hundred millilitres) before collection of the final sample.
LG/SS 2012. Total particulate Hg was determined on ten samples taken at LG intermittently from 16 May to 11 July 2012. Meltwater was collected in a borosilicate glass bottle, taken to a laboratory tent and immediately filtered through a 0.22 µm membrane filter (PES) with the sediment retained on the filter. Filters were stored air-dried and cool before analysis. Three samples for total dissolved Hg (dHg) were also taken downstream of LG from a salinity gradient in surface waters of SS (ca. 25 km downstream of LG) on 16 June, using the same protocols as particulate Hg samples above, with the 0.22 µm filtrate stored in 250-mL acid-washed borosilicate glass bottles and stored cool in the dark before analysis (Extended Data Fig. 5).
LG 2015. Samples were collected from 1 May to 28 July 2015. Samples from 80 time points were size fractionated (0.02 + 0.45 µm) for total filterable Hg concentrations, and an additional sample was taken from a geothermal spring located 400 m in front of the glacier terminus on 28 July 46 . Bulk water samples were collected in 250-mL or 1000-mL LDPE bottles, zip-lock bagged and taken immediately to a designated 'clean' laboratory tent (no shoes, hairnets and strict cleaning protocols observed) for processing. Bulk samples were filtered through a 0.45 µm syringe filter (12 mL to waste and rinse, with final 10 mL collected), then through a stack of 0.45 µm/0.02 µm syringe filters (12 mL to waste and rinse with final 10 mL collected; see 'Trace metal clean sampling'). These filtered samples were immediately preserved with Optima grade HNO 3 to a final acidity of 0.074 M according to Environmental Protection Agency (EPA) method 245.1 and stored individually double zip-lock bagged in the dark and refrigerated as soon as possible (maximum of 2 weeks).
Although Hg storage in LDPE bottles is not advisable for low-level Hg analysis 62,63 , short-term storage of samples (~weeks) for total Hg analysis with elevated concentrations when acidified is likely to be acceptable 64 . Loss or enrichment of Hg due to transmission through plastic of volatile Hg 0 is the most likely storage artefact 63 due to the porosity of unfluorinated polyethylene. A study has shown a >10 ng L −1 (~50 pM) enrichment in samples stored in LDPE over a period of 5 weeks (refs. 63 ). This corresponds to ~15% of the discharge-weighted mean concentration of sHg and 3.2% of the discharge-weighted mean concentration of cHg from LG in 2015, assuming a worst-case scenario of 200 pM contamination over the storage to measurement period of ~120 days. Furthermore, any contamination effect should be reflected in the procedural blanks.
LG/RG/IS 2018. Samples were collected concurrently at three glaciers over a 25-day period in 2018 from 19 June to 14 July to determine dissolved inorganic Hg (dIHg, as Hg 2+ ) and dissolved methylmercury (dMeHg) concentrations. Sample collection took place at 11:00 at LG and RG, and 15:00 at IS on the same days, the latter delayed by 4 h to account for the water transit time from the glacier terminus to the sampling site ~26 km downstream (that is, to capture a body of water similar to the two other sampling locations; Fig. 1). Bulk water samples were collected in acid-cleaned 1000-mL Fl-HDPE bottles and taken immediately for processing in a laboratory tent or field laboratory. Water was filtered through a 0.45 µm AcroPak 500 capsule filter using a peristaltic pump and acid-cleaned C-Flex tubing into 250-mL Fl-HDPE bottles after rinsing with three aliquots of filtered water. Samples were immediately acidified using ultra trace metal grade (Optima) HCl to a final concentration of 0.048 M (0.4% v/v concentrated HCl). Sample bottles were double zip-lock bagged and stored in the dark until returning to the University of Bristol, where they were transferred to refrigerated storage.
NK/AF 2018. Samples were collected from the surface waters of two fjord systems during a cruise in July 2018 (5 July-9 July). Trace metal clean sampling was undertaken using a towfish system coupled to acid-cleaned polyvinylchloride tubing, a Teflon diaphragm pump (A2CH1 F8 AstiPure II, Saint Gobain) and final filtration through a 0.45 µm AcroPak 500 capsule filter. The towfish comprised a plastic fin held ~1.5 m underwater by ~25 kg of epoxy-covered dive weights. Samples were filtered into 500-mL Fl-HDPE bottles after rinsing three times with filtered sample water, then acidified using ultra trace metal grade (Optima) HCl to a final concentration of 0.048 M (0.4% v/v concentrated HCl) and stored refrigerated until analysis. Geochemical analysis. Particulate and total dissolved mercury (2012). Particulate and total dissolved mercury concentrations were determined at the Woods Hole Oceanographic Institution. Filters were first oven dried at 55 °C overnight before being digested in ultra trace metal grade 2 M HNO 3 (refs. 65,66 .) Briefly, filters were placed in pre-cleaned 15-mL polypropylene centrifuge tubes, and 8 mL of acid was added. Tubes were inverted to mix then placed in an ultrasonic bath at 55 °C for 15 min. Samples were then centrifuged for 5 min at 3,000 rpm, and supernatant retained in a new, clean tube for total Hg analysis.
Total Hg in leach solutions and in filtered fjord samples was determined by cold vapour atomic fluorescence spectroscopy (CV-AFS) following EPA method 1631. Samples were digested with bromine monochloride (BrCl) and then reduced with hydroxylamine hydrochloride (NH 2 OH⋅HCl). Hg present was reduced to Hg 0 with stannous chloride (SnCl 2 ) and quantified by dual Au-amalgamation CV-AFS using an external standard calibration curve. Filter blanks processed as per samples were all below the instrumental limit of detection (<0.5 pM).
Total dissolved mercury (2015). Total dissolved Hg in 0.02 and 0.45 µm filtered samples were determined using a Thermo Scientific XSERIES 2 quadrupole ICP-MS with collision/reaction cell after acidification of samples to 1% (v/v) HCl (laboratory distilled ultra trace metal grade) at the National Oceanography Centre in Southampton, within 3 months of sample collection. This method is sub-optimal for analysis of Hg due to a memory effect and very high first ionization potential in ICP-MS analysis, leading to low sensitivity and prolonged sample washout. Several measures of analytical quality were noted to ensure reliable results. For example, the analytical precision over the concentration range observed in samples was ±23.4%. Sample memory effects (a common issue in Hg ICP-MS analysis), determined by running a high-concentration sample followed by five analytical blanks, was high at ~400 pM in the first blank to <100 pM in the final blank. Procedural blanks (n = 16) taken in the field using transported ultra-pure water and the same procedure as samples, run intermittently throughout the run, averaged 522 ± 255 pM, and the mean was subtracted from measured values before reporting. It is expected that a large fraction, if not all, of this (~400 pM), is due to the system memory effect or storage contamination in LDPE bottles, further highlighting that this methodology is not suitable for low-level Hg analysis. The detection limit (mean of ten blanks + standard deviation × 3) was therefore high at 800 pM, and all blanks fell under this. A subset of samples were also analysed by CV-AFS (as per section 4.1) in a separate Hg laboratory (Greg Carling at BYU) with blank values <200 pM and similarly high sample concentrations that correlated to those reported via ICP-MS (R 2 = 0.58) with a positive offset for CV-AFS samples (~1,000 pM), despite 4 years of storage time, which would otherwise be deemed unacceptable. The higher values reported for the CV-AFS method are likely either due to underreporting by ICP-MS, or sample contamination over the prolonged storage period 63  and MeHg spike recoveries of 103% (77-111%). Precision values based on replicate analysis (n = 5) of 10 pM Hg 2+ standards and 2 pM MeHg standards placed randomly throughout the analytical run were ±3.1% and ±9.1% respectively. Blanks were below instrumental detection limits for both Hg 2+ (~<1 pM) and MeHg (~<0.1 pM).
Fjord samples were analysed using techniques similar to those above. MeHg and Hg 2+ concentrations were determined following Mansfield and Black 68 with addition of NaCl before direct aqueous ethylation using NaBEt 4 . Samples were diluted 1:2 in Hg-free UPW (Milli-Q, 18.2 MΩ cm −1 ) and buffered to pH 4.1 (using 3 M sodium acetate/acetic acid buffer), followed by NaCl addition to a final concentration of 1.2 M and direct ethylation with addition of NaBEt 4 to a final concentration of 40 µM. Spiked samples (n = 12, spike values: Hg 2+ concentration of 10 pM, MeHg concentration of 1 pM) displayed mean Hg 2+ recovery values of 117% (78-153%) and MeHg recoveries of 104% (88-119%). MeHg precision was ±8.3%, calculated from five replicate standards (1 pM) run randomly throughout the analytical run, while Hg 2+ precision using the direct ethylation technique was ±11.7%. Three NaCl blanks (that is, blanks treated exactly the same as samples with 1.2 M NaCl addition) were below the detection limit for MeHg (<0.1 pM), and were ~1 pM for Hg 2+ . Total dissolved Hg was also determined independently in two other Hg laboratories using CV-AFS and EPA method 1631 (Greg Carling and Carl Lamborg) as in 'Particulate and total dissolved mercury (2012)' . Differences between the total dissolved Hg concentration in the three laboratories were ±23% (±8-42%). Mean total dissolved concentration from all three measurements are used as the reported value for TdHg (Fig. 4) along with individual error based on the standard deviation of the three concentrations reported by individual laboratories.
Dissolved organic carbon (2018). Samples for DOC concentrations were filtered using 0.45 µm Whatman Aqua 30 syringe filters into acid-washed HDPE bottles (10% HCl v/v for 48 h) and kept frozen in the dark until analysis. DOC concentrations were measured after acidification using a Shimadzu TOC-L CHN analyser with high-sensitivity catalyst. Sample concentrations were calculated from a manual calibration curve based on standards (0.05-2 mg C L −1 ), gravimetrically prepared from a certified TOC standard of potassium hydrogen naphthalate (1,000 ± 10 mg C L −1 ; Sigma Aldrich). The limit of quantification (LOQ) was 0.026 mg L −1 (where LOQ = LOB + 5 × s.d. of low-concentration sample, where LOB is the limit of baseline) 69 . Precision was calculated from repeat standard measurements of 0.10 and 0.50 mg C L −1 and was better than ±5%.

Microbiological analysis (2015)
. Samples for microbiological analysis were taken at the portal of LG during May-July 2015 concurrently with samples taken for geochemical analysis to identify genes associated with Hg cycling in the subglacial microbial community. Water samples were passed through a Sterivex syringe filter (0.22 µm; Millipore) using a sterile syringe until clogged (~300 mL). Water was flushed from the filter chamber, filled with nucleic acid preservation buffer (LifeGuard, QIAGEN) and immediately frozen in the field at −20 °C. Total DNA in suspended sediment retained in the Sterivex filter was extracted using the PowerWater Sterivex DNA isolation kit (MO BIO). A full description of microbiological sample collection and DNA extraction can be found in refs. 22,70 , as the extracts sequenced for metagenomes here consisted of subsamples of those used for 16S rRNA gene analysis in the cited studies.
Four pooled samples were used for metagenomic sequencing due to the low DNA content of individual extracts and cost. Pool 1 represented early-season assemblages (three daily samples in the period from 4 May to 7 June), pool 2 represented 'early' outburst events (five daily samples in the period from 20 June to 2 July), pool 3 represented 'late' outburst events (four daily samples in the period from 5 to 10 July) and pool 4 represented late season (four daily samples in the period from 13 July onward) when an efficient drainage system was established. A total of at least 20 ng of treated DNA from each pooled sample served as the input for the TruSeq Nano DNA library preparation kit (Illumina). Insert length ranged between 295 and 409 bp. DNA libraries were prepared using a TruSeq Nano DNA LT kit and sequenced (2 × 150 bp) in four lanes of Illumina NextSeq500 using a High-Output V2-300 kit at the Bristol Genomics Facility. The four metagenomes obtained were assembled as follows: reads were trimmed by removing adapters with Trimmomatic using Illumina TruSeq2-PE adapters with a seed mismatch threshold, palindrom clip threshold and simple clip threshold set at 2, 30 and 10, respectively. Sequencing reads were quality filtered by base call quality using the FASTX-Toolkit, specifically fastq_quality_ filter, with the following parameters: -Q33 -q 30 -p 50. The paired-end assembly of the remaining reads was performed with the MEGAHIT assembler (v. 1.1.3) with the following parameters: megahit −1 all.pe.qc.fq.1 -2 all.pe.qc.fq.2 -r all.se.qc.fq -o all.Megahit.assembly. Reads from all samples were assembled together. A total of 9,168,354 contigs were obtained, with an average length of 618.7 bp; N50 was 712 bp. The total length of the assembled metagenome reached 5,672,583,505 bp, while the longest contig had 297,443 bp. Assembled data were annotated using the MG-RAST automated pipeline against the RefSeq database. For phylogenetic analyses, homologues of translated merA sequences were found using BLASTP against the National Center for Biotechnology Information nr database 71 . Representatives of known main lineages were also added. Sequences were aligned in MAFFT v. 7 72 and trimmed manually. The phylogenetic tree (Extended Data Fig. 4) was determined by Phyml v. 3.1 73 under the LGgamma model and non-parametric bootstrap analysis with 100 replicates.
Statistical analyses. All statistical analyses were performed in JMP Pro 15 (SAS). Student's t test was used to determine whether population means were statistically different. Linear correlations were used to determine the strength of relationship between two variables and the significance of the correlation. P values are presented in the main text with P < 0.05 assumed to be statistically significant.
Flux estimations. Discharge-weighted mean method. The first-order flux estimates presented in Table 1 are calculated using the discharge-weighted mean concentrations from the three 2018 sampling sites, where the most complete IHg and MeHg record is present. The discharge-weighted concentration was then multiplied by the mean annual modelled meltwater discharge from the southwestern margin of the GrIS (2009-2018; Extended Data Fig. 2) from ref. 50 to generate dIHg and dMeHg flux estimates for southwest Greenland.
LOADEST method. Catchment fluxes were also estimated using the constituent-load modelling software LOADEST 74 , with measured 2018 mercury concentrations and mean daily discharge. Regression models for concentrations were generated using the adjusted maximum-likelihood estimator and automatic best model selection (model 0) to allow for model selection based on the Akaike information criterion. Total annual fluxes were estimated using the mean concentration output of the model with the modelled meltwater discharge detailed in Table 1, with uncertainty defined as the modelled concentrations at the 25th and 75th percentiles. These data are presented in Supplementary Table 2. Mean flux estimates deviated ±0.7-7.6% from those calculated using the discharge-weighted mean concentrations in Table 1.
Representativeness of the study area. There are obvious limitations to using three glaciers to estimate Hg export from the southwestern margin of the GrIS. We acknowledge that IS, LG and RG collectively cover a small proportion of the ice sheet (~0.3% by surface area) and southwestern Greenland (~4.1%). However, the area is significantly larger (by almost two orders of magnitude) than any other glacier study area investigated thus far in the literature (both in Greenland and worldwide), and to our knowledge, this is the only study to date that has measured Hg concentrations in ice sheet runoff (rather than disconnected ice caps in Greenland). Furthermore, the annual discharge of these three outlets equates to nearly 8% of the runoff from the southwestern sector of the GrIS. The Kangerlussuaq area is an excellent analogue study area for the southwestern region of the GrIS for several other reasons that have been documented in previous studies, including the underlying geology, hydrology and because its discharge is proportional to total modelled ice sheet runoff 49,51,[75][76][77][78] . The underlying debris and morphology of the catchment are similar to other catchments of southwestern GrIS 79 (Extended Data Fig. 2), and bedrock geology is predominantly Neoproterozoic gneiss/granite, which is typical of large areas of the crystalline rocks that dominate the Precambrian Shield on which Greenland lies 75,80 (Extended Data Fig. 1). The catchment hydrology here is also believed typical of the large Greenland outlet glaciers that dominate discharge of meltwaters from the GrIS 76-78,81 .

Data availability
All data contained in this article that are not included in tables are freely available (https://doi.org/10.5281/zenodo.4634280). The metagenomic data are accessible in MG-RAST under the number 4765056.3.