Proglacial groundwater spring sampling and analysis
Field samples for groundwater methane concentrations were collected between February and May of both 2021 and 2022. Sampling locations are detailed in Fig. 3a and the Supplementary Data and represent winter icings of 78 glaciers across central Svalbard. Proglacial groundwater springs were reached by drilling into ice blisters (large, upheaved dome features of ice formed by pressurized water flow beneath or within the icing; Fig. 2b) on winter icings with a 7-cm-diameter auger until pressurized water was released. Unfiltered samples for methane analysis were taken directly from the spring outlet in 22 ml gas-tight serum vials, ensuring they were bubble-free before sealing shut. Samples were spiked within 24 h with 1 ml of 1 M NaOH while concurrently removing 1 ml of sample to stop microbial activity. Samples were stored upside down in the dark at 4 °C until analysis.
A headspace of 4 ml pure N 2 gas was created in each vial, and samples were left to equilibrate at room temperature for at least 24 h. The concentration of methane in the headspace was measured at Queen Mary University of London by gas chromatography fitted with a flame ionization detector (Agilent Technologies UK Ltd). Calibration gas standards (100 ppm CH 4 ; BOC Limited) and ambient air were analysed in each analytical sequence and repeated at regular intervals to check for drift. The total amount of methane in each sample (headspace and water phase) was calculated using the ideal gas law and solubility coefficients from ref. 50 (Supplementary Information). Stable carbon isotope ratios of methane (δ13C–CH 4 ) in samples with sufficiently high methane concentrations were analysed at the University of Cambridge in the LASER-ENVI facility using a cavity ringdown spectrometer (Picarro G2201-I, Picarro Inc.) and reported relative to the Vienna Pee Dee Belemnite standard.
Icing depth measurements
Radio echo sounding with GPR equipment was used to obtain icing thickness measurements, which were further validated by manual auger drillings. The GPR surveys were performed using a Malå ProEx control unit, 800 MHz, 500 MHz and 100 MHz antennas (changed due to equipment availability over the field season) and a generic USB (universal serial bus) GNSS (global navigation satellite system) receiver. A total of 29 km of lines was surveyed over 11 discrete proglacial icings during spring 2022, and 25 validation points along the GPR track were drilled manually. Post-processing of the GPR data included time-zero correction, horizontal and vertical resolution standardization, topographic Kirchhoff migration51 with a medium velocity of clean ice (0.168 m ns–1) and automatic gain control. The correlation between the digitized radargrams and the manual auger drillings was significant (R² = 0.74), with no apparent bias, which validates that the correct reflection was digitized and the medium velocity of clean ice is accurate.
Methane emission calculations
Icing volume calculations
Icing thickness measurements obtained by GPR surveys were used to determine an average ice depth across all icings. Areas of each glacial icing were measured on Sentinel-2 satellite images from May 2021 and April 2022 using the polygon feature in the geographic information system application QGIS. Measurements were repeated three times on separate days to determine an error for the measurement. Icing areas were multiplied by the average icing depth to obtain an icing volume at each glacial site.
Monte Carlo simulation
A Monte Carlo simulation, run 10,000 times for each icing, was used to calculate a range of potential emissions from each glacial site. The Monte Carlo sampling was taken from a normal distribution according to the defined errors of each parameter (Supplementary Data). Maximum and minimum emission scenarios were considered for each glacier, in which the maximum measured methane concentration and icing area were used for the highest emission scenario and the minimum concentration and icing area were used in the lowest emission scenario. These emissions were summed to determine maximum and minimum total methane fluxes from the full study region.
Annual methane emissions for each glacial icing were estimated using the following equation:
$${f}_{{{{\mathrm{CH}}}}_{4},i}=\frac{{A}_{{{\mathrm{icing}}},i} \times {d}_{{{\mathrm{icing}}}({{\mathrm{avg}}})}\times {\rho }_{{{\mathrm{ice}}}}\times \left({[{{{\mathrm{CH}}}}_{4}]}_{i}-{[{{{\mathrm{CH}}}}_{4}]}_{({{\mathrm{eq}}})}\right)}{{n}_{{{\mathrm{days}}}{{\mathrm{of}}}{{\mathrm{icing}}}{{\mathrm{formation}}},i}}\times 365.24$$
where f CH4,i is the annual flux of methane from the glacial icing i, A icing,i is the measured area of the icing i, d icing(avg) is the average thickness of all icings determined by GPR, ρ ice is the density of ice (0.917 g cm–3), [CH 4 ] i is the aqueous methane concentration of the formation water of the icing i, [CH 4 ] (eq) is the equilibrium concentration of methane in water at 0 °C and 1 atm (4 nM) and n days of icing formation,i is the number of days between the onset of constant freezing air temperatures in the previous year to the date of the satellite image used to calculate the icing area.
Assumptions made for the Monte Carlo simulation
In the calculations, we assumed that the total discharge from the proglacial groundwater spring since the onset of freezing temperatures has been trapped frozen in the icing, and thus the icing volume can be used as a proxy of total discharge over this period. We also made the conservative assumption in our calculations that this discharge rate is constant year-round. Finally, it is assumed that the aqueous methane concentrations of the proglacial springs are constant throughout the year. Limited data collected in summer from a subset of the sites support this simplification (Supplementary Data). All assumptions made are listed in the following, and all errors considered in the calculations are provided in the Supplementary Data.
Groundwater spring discharge rate is constant throughout the year.
All methane above the ~4 nM equilibrium concentration will be emitted to the atmosphere.
Measured concentrations represent the methane concentrations throughout the year for each spring sampled.
Measured methane concentrations represent the concentrations of all groundwater that has formed the icing.
All icings have a similar average depth, within the range of error.
All discharge from the groundwater spring since the onset of freezing until the date of area measurement has been trapped within the icing.
Any snow incorporated within the icing is negligible.
Any saturated sediments below the icing are negligible in terms of discharge water volume.
Samples taken from the same location but on different days are from the same groundwater source, and their methane concentrations can be averaged.
Samples taken from different locations within the icing represent different groundwater springs if their water chemistries are different.
Water samples have not been freeze-concentrated or diluted by inclusion of snow around the icing.
The addition of NaOH to each sample has ceased any microbial activity that may otherwise have altered the methane concentrations during sample storage.
Satellite imagery captures the full icing.
Density of all icing ice is 0.917 g cm −3 .
All glacial icings within the selected region have been identified and are represented in the calculations.
The entirety of the previous year’s icing melted before the new icing began to form (satellite imagery confirms that this is largely the case).
Handling of icings that were not sampled
Twenty-two glacial icings that were not sampled were identified in the sampling region by satellite imagery (Fig. 3a). These icings were included in the Monte Carlo simulation to consider emissions from all proglacial groundwater springs in the study region. Icing volumes of these sites were determined following the preceding approach. An average methane concentration from all sampled proglacial groundwater springs was used in place of a measured methane concentration for each of the unsampled sites to calculate a methane flux and estimate their potential contribution to a regional methane flux.
Limitations in the method
Our method for calculating methane emissions is likely to underestimate the discharge rates of the groundwater springs due to challenges in assessing the icing volume. As satellite imagery is used to measure the icing area, it is not possible to determine whether sections are snow covered, and therefore large areas of the icings may be left out of the measurements. In addition, the GPR surveys used to determine an average icing thickness were conducted between late March and early May, before growth of the icing was complete. This means that the measurements do not represent the maximum depth reached by the icing during the winter season and therefore underestimate total discharge.
Sampling and analysis of fjord surface waters
Sampling of fjord surface waters adjacent to marine-terminating glaciers was completed in April and May of 2021 while the fjords were still covered by sea ice. Samples were taken in a high-spatial-resolution grid pattern in front of the glaciers and were acquired by drilling holes through the sea ice and using a bailer to retrieve surface water from beneath the ice. Samples for methane concentration analysis were transferred into 100 ml gas-tight serum vials using a syringe, ensuring they were bubble-free before sealing. Electrical conductivity measurements were taken in the field with a Hach probe. Upon return from the field, samples for methane concentration measurements were spiked with 1 ml of 1 M NaOH and stored upside down in the dark at 4 °C until analysis. Before analysis, a headspace of 5 ml N 2 was created, replacing 5 ml of sample, and left to equilibrate for more than 24 h at room temperature. Measurements of methane in the headspace were made at the University of Tromsø with a ThermoScientific GC Trace gas chromatograph fitted with a flame ionization detector and an MSieve 5 A column (ThermoScientific). Calculations of methane concentrations were handled as described in the Supplementary Information, while also accounting for salinity.
Statistical analysis