Heterogeneous CO₂ and CH₄ content of glacial meltwater from the Greenland Ice Sheet and implications for subglacial carbon processes

Accelerated melting of the Greenland Ice Sheet has increased freshwater delivery to the Arctic Ocean and amplified the need to understand the impact of Greenland Ice Sheet meltwater on Arctic greenhouse gas budgets. We evaluate subglacial discharge from the Greenland Ice Sheet for carbon dioxide (CO₂) and methane (CH₄) concentrations and δ¹³C values and use geochemical models to evaluate subglacial CH₄ and CO₂ sources and sinks. We compare discharge from southwest (a sub-catchment of the Isunnguata Glacier, sub-Isunnguata, and the Russell Glacier) and southern Greenland (Kiattut Sermiat). Meltwater CH₄ concentrations vary by orders of magnitude between sites and are saturated with respect to atmospheric concentrations at Kiattut Sermiat. In contrast, meltwaters from southwest sites are supersaturated, even though oxidation reduces CH₄ concentrations by up to 50 % during periods of low discharge. CO₂ concentrations range from supersaturated at sub-Isunnguata to undersaturated at Kiattut Sermiat. CO₂ is consumed by mineral weathering throughout the melt season at all sites; however, differences in the magnitude of subglacial CO₂ sources result in meltwaters that are either sources or sinks of atmospheric CO₂. At the sub-Isunnguata site, the predominant source of CO₂ is organic matter (OM) remineralization. However, multiple or heterogeneous subglacial CO₂ sources maintain atmospheric CO₂ concentrations at Russell but not at Kiattut Sermiat, where CO₂ is undersaturated. These results highlight a previously unrecognized degree of heterogeneity in greenhouse gas dynamics under the Greenland Ice Sheet. Future work should constrain the extent and controls of heterogeneity to improve our understanding of the impact of Greenland Ice Sheet melt on Arctic greenhouse gas budgets, as well as the role of continental ice sheets in greenhouse gas variations over glacial–interglacial timescales.

1 Introduction

Glaciers play an important role in global chemical cycles due to the production of fine-grained sediments that participate in carbonate and silicate mineral weathering reactions (Table 1), which are the principal sink of atmospheric CO₂ over geologic timescales (Berner et al., 1983; Walker et al., 1981). Variations in the weathering intensity of comminuted sediments may contribute to glacial–interglacial atmospheric CO₂ variations as sediments are alternately covered by ice and exposed following ice retreat. However, the importance of CO₂ consumption by mineral weathering is poorly understood, including effects from the advance and retreat of continental ice sheets (Ludwig et al., 1999). Recent evaluations of carbon budgets in proglacial environments indicate that mineral weathering results in net sequestration of atmospheric CO₂, suggesting that proglacial systems are underrecognized as Arctic CO₂ sinks (St Pierre et al., 2019); however, alternate processes could lead to the production of greenhouse gases in glacial systems. For instance, CH₄ production in anaerobic subglacial environments driven by the remineralization of organic matter (OM) contained in soils and forests covered during glacial margin fluctuations has been suggested as a potential carbon feedback to drive warming (Sharp et al., 1999; Wadham et al., 2008). Because the global warming potential of CH₄ is 25 times greater than CO₂, even limited subglacial methanogenesis has the potential to strongly impact the greenhouse gas composition of glacial meltwater. Combined inorganic and organic subglacial processes may therefore produce glacial meltwater that is a source or sink of greenhouse gas. While the net impact of these processes on modern carbon fluxes is poorly constrained, determining these impacts will improve modern carbon budgets as well as depictions of how fluxes may have evolved during the advance and retreat of continental ice sheets.

Table 1 Mineral weathering reactions and impacts on dissolved CO₂ concentrations.

^* Abbreviations are based first on the mineral class (carbonate: carb; silicate: sil) and then on the acid (carbonic acid: CA; sulfuric acid: SA).

In subglacial environments where remineralization is limited by low OM availability, the major element solute load of glacial meltwater is typically dominated by products of mineral weathering reactions (Tranter, 2005). The extent of mineral weathering in subglacial environments depends in part on the availability of acids to drive reactions, namely sulfuric and carbonic acids (Table 1). Sulfuric acid is derived from the oxidation of reduced sulfur species, which largely occur as iron sulfide minerals including pyrite (Tranter, 2005). Sulfide oxidation may occur abiotically; however, the kinetics of microbially mediated sulfide oxidation is several orders of magnitude faster and may lead to local depletion of oxygen given a sufficient supply of sulfide minerals (Sharp et al., 1999). In contrast, carbonic acid may be derived from multiple external or in situ sources of CO₂ to the system. The dominant external source is supraglacial meltwater that flows to the subglacial system through moulins following equilibration with atmospheric CO₂ (Fig. 1). Unlike proglacial environments where free exchange between water and the atmosphere may resupply CO₂ consumed by weathering, subglacial environments may be partially or fully isolated from the atmosphere, limiting further atmospheric CO₂ invasion and thus the extent of mineral weathering with carbonic acid. However, additional atmospheric CO₂ may be delivered in open portions of the subglacial environment though exchange in fractures or moulins along subglacial flow paths or in partially air filled conduits, allowing a much greater magnitude of carbonic acid mineral weathering (Graly et al., 2017). CO₂ may also be derived from in situ sources, such as gaseous CO₂ contained in ice bubbles of basal ice or fluid inclusions in rocks that release volatiles (including CO₂) following mechanical grinding (Macdonald et al., 2018). When OM is available, its remineralization also generates CO₂ (and potentially CH₄) along with nutrients, but low OM availability in many subglacial systems limits remineralization as a CO₂ source (Fig. 1).

Figure 1 Conceptual diagram of subglacial sources and sinks of CO₂ and CH₄. Arrows indicate the direction of fluxes. Boxes represent processes, and sources of gases to subglacial meltwaters are indicated by blue text, while sinks of gases to subglacial meltwater are indicated by red text. Gas bubbles, mechanical grinding, and OM remineralization are grouped because all are CO₂ and CH₄ sources.

The role of subglacial carbon processes may play an increasingly important role in modern Arctic carbon budgets as disproportionate warming increases glacial meltwater and sediment fluxes to the ocean, particularly from the Greenland Ice Sheet (Wadham et al., 2019). The is the last remaining ice sheet in the Northern Hemisphere following collapse of other ice sheets since the Last Glacial Maximum (∼ 20 ka). It has been losing mass at increasing rates that averaged 286 ± 20 Gt/yr between 2010–2018, representing a 6-fold increase since the 1980s (Mouginot et al., 2019). While mineral weathering significantly modifies the chemical composition of Greenland Ice Sheet subglacial discharge (e.g., Hindshaw et al., 2014; Deuerling et al., 2018; Urra et al., 2019) and should consume CO₂ similar to other glacial and proglacial environments, the recent identification of microbially driven reactions (including methanogenesis) in subglacial environments of the Greenland Ice Sheet indicates that organic processes may also play a role (Dieser et al., 2014; Lamarche-Gagnon et al., 2019). The relative importance of subglacial greenhouse gas sinks (CO₂ consumption through mineral weathering) and sources (such as OM remineralization) determines the greenhouse gas composition of subglacial discharge, which may then serve as a source or a sink of atmospheric greenhouse gases. Constraining the relative impacts and variability of these processes underneath the Greenland Ice Sheet will provide important information regarding the current and future impact of Greenland Ice Sheet loss on Arctic carbon budgets, as well the role of continental ice sheets on carbon cycle feedbacks.

To evaluate the net impact of carbon processes on the greenhouse gas composition of subglacial discharge of the Greenland Ice Sheet, we compare water chemistry, dissolved CO₂ and CH₄ concentrations, and gas stable isotopic compositions between three subglacial discharge sites draining land-terminating glaciers of the Greenland Ice Sheet over the melt seasons of 2017 and 2018 (Fig. 2). We employ mass balance models utilizing the concentrations of major cations and anions to determine the magnitude of the impact on CO₂ concentrations from mineral weathering reactions (Table 1). These results are combined with measured gas concentrations and δ¹³C to determine the relative importance of mineral weathering compared to OM remineralization on the CH₄ and CO₂ content of subglacial discharge. We also assess the temporal and spatial variability of these processes under the Greenland Ice Sheet to improve our understanding of carbon cycling in Greenland subglacial environments and the implications of Greenland Ice Sheet mass loss on Arctic carbon budgets.

Figure 2 Google Earth satellite images of study locations in (a) Greenland including (b) locations near the town of Kangerlussuaq, including the sub-Isunnguata water sampling location (dark blue circle) and Russell water sampling location (light blue square). The gold star represents the location of AK4, where proglacial river discharge records were collected. (c) Location of Kiattut Sermiat site (orange triangle) near the town of Narsarsuaq in southern Greenland where water samples were collected. © Google Earth.

2 Methods

2.1 Study locations

Our three subglacial discharge locations are located in southwest (Fig. 2a, b) and southern (Fig. 2a, c) Greenland (pictures given in the Supplement). Our sub-Isunnguata watershed (IS; $\mathrm{67}{}^{\circ }{\mathrm{09}}^{\prime }{\mathrm{27.1}}^{\prime \prime }$ N, $\mathrm{50}{}^{\circ }{\mathrm{03}}^{\prime }{\mathrm{25.0}}^{\prime \prime }$ W) and Russell Glacier watershed (RU; $\mathrm{67}{}^{\circ }{\mathrm{05}}^{\prime }{\mathrm{22.1}}^{\prime \prime }$ N, $\mathrm{50}{}^{\circ }{\mathrm{14}}^{\prime }{\mathrm{18.7}}^{\prime \prime }$ W) drain to the Akuliarusiarsuup Kuua, which is a tributary to the Qinnguata Kuussua. The short stretch of river downstream of the confluence of the Akuliarusiarsuup Kuua and the main stem of the Qinnguata Kuussua near the town of Kangerlussuaq is also known as the Watson River (Fig. 2b). The majority of drainage from the Isunnguata Glacier drains to the northern Isortoq River (Fig. 2); however, a sub-catchment of the Isunnguata drains to a stream that directly feeds the Akuliarusiarsuup Kua, which we refer to as the Northern Tributary, while the Akuliarusiarsuup Kuua refers to river that flows between the outlet of the Russell Glacier and the confluence with the Qinnguata Kuussua (Fig. 2b). Watson River discharge is monitored by PROMICE (Programme for Monitoring of the Greenland Ice Sheet; van As et al., 2018), and total discharge was 4.3 and 3.6 km³ of water in 2017 and 2018, respectively (van As et al., 2018). The total catchment size for the Isunnguata is 15 900 km², though the size of the sub-catchment draining to the Northern Tributary is much smaller with a drainage area of approximately 40 km² (Lindbäck et al., 2015; Rennermalm et al., 2013); therefore, we hereby refer to this site as the sub-Isunnguata watershed. The total drainage area for the Russell Glacier is not precisely known; however, the catchment draining both the Russell and Leverett glaciers has been estimated at approximately 900 km² (Lindbäck et al., 2015), while the Leverett drainage area alone is estimated at approximately 600 km² (Hawkings et al., 2016). We therefore estimate the Russell drainage area at approximately 300 km², though it may be considerably smaller (van de Wal and Russell, 1994). Discharge from the third site in southern Greenland, Kiattut Sermiat (KS; $\mathrm{61}{}^{\circ }{\mathrm{12}}^{\prime }{\mathrm{13.5}}^{\prime \prime }$ N, $\mathrm{45}{}^{\circ }{\mathrm{19}}^{\prime }{\mathrm{49.1}}^{\prime \prime }$ W), drains to the Kuusuaq River near the town of Narsarsuaq. While Kuusuaq discharge is not monitored, a previous study using dye tracing techniques estimated approximately 0.22 km³ of discharge in 2013, and its catchment size was estimated at 36 km² (Hawkings et al., 2016).

Underlying lithologies differ between sites. Southwest Qinnguata Kuussua sites are located near the boundary between the Archean craton to the south and the southern Nagssugtoqidian Orogen to the north (Henriksen et al., 2009). The Archean block is composed of granites and granulite facies orthogneisses that were intruded by mafic dikes during Paleoproterozoic rifting. These rocks were deformed and modified during subsequent continent-to-continent collision in the Paleoproterozoic to create the amphibolite facies gneisses of the southern Nagssugtoqidian Orogen (van Gool et al., 2002). Kiattut Sermiat lies within the Paleoproterozoic Ketilidian fold belt (Henriksen et al., 2009). Lithologies in this region include the Julianehåb granite and associated basic intrusions and the sedimentary and volcanic rocks of the Mesoproterozoic Gardar Province that include a suite of alkaline igneous rocks and basaltic dikes with interbedded sandstones (Kalsbeek and Taylor, 1985; Upton et al., 2003).

Previous studies have characterized chemical weathering reactions in subglacial discharge to the Akuliarusiarsuup Kua and Qinnguata Kuussua (Deuerling et al., 2018; Hasholt et al., 2018; Yde et al., 2014), in the Kuusuaq that drains Kiattut Sermiat (Hawkings et al., 2016), and in comparison between these two sites (Urra et al., 2019). There has been extensive work regarding ice sheet dynamics and hydrology in the southeast glaciers draining to the Akuliarusiarsuup Kua and Qinnguata Kuussua (Van As et al., 2017, 2018; Lindbäck et al., 2015) as well as southern Kuusuaq catchments (Warren and Glasser, 1992; Winsor et al., 2014). Subglacial permafrost has been identified near the sub-Isunnguata site (Ruskeeniemi et al., 2018) and mostly likely formed during Holocene fluctuations in the ice sheet margin. While a similar Holocene ice retreat and re-advance may have occurred in southern Greenland (Larsen et al., 2016), it is unknown whether this retreat led to the organic deposits below the ice sheet.

2.2 Sample collection

We collected water samples from subglacial discharge sites in spring and fall of 2017 and the summer of 2018 to observe seasonal variations in water chemistry. To minimize atmospheric influence, samples were collected as close as possible to the glacier front where subglacial waters emerge, which was less than 10 m for the sub-Isunnguata site and approximately 100 m for the Russell Glacier site. Subglacial discharge from Kiattut Sermiat site flowed through a glacial meltwater lake prior to arriving at the sampling location, which was approximately 1.1 km from the glacial outlet (Fig. 2) and therefore may experience a more interaction with atmospheric gases during transport from the subglacial discharge site to our sampling location. We collected water samples by pumping water through a 0.5 cm flexible PVC tube that was placed in flowing water as far as possible from shore (approximately 1–2 m). A YSI Pro Plus sensor that was calibrated daily was installed in an overflow cup filled from the bottom to measure specific conductivity (Sp.C), temperature, pH, dissolved oxygen, and oxidation–reduction potential (ORP). These parameters were monitored until stable, between about 10 and 30 min, after which samples were collected and preserved in the field according to the solute to be measured after being filtered through a 0.45 µm trace-metal grade Geotech high-capacity disposable canister filter. Samples for cations and anions were collected in HDPE bottles; cation samples were preserved with Optima-grade ultrapure nitric acid (pH < 2), while no preservative was added to anion samples. Dissolved inorganic carbon (DIC) samples were filtered through 0.2 µm filters directly to the bottom of 20 mL Qorpak glass vials and allowed to overflow until sealed tightly with no headspace.

Gas samples were collected in duplicate via headspace extractions according to methods outlined in Repo et al. (2007) and Pain et al. (2019b). Unfiltered water was pumped into the bottom of 500 mL bottles until they overflowed. Bottles were immediately capped with rubber stoppers fitted with two three-way inlet valves. A total of 60 mL of water was extracted from one inlet and replaced with 60 mL of atmospheric air (for spring and fall 2017 sampling trips) or ultrapure N₂ gas in a gas bag (summer 2018 sampling trip). Bottles were shaken for 2 min to equilibrate headspace gas with water, and headspace gas was extracted and immediately injected into 60 mL glass serum bottles that had been evacuated immediately prior to sample introduction. Samples were stored at room temperature until analysis, which occurred within 1 week of collection. Measured headspace concentrations were converted to dissolved concentrations using methods outlined in Pain et al. (2019b). When atmospheric air was used for headspace extractions, atmosphere samples were collected in tandem and analyzed to correct the calculated dissolved CO₂ and CH₄ concentrations and isotopic compositions for atmospheric CO₂ and CH₄. This correction altered CH₄ concentrations by up to 22 % for one sample from the Russell Glacier, though less than 5 % for all other samples, and resulted in a correction of δ¹³C-CH₄ of up to 1.3 ‰. For CO₂, the correction altered concentrations by up to 15 % for one sample collected at Kiattut Sermiat, though less than 10 % for all other samples, and resulted in a correction of δ¹³C-CO₂ of up to 0.4 ‰.

In samples collected in fall 2017 and summer 2018, alkalinity was measured in the field laboratory within 3 d of collection by titration with 0.01 N HCl using the Gran method. Because alkalinity measurements were not available for the spring 2017 sampling trip, we estimate alkalinity with PHREEQC modeling and the phreeqc.dat database (Parkhurst, 1997) using major cations and anions, pH, temperature, and DIC concentrations as model inputs.

2.3 Laboratory analysis

Gas samples were analyzed for CO₂ and CH₄ concentrations as well as δ¹³C-CO₂ and δ¹³C-CH₄ on a Picarro G2201-i cavity ring-down spectrometer in the field within a few days of collection. Carbon isotopic compositions are reported in reference to Vienna Pee Dee Belemnite (VPDB). Check standards of known CO₂ and CH₄ concentrations and isotopic compositions were measured during each sample run and were accurate within 10 %. Anion and cation concentrations were measured on an automated Dionex ICS-2100 and ICS-1600 ion chromatograph, respectively. Error on replicate analyses was less than 5 %. DIC concentrations were measured on a UIC (coulometrics) 5011 CO₂ coulometer coupled with an AutoMate Preparation Device. Samples were acidified, and the evolved CO₂ was carried through a silver nitrate scrubber to the coulometer where total C was measured. Accuracy was calculated to be ±0.1 mg/L based on measurement of check standards.

2.4 Methane modeling

To assess CH₄ sources and sinks, we calculate ε_c, or the carbon isotopic fractionation factor between CO₂ and CH₄ as defined in Whiticar (1999):

$$\begin{array}{}\text{(5)}& {\mathit{\epsilon }}_{\mathrm{c}}={\mathit{\delta }}^{\mathrm{13}}{C}_{{\mathrm{CO}}_{\mathrm{2}}}-{\mathit{\delta }}^{\mathrm{13}}{C}_{{\mathrm{CH}}_{\mathrm{4}}}.\end{array}$$

Values of ε_c reflect methanogenesis pathways (acetoclastic or CO₂ reduction) as well as the extent of oxidation. Values of ε_c between approximately 40 ‰ and 55 ‰ are produced for CH₄ generated via acetoclastic methanogenesis, while CO₂ reduction produces values between approximately 55 ‰ and 90 ‰. Lower values (ε_c between 5 and 30) result when CH₄ oxidation predominates. Modern atmospheric input without additional alteration of CO₂ or CH₄ isotopic systematics results in a ε_c value of approximately 40 (Whiticar, 1999).

We calculated CH₄ oxidation using the isotopic method outlined in Mahieu et al. (2008) and Preuss et al. (2013). The fraction of oxidized methane (f_ox) in an open system is given by

$$\begin{array}{}\text{(6)}& f_{\mathrm{ox}}={\displaystyle \frac{{\mathit{\delta }}_{\mathrm{E}}-{\mathit{\delta }}_{\mathrm{P}}}{\mathrm{1000}\times ({\mathit{\alpha }}_{\mathrm{ox}}-{\mathit{\alpha }}_{\mathrm{trans}})}},\end{array}$$

where δ_E is the measured δ¹³C-CH₄ value for each water sample, δ_P is δ¹³C-CH₄ of produced methane, α_ox is the oxidation fractionation factor, and α_trans is a fractionation factor resulting from diffusive transportation of CH₄. While the exact value of δ_P is unknown, diagenetic alteration of δ¹³C-CH₄ values through oxidation or transport only enrich δ¹³C-CH₄ signatures; therefore, the value of δ_P is taken as the most depleted δ¹³C-CH₄ signature assuming it is the least impacted by diagenetic alteration. Literature-reported values for α_ox range between 1.003 and 1.049. We calculate the fraction of oxidized methane with the largest fraction factor (α_ox=1.049; Mahieu et al., 2008), which yields the minimum amount of CH₄ oxidation required to explain the observed variations in δ¹³CH₄ and thus is a conservative estimate for CH₄ oxidation, and actual oxidation ratios may be higher. Literature-reported values for α_trans vary from 1 for advection-dominated systems to 1.0178 for diffusion-dominated porous media (de Visscher et al., 2004; Mahieu et al., 2008; Preuss et al., 2013). We assume that transport is advection dominated and thus assume α_trans=1; however, diffusive transport of CH₄ may result in fractionation of CH₄ in the subglacial environment and lead to relatively lower estimates of f_ox. Because hydrologic connectivity between subglacial methanogenic meltwater pockets and drainage features is not well described, the relative importance of advective compared to diffusive CH₄ transport within the subglacial drainage system is not well understood; however, it is presumed to be an advection-dominated system in which expanding drainage networks access and drain methanogenic meltwater pockets throughout the melt season.

2.5 Mineral weathering and carbonate modeling

We used major cation and anion concentrations and alkalinity to partition solutes into the four mineral weathering reactions in Table 1 after correcting solute concentrations for marine aerosol deposition using measured chloride concentrations and standard seawater element ratios. The mass balance model followed the methods of Deuerling et al. (2019). After apportioning solutes to mineral weathering reactions, we used the stoichiometries of reactions to calculate the impact of each reaction on dissolved CO₂ concentrations (Table 1). The mineral weathering model apportions solutes to reactions in Table 1 based on the ratios of $\mathrm{Ca}/\mathrm{Na}$ and $\mathrm{Mg}/\mathrm{Na}$ in silicate minerals in stream bed load samples, which were taken to be 0.54 and 0.38, respectively, for sub-Isunnguata and Russell Glacier samples (Deuerling et al., 2019; Hindshaw et al., 2014; Wimpenny et al., 2010, 2011) and 0.39 and 0.27, respectively, for Kiattut Sermiat samples (Da Prat and Martin, 2019). Because mineral weathering reactions may both add and remove CO₂, we discuss both the net impact of mineral weathering on CO₂ concentrations (Net CO_2-MW), which may have a positive or negative value,

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \left[\text{Net}{\mathrm{CO}}_{\text{2-MW}}\right]=\left[{\mathrm{CO}}_{\text{2-CarbCA}}\right]+\left[{\mathrm{CO}}_{\text{2-CarbSA}}\right]\\ & +\left[{\mathrm{CO}}_{\text{2-SilCA}}\right],\end{array}\end{array}$$

as well as the total impact of mineral weathering on CO₂ concentrations (Total CO_2-MW),

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \left[{\text{Total CO}}_{\text{2-MW}}\right]=\left|\right[{\mathrm{CO}}_{\text{2-CarbCA}}\left]\right|+\left|\right[{\mathrm{CO}}_{\text{2-CarbSA}}\left]\right|\\ & +\left|\right[{\mathrm{CO}}_{\text{2-SilCA}}\left]\right|,\end{array}\end{array}$$

where changes in the concentrations of CO₂ are defined by their absolute values. To discuss the relative importance of individual reactions, we define proportional contributions of each reaction as follows:

$$\begin{array}{}\text{(9a)}& {\displaystyle }\mathit{\%}{\text{CO}}_{\text{2-CarbCA}}={\displaystyle \frac{\mathrm{|}\left[{\mathrm{CO}}_{\text{2-CarbCA}}\right]\mathrm{|}}{\left[{\text{Total CO}}_{\text{2-MW}}\right]}}\times \mathrm{100},\text{(9b)}& {\displaystyle }\mathit{\%}{\text{CO}}_{\text{2-CarbSA}}={\displaystyle \frac{\mathrm{|}\left[{\mathrm{CO}}_{\text{2-CarbSA}}\right]\mathrm{|}}{\left[{\text{Total CO}}_{\text{2-MW}}\right]}}\times \mathrm{100},\text{(9c)}& {\displaystyle }\mathit{\%}{\text{CO}}_{\text{2-SilCA}}={\displaystyle \frac{\mathrm{|}\left[{\mathrm{CO}}_{\text{2-SilCA}}\right]\mathrm{|}}{\left[{\text{Total CO}}_{\text{2-MW}}\right]}}\times \mathrm{100}.\end{array}$$

We combine measured CO₂ concentrations with Net CO_2-MW in order to determine the magnitude of CO₂ production or consumption in the subglacial environment due to processes besides mineral weathering. This analysis assumes that the concentration of CO₂ measured at the subglacial outlet is equivalent to the net change in CO₂ due to mineral weathering plus the sum of all other subglacial CO₂ sources and sinks. We refer to the sum of all other subglacial CO₂ sources and sinks as CO_2-total, which represents the amount of CO₂ that must have been supplied to the subglacial environment to balance the mineral weathering CO₂ sink:

$$\begin{array}{}\text{(10)}& {\text{CO}}_{\text{2-measured}}={\text{Net CO}}_{\text{2-MW}}+{\mathrm{CO}}_{\text{2-total}}.\end{array}$$

The sources of CO₂ to CO_2-total may be evaluated through the use of Keeling plots, which are constructed as the inverse of CO₂ concentrations ([CO₂]⁻¹) versus stable isotopic composition (δ¹³C-CO₂). If variations in the concentration and isotopic composition of CO₂ arise from the mixing of two CO₂ reservoirs with constant isotopic compositions and concentrations (Keeling, 1958), a linear relationship is expected between [CO₂]⁻¹ and δ¹³C-CO₂. The y intercept of a regression between these variables represents the isotopic composition of the high-CO₂ end-member. Because measured CO₂ concentrations include both subglacial CO₂ sources and sinks, which may include considerable consumption through mineral weathering reactions, the magnitude of the total subglacial CO₂ source is taken as CO_2-total. We therefore construct Keeling plots between [CO_2-total]⁻¹ and measured δ¹³C-CO₂ values because while mineral weathering impacts the concentration of CO₂, its isotopic composition is not appreciably altered (Myrttinen et al., 2012) compared to the range of isotopic compositions of potential CO₂ end-members, namely OM remineralization, atmospheric CO₂, and lithogenic CO₂ sources due to mechanical grinding (Fig. 1).

2.6 Discharge relationships

We evaluate the relationship between subglacial CH₄ and CO₂ dynamics and glacial meltwater river discharge records collected downstream of the sub-Isunnguata and upstream of the Russell sampling sites. Proglacial river discharge was collected in the Akuliarusiarsuup Kuua (AK) River at the AK4 site, 2 km downstream of the sub-Isunnguata sampling site (Fig. 2b). The river discharge dataset is an updated and extended version of Rennermalm et al. (2012) using reference and regression models to correct Solinst level logger drift in water stage (Solinst, 2017) and a total of 57 discharge measurements to convert continuous water stage to discharge. The standard uncertainty (i.e., the 68th percent confidence interval or 1 standard deviation) was determined to be 17 % using methods and recommendations provided by Herschy (1999), ISO Guide 98-3 (2008), and WMO (2010).

Because diurnal fluctuations in river discharge can be large, and differing water travel times from subglacial outlet sites to the discharge monitoring site induces a lag between maximum daily discharge at subglacial discharge sites and the AK4 site outlet, we compare subglacial CH₄ and CO₂ concentrations to average daily discharge, calculated as the average of hourly discharge estimates over the days on which subglacial discharge water samples were collected. We use the AK4 site discharge records for evaluating concentration–discharge relationships for both sub-Isunnguata and Russell sites. Although this site is upstream of the Russell Glacier, its close proximity to the Russell Glacier suggests it is more likely to reflect local melting patterns similar to those that would be controlling discharge dynamics at the Russell than discharge records collected at the Watson River outlet. While Watson River discharge records are also available through PROMICE (van As et al., 2018), which includes some contributions from the Russell Glacier, the Watson River includes discharge from the Akuliarusiarsuup Kuua (draining sub-Isunnguata, Russell, and Leverett catchments) as well as the much larger Qinnguata Kuussua catchment, and therefore Watson River discharge records are not likely to be representative of the temporal changes in the magnitude and variability of discharge from the much smaller Russell Glacier catchment.

3 Results

3.1 Temporal variability in water chemistry and gas concentrations

Chemical parameters differ spatially between subglacial discharge sites as well as temporally through the 2017–2018 melt seasons. Comparing the means and standard deviations of water samples collected throughout 2017 and 2018, specific conductivity (Sp.C; Fig. 3a) is typically highest at Kiattut Sermiat (26 ± 8 $\mathrm{\mu }\mathrm{S}/\mathrm{cm}$), followed by Russell (22 ± 5 $\mathrm{\mu }\mathrm{S}/\mathrm{cm}$) and sub-Isunnguata sites (13 ± 9 $\mathrm{\mu }\mathrm{S}/\mathrm{cm}$; Fig. 3a). All sites show variability over time, with lowest values occurring in the summer for sub-Isunnguata and Russell, while Sp.C drops continuously with days of the year for Kiattut Sermiat. Sites differ in pH, and average values at Kiattut Sermiat (8.2 ± 0.4) are higher than both Russell (7.2 ± 0.2) and sub-Isunnguata (6.6 ± 0.6; Fig. 3b), and while values vary over time, no consistent trend is identified between sites. The saturation of dissolved oxygen (D.O.) with respect to atmospheric concentrations is similar between sites, though sub-Isunnguata (98 ± 8 %) values fall below Russell (115 ± 16 %) and Kiattut Sermiat (117 ± 11 %) during all sampling times and exhibit undersaturation in the mid-summer samples, while Russell and Kiattut Sermiat are consistently supersaturated (Fig. 3c). Alkalinity is similar at Russell (93 ± 31 $\mathrm{\mu }\mathrm{eq}/\mathrm{L}$) and Kiattut Sermiat (93 ± 26 $\mathrm{\mu }\mathrm{eq}/\mathrm{L}$), which are higher than at sub-Isunnguata (39 ± 25 $\mathrm{\mu }\mathrm{eq}/\mathrm{L}$), but all reach minimum values in summer (Fig. 3d). CH₄ concentrations differ by orders of magnitude between sites (Fig. 3e) and are consistently supersaturated with respect to atmospheric concentrations at sub-Isunnguata (648 ± 411 ppm or 1575 ± 997 nM) and Russell (58 ± 33 ppm or 110 ± 78 nM) sites, while they are close to atmospheric equilibrium at Kiattut Sermiat (4 ± 2 ppm or 9 ± 5 nM). Mean δ¹³C-CH₄ values (Fig. 3f) are similar between sub-Isunnguata (−54.7 ± 7.5 ‰), Russell (−52 ± 7.3 ‰), and Kiattut Sermiat (−57.6 ± 14.2 ‰). Measured CO₂ concentrations (Fig. 3g) are consistently supersaturated with respect to atmospheric concentrations for sub-Isunnguata (685 ± 230 ppm or 58 ± 18 µM), near atmospheric equilibrium for Russell (442 ± 31 ppm or 29 ± 4 µM), and undersaturated for Kiattut Sermiat (263 ± 33 ppm or 19 ± 2 µM). Mean δ¹³C-CO₂ values (Fig. 3h) are lower in spring and fall for sub-Isunnguata (−16.6 ± 4.0 ‰) compared to Russell (−13.7 ± 2.3 ‰) and Kiattut Sermiat (−16.1 ± 1.6 ‰) sites, though similar seasonal variation occurs for all sites with relatively more depleted values in the spring and fall compared to summer.

Figure 3 Chemical parameters at sub-Isunnguata (sub-IS), Russell (RU), and Kiattut Sermiat (KS) subglacial water sampling sites versus day of year for (a) specific conductivity (Sp.C), (b) pH, (c) dissolved oxygen (D.O.) percent saturation, (d) alkalinity (Alk), (e) measured CH₄ concentrations (left y axis in ppm and right y axis in nM), (f) δ¹³-CH₄ values, (g) measured CO₂ concentrations (left y axis in ppm and right y axis in µM), and (h) δ¹³C-CO₂ values. Atmospheric equilibrium concentrations are indicated by dashed lines and taken as 1.9 ppm for CH₄ and 410 ppm for CO₂. Error bars on CH₄ and CO₂ concentrations and stable isotopic compositions represent the standard deviation of replicates and are smaller than symbols for some data points.

3.2 Methane oxidation and relationship with discharge

Values of ε_c are similar over time for sub-Isunnguata (38 ± 10 ‰) and Russell (38 ± 9 ‰) and are relatively higher in the summer sampling period, while Kiattut Sermiat ε_c values are higher on average (42 ± 13 ‰) with lowest values in the summer (Fig. 4a). Estimates of f_ox are similar between sub-Isunnguata (17 ± 15 %), Russell (23 ± 15 %), and Kiattut Sermiat sites (25 ± 22 %; Fig. 4b). However, f_ox values are higher in the spring and fall sampling times compared to summer for sub-Isunnguata and Russell and approach 50 % in the spring, while Kiattut Sermiat values decrease throughout the melt season.

Figure 4 CH₄ dynamics over the course of the 2017 and 2018 melt seasons including (a) the carbon fractionation factor (ε_c) between dissolved CO₂ and CH₄ and (b) the fraction of CH₄ oxidized (f_ox) for sub-Isunnguata (sub-IS), Russell (RU), and Kiattut Sermiat (KS) samples. Fields of ε_c representing methanogenesis and oxidation values are based on Whiticar (1999). Values of ε_c between approximately 40 and 55 (gray shaded region in panel a) are produced for methanogenesis via acetate fermentation, while CO₂ reduction produces values between approximately 50 and 90. Lower values result from a predominant isotopic signature of CH₄ oxidation. Atmospheric input without additional alteration of CO₂ or CH₄ isotopic systematics results in a ε_c value of approximately 40.

CH₄ concentrations, δ¹³C-CH₄ values, and f_ox are weakly negatively correlated to average daily discharge for both sub-Isunnguata and Russell sites (Fig. 5a, b and d), while ε_c is weakly positively correlated with discharge for both sub-Isunnguata and Russell (Fig. 5c).

Figure 5 Relationships between average daily discharge and CH₄ dynamics including (a) CH₄ concentrations, (b) δ¹³C-CH₄, (c) f_ox, and (d) ε_c for sub-Isunnguata (sub-IS) and Russell (RU) samples. Regressions are shown by dotted lines for Isunnguata and dashed lines for Russell samples. Horizontal error bars represent the standard deviation of average daily discharge for days samples were collected and are smaller than symbols for some data points.

3.3 Mineral weathering impacts on CO₂ and relationship with discharge

Mineral weathering leads to net sequestration of CO₂ at all three sites (Fig. 6a). The magnitude of net ΔCO₂ differs between sites with the lowest average values at sub-Isunnguata (−39 ± 37 µM) followed by Russell (−65 ± 32 µM) and Kiattut Sermiat (−98 ± 17 µM) sites. Individual mineral weathering reactions produce differing contributions between sites and over time, with notable differences between southwest sites (sub-Isunnguata and Russell) and the southern Kiattut Sermiat site (Fig. 6b). For instance, the proportional contribution of Carb_SA is similar between sub-Isunnguata (17 ± 11 %) and Russell (15 ± 6 %) but lower at Kiattut Sermiat (8 ± 1 %; Fig. 6b). Kiattut Sermiat has a relatively greater contribution from Carb_CA (62 ± 2 %) compared to sub-Isunnguata (41 ± 10 %) and Russell (38 ± 6 %), while Sil_CA is lower at Kiattut Sermiat (28 ± 1 %) compared to sub-Isunnguata (41 ± 17 %) and Russell (47 ± 11 %). Kiattut Sermiat additionally exhibits low seasonal variability in the proportional contributions of individual mineral weathering reactions compared to sub-Isunnguata and Russell sites.

Figure 6 Mineral weathering model results in (a) net impact of mineral weathering reactions on CO₂ (Net CO_2-MW; Eq. 7) for Isunnguata, Russell, and Kiattut Sermiat subglacial discharge sites (where negative values of Net CO_2-MW indicate net sequestration of CO₂ due to mineral weathering) and (b) the proportional contribution of each mineral weathering reaction to the total change in CO₂ (% Total CO_2-MW Eqs. 9a–9c).

CO_2-total represents CO₂ concentrations in the subglacial environment prior to addition and/or consumption of CO₂ through mineral weathering (Eq. 10; Fig. 7). Because the Net CO_2-MW is always negative (more consumption than production), the value of CO_2-total is always greater than measured concentrations (CO_2-measured). Regardless of differences in CO_2-measured between sites, the average CO_2-total values are similar between sites and average 91 ± 47 µM for sub-Isunnguata, 94 ± 33 µM for Russell, and 117 ± 16 µM for Kiattut Sermiat.

Figure 7 Calculated CO_2-total values for (a) sub-Isunnguata, (b) Russell, and (c) Kiattut Sermiat subglacial discharge against day of the year.

For both sub-Isunnguata and Russell sites, average daily discharge is negatively correlated with CO₂ concentrations (Fig. 8a), while it is positively correlated with δ¹³C-CO₂ (Fig. 8b).

Keeling plots between [CO_2-total]⁻¹ and δ¹³C-CO₂ for each site indicate no linear relationship for Russell or Kiattut Sermiat samples; however, a strong linear correlation is observed for sub-Isunnguata (r²=0.99; p<0.001) samples with the removal of one outlier, which also had the lowest CO_2-total value (Fig. 9).

Figure 8 Relationships between average daily discharge (a) CO₂ concentrations and (b) δ¹³C-CO₂. Regressions are shown by dotted lines for sub-Isunnguata and dashed lines for Russell samples. Horizontal error bars represent the standard deviation of average daily discharge for days samples were collected and are smaller than symbols for some data points.

4 Discussion

We observe orders of magnitude variability in dissolved CH₄ and CO₂ concentrations in subglacial discharge of the Greenland Ice Sheet, indicating significant differences in the magnitudes of the sources and sinks of these gases across time and space. Supersaturation of both CO₂ and CH₄ with respect to atmospheric concentrations indicates that sub-Isunnguata discharge is a source of both gases to the atmosphere, the neighboring Russell Glacier discharges water that is a source of CH₄ but near equilibrium with respect to CO₂, and Kiattut Sermiat in southern Greenland is a sink of atmospheric CO₂ but near equilibrium with respect to CH₄ (Fig. 3e, g). Because CH₄ dynamics may be largely microbially driven, while CO₂ dynamics include microbial as well as abiotic mineral weathering processes, we first discuss CH₄ dynamics including a comparison of concentrations, isotopic compositions, and extent of oxidation between sites and over the melt season. We then discuss CO₂ concentrations, impacts of mineral weathering reactions (Table 1), and an assessment of subglacial CO₂ sources, including OM remineralization. These assessments will contribute to our understanding of the variability and controls of CH₄ and CO₂ concentrations in subglacial discharge from the Greenland Ice Sheet and may improve predictions of the impact of future ice melt on Arctic carbon budgets.

4.1 Sources and sinks of CH₄

Differences in CH₄ concentrations and relationships with discharge between sites imply heterogeneity in both the extent and controls of subglacial methanogenesis under the Greenland Ice Sheet. CH₄ supersaturation occurs at the two subglacial discharge sites that flow to the Akuliarusiarsuup Kuua (sub-Isunnguata and Russell), and concentrations are similar to the ranges reported in discharge of the Leverett Glacier (up to 600 nM; Lamarche-Gagnon et al., 2019), located near the Russell Glacier in this study (Fig. 2b). However, CH₄ concentrations are near atmospheric equilibrium for the Kiattut Sermiat site (Fig. 3e). Because methanogenesis is an anaerobic OM remineralization pathway, it is more likely to occur in subglacial environments isolated from atmospheric O₂ sources. Widespread observations of methanogenesis in glacial meltwater of southwest Greenland from this and other studies (Christiansen and Jørgensen, 2018; Dieser et al., 2014; Lamarche-Gagnon et al., 2019), and limited observations of CH₄ in subglacial discharge in southern Greenland, suggests heterogeneity in subglacial conditions that support methanogenesis. Methanogenesis fueled by organic material overridden during ice sheet growth has been suggested as a potential climate feedback over glacial interglacial timescales (Wadham et al., 2008) and may contribute to observed variations in CH₄ concentrations.

Subglacial methane concentrations may additionally be controlled by hydrologic factors as the subglacial hydrological network develops throughout the melt season and channelization of meltwater conduits increases subglacial drainage efficiency (Andrews et al., 2015; Cowton et al., 2013). Drainage efficiency impacts both subglacial water residence time as well the transport of aerobic supraglacial meltwater to the ice bed. Both residence time and oxygen delivery may impact subglacial redox status and methanogenesis potential and favor methanogenesis when oxygen supply rates are low compared to OM remineralization rates. This condition is most likely to be met in distributed subglacial systems that are hydrologically isolated with limited inputs from aerobic supraglacial meltwater. Such a hydrologic control on methanogenesis is supported by the negative correlation between CH₄ concentrations and average daily discharge at both sites (Fig. 5a). This correlation would suggest that either CH₄ production occurs predominantly during periods of low discharge and greater residence time or higher discharge results in the dilution of a CH₄-laden subglacial water source. While both mechanisms would lead to a similar relationship between discharge and CH₄ concentrations, they carry different implications for subglacial methanogenesis. If limited by residence time, a hydrologic link between glacial hydrology and subglacial biogeochemistry would be established because supraglacial discharges deliver terminal electron acceptors to the ice bed and would limit methanogenesis. If predominantly controlled by dilution, however, active methanogenesis would not be required – only the existence of meltwater pockets containing CH₄ that may or may not have been recently produced. Further discussion of these mechanisms is outside the scope of this study. While we have limited data to make further inferences about hydrologic controls of methanogenesis, the presence of several outliers at the sub-Isunnguata site in particular (Fig. 5a) highlights the possibility for additional controls including stochastic drainage events or heterogeneity in subglacial CH₄ concentrations that result in variability in the relationship between concentration and discharge, as was observed in Lamarche-Gagnon et al. (2019).

Our results suggest heterogeneity in the extent of methanogenesis between outlet glaciers but homogeneity of the microbial methanogenesis pathway as well as CH₄ oxidation dynamics between sites. Methanogenesis pathways may be evaluated by δ¹³C-CH₄ as well as ε_c values because they impart distinct δ¹³C signatures to CH₄ and CO₂ (Whiticar and Schoell, 1986). Dieser et al. (2014) measured a microbial δ¹³C-CH₄ production signal at the Russell Glacier with values between −63 ‰ and −64 ‰, which was interpreted to reflect a possible combination of CH₄ produced through both acetoclastic and CO₂ reduction pathways. The most depleted δ¹³C-CH₄ value measured at the sub-Isunnguata in this study was −62.7 ‰, close to values measured by Dieser et al. (2014) (Fig. 3f) and similar to values reported by Lamarche-Gagnon et al. (2019) for the Leverett Glacier. The similar isotopic ratio between our samples and that measured in active methanogenic communities could indicate that similar methanogenesis pathways occur across this region or that the δ¹³C-CH₄ of stored subglacial CH₄ has not been fractionated by oxidation or transport in the peak melt season when we observe these depleted δ¹³C-CH₄ values.

While the exact contributions from each methanogenesis pathway cannot be inferred from isotopic information alone, the range of ε_c values at outlet glaciers are consistent with predominantly acetoclastic methanogenesis during the peak melt season (Fig. 4a). However, ε_c values fall below the expected range from acetoclastic methanogenesis during the early and late melt seasons and may result from variations in the extent of subglacial CH₄ oxidation. Seasonal variation in CH₄ oxidation is supported by consistency between ε_c and f_ox values, which both indicate the greatest impact of oxidation (approaching 50 %) in the early melt season compared to peak melt season (Fig. 4a, b), with additional evidence of elevated CH₄ oxidation in the late melt season at both sub-Isunnguata and Russell sites. Because our water sampling locations were slightly downstream of glacial discharge outlets, there is also the possibility that outgassing in between the outlet and our sampling location reduced dissolved CH₄ concentrations and led to more enriched isotopic compositions of remaining dissolved CH₄. It is likely that some outgassing did occur; however, it is unlikely that the extent of outgassing between the glacial outlet and our sampling location would vary significantly between sampling times, and thus outgassing would not fully explain temporal differences in concentration, δ¹³C-CH₄, f_ox, or ε_c. While our measured gas concentrations and isotopic compositions likely reflect some modification of CH₄ and CO₂ isotopic compositions due to outgassing, the differences over time are more likely due to changes in subglacial CH₄ dynamics than outgassing.

The extent of CH₄ oxidation may be controlled by multiple factors including oxygen availability, subglacial residences time, and subglacial hydrology, similar to methanogenesis. A hydrologic control of CH₄ oxidation is supported by relationships between f_ox and ε_c with average daily discharge (collected at site AK4; Fig. 2b) at both sub-Isunnguata and Russell sites: f_ox is negatively correlated with discharge for both sites (Fig. 5b) while ε_c is positively correlated with discharge (Fig. 5c). While weak, these correlations suggest that CH₄ oxidation is greatest during periods of low flow, which may be associated with longer residence times to allow subglacial CH₄ oxidation; however, this relationship could also result from differences in CH₄ sources throughout the melt season as the subglacial drainage network expands. Assuming the former, the delivery of oxygen to the subsurface by supraglacial melting does not appear to be a limiting factor in subglacial CH₄ oxidation, which should increase f_ox as more oxygenated supraglacial water is delivered to the subglacial system. Instead, the observed greater CH₄ oxidation during periods of low discharge may reflect mixing of methane-rich subglacial meltwater pockets and aerobic subglacial meltwater leading to CH₄ oxidation. Longer transit times during periods of low flow may allow more subglacial methane oxidation to occur than during peak discharge, when the development of channelized flow paths reduces meltwater residence time in the subglacial environment.

Our results indicate a high degree of heterogeneity in subglacial methanogenesis under the Greenland Ice Sheet, as well as a significant impact of CH₄ oxidation, which serves to reduce atmospheric CH₄ fluxes. Given the observed spatial and temporal heterogeneity of CH₄ concentrations and processes, further investigation of the spatial variability in outlet glacier CH₄ concentrations is needed to determine the impact of Greenland Ice Sheet loss on Arctic and global CH₄ budgets, while a better understanding of the controls of these differences will improve models of how CH₄ fluxes from subglacial discharge will change with continued warming.

4.2 Sources and sinks of CO₂

Dissolved CO₂ concentrations in subglacial discharge are consistently supersaturated with respect to atmospheric concentrations at the sub-Isunnguata site, near atmospheric equilibrium at the Russell Glacier, and undersaturated at Kiattut Sermiat, indicating that glacial meltwater from the Greenland Ice Sheet can serve as either a source or sink of CO₂ to the atmosphere. Similar to CH₄, differences in dissolved CO₂ dynamics (Fig. 3g) imply variability in carbon processes under the Greenland Ice Sheet. We first discuss potential subglacial CO₂ sources, including OM remineralization, followed by a discussion of CO₂ consumption due to mineral weathering.

4.2.1 Subglacial CO₂ sources

There are many potential sources of CO₂ in the subglacial environment including dissolution of atmospheric gases in air-filled conduits or fractures in ice, CO₂ contained in ice bubbles (Fig. 1; Anklin et al., 1995; Graly et al., 2017), mechanical grinding and volatilization of fluid inclusions in bedrock (Macdonald et al., 2018), and OM remineralization. While previous studies have indicated that additional atmospheric CO₂ input through fractures and air-filled conduits may supply sufficient CO₂ to drive mineral weathering observed in many subglacial environments, including several sites in Greenland (Graly et al., 2017), CO₂ is also a product of OM remineralization, which is believed to account for CH₄ concentrations elevated above atmospheric equilibrium at the two southwest sites in this study. Both CO₂ and CH₄ exhibit negative correlations with average daily discharge for both sub-Isunnguata and Russell sites and could suggest a common OM remineralization source (Fig. 8a). While the magnitude of this source and its relative importance compared to other subglacial CO₂ sources is currently unknown, differing sources of carbonic acid for mineral weathering reactions carry different implications for subglacial CO₂ budgets. For instance, carbonic acid weathering driven by invasion of atmospheric CO₂ would represent a sink of atmospheric CO₂, but carbonic acid weathering driven by OM remineralization would instead serve to consume CO₂ from in situ sources and limit the potential for subglacial meltwater to be an atmospheric CO₂ source once discharged from the glacier. Determining the sources of carbonic acid to subglacial weathering reactions is therefore critical to understand the controls of mineral weathering in subglacial environments and its role in atmospheric CO₂ sequestration.

Comparisons between measured δ¹³C-CO₂ in subglacial discharge samples and likely δ¹³C-CO₂ values of CO₂ sources indicate that CO₂ sources differ between sites, with OM remineralization as the most important CO₂ source at the sub-Isunnguata but likely not the predominant or sole source at Russell or Kiattut Sermiat glaciers. Keeling plots of [CO_2-total]⁻¹ versus δ¹³C-CO₂ indicate that CO_2-total may be represented by a two-end-member mixing model for sub-Isunnguata discharge but not for discharge from the Russell and Kiattut Sermiat glaciers (Fig. 9). Mixing model end-members include a ¹³C-enriched, lower-concentration CO₂ source and a ¹³C-depleted, higher-concentration CO₂ source. The y intercept of the regression between [CO_2-total]⁻¹ versus δ¹³C-CO₂, which represents the isotopic signature of the high-CO₂ end-member, is −27.4 ‰. This value is close to what would be expected from OM remineralization as indicated by remineralized OM in Greenlandic heath soils that produced δ¹³C-CO₂ of approximately −27‰ to −25 ‰ (Ravn et al., 2020) and thawed Alaskan permafrost soils that produced δ¹³C-CO₂ of between −20 ‰ to −30 ‰ (Mauritz et al., 2019), both of which may be similar to subglacial organic matter. The low-CO₂ end-member could reflect atmospheric CO₂ input, with an associated δ¹³C-CO₂ value of approximately −8 ‰. While the δ¹³C-CO₂ value of the lowest-CO_2-total samples in the sub-Isunnguata Keeling plot (e.g., highest [CO_2-total]⁻¹ not including the outlier) is slightly depleted compared to atmospheric values at −12.1 ‰, even the lowest CO₂ concentrations measured at sub-Isunnguata are supersaturated with respect to atmospheric concentrations (Fig. 3g). Consistent CO₂ supersaturation suggests that OM remineralization contributes CO₂ even for low-CO₂-concentration samples and isotopically depletes the subglacial CO₂ reservoir.

Figure 9 Keeling plot indicating correlations between the magnitude of CO_2-total and δ¹³C-CO₂ and for sub-Isunnguata (sub-IS), Russell (RU), and Kiattut Sermiat (KS) samples. Asterisk denotes the outlier not included in the regression between CO_2-total and δ¹³C-CO₂ for sub-Isunnguata samples. The plotted regression line was constructed using the sub-IS samples only.

While δ¹³C-CO₂ values of Russell and Kiattut Sermiat samples are within the range of sub-Isunnguata samples, suggesting possible contributions of CO_2-atm and CO_2-OM, scatter in the Keeling plots indicates variability in the CO₂ concentration and/or isotopic composition of end-members, or significant contributions of at least one other major subglacial CO₂ source. We address both possibilities here. While atmospheric CO₂ concentrations and δ¹³C values should be relatively invariable, CO_2-OM may vary both in concentration and isotopic composition, depending on variability in the quantity and composition of organic deposits as well as remineralization rates. For instance, if remineralization largely occurs in hydrologically isolated subglacial meltwater pockets, some variability in the concentration and δ¹³C-CO₂ of CO_2-OM is likely. While no data yet exist to characterize the variability in subglacial OM reservoirs, variability in either concentration or isotopic composition of CO_2-OM could plausibly result in the scatter shown in Fig. 9. Additional subglacial CO₂ sources from ice bubbles, or lithogenic CO₂ liberated by mechanical grinding, would be expected to enrich rather than deplete the δ¹³C-CO₂ values of the samples relative to modern atmospheric δ¹³C-CO₂ values. Ice bubbles contain gaseous CO₂ at concentrations and isotopic compositions reflecting atmospheric conditions during ice formation. While heterogeneity may result from gas bubbles recording changes in atmospheric CO₂, variability in δ¹³C-CO₂ of gas bubble CO₂ should be only a few per mill, which is small compared to the variation observed in Russell and Kiattut Sermiat samples (e.g., Tipple et al., 2010; Fig. 9). Gas bubble CO₂ should also be ¹³C-enriched compared to modern atmospheric CO₂ due to fossil fuel contributions and thus would be unlikely to cause the observed depletion of ¹³C in the subglacial water.

Recent work has also highlighted the potential for subglacial mechanical grinding to result in CO₂ production through the volatilization of bedrock fluid inclusions (Macdonald et al., 2018). While this process was found to produce sufficient CO₂ to drive approximately 20 % of mineral weathering in Svalbard subglacial environments, the expected isotopic composition of lithogenic CO₂ is more ¹³C-enriched than our measured δ¹³C-CO₂ values. For instance, estimates of δ¹³C for bulk hydrocarbons in fluid inclusions in the Ilímaussaq alkaline complex of South Greenland have values of −4.5 ± 1.5 ‰ (Madsen, 2001), which is close to the δ¹³C-CO₂ of CO₂ in fluid inclusions in the Bamble granulite sector of South Norway ($\sim -$6 ‰; Newton et al., 1980). Because mechanical grinding should not fractionate the δ¹³C-CO₂ values (Lüders et al., 2012), our low δ¹³C-CO₂ values suggest this source is small relative to other sources.

One additional source or sink of CO₂ to some of our samples is atmospheric exchange as water flows from the subglacial outlet site to our sampling sites. However, atmospheric CO₂ exchange after discharge would have the same impact on Keeling plots as atmospheric CO₂ exchange prior to discharge. Incorporation of an atmospheric source between the outlet and sample site would be most likely at Kiattut Sermiat, where CO₂ concentrations are undersaturated with respect to atmospheric concentrations, which would promote invasion of atmospheric CO₂. However, the measured δ¹³C-CO₂ values are more ¹³C-depleted than modern atmospheric CO₂ and are not consistent with atmospheric CO₂ as the sole or dominant source of CO₂ to these samples (Fig. 9).

While more information is needed to determine all possible sources of CO₂ to Russell and Kiattut Sermiat samples, δ¹³C-CO₂ values of samples from both sites imply mixing between a ¹³C-depleted CO₂ source, such as OM remineralization, and one or more ¹³C-enriched CO₂ sources, such as atmospheric or lithogenic CO₂. Similar to CH₄, concentrations and isotopic compositions of gases may be impacted by atmospheric exchange between the glacial outlet and our sampling locations, which would alter dissolved CO₂ concentrations and δ¹³C-CO₂ compositions to values more similar to atmospheric values. Therefore, we are unable to distinguish the impacts of the atmospheric exchange that occurs prior to discharge from the exchange that occurs between discharge and our sampling locations. However, the impact of this exchange should be relatively constant between sampling times and sampling locations; therefore, outgassing would not account for temporal or spatial variability in CO₂ concentrations or isotopic compositions between sites.

4.2.2 Subglacial CO₂ sink: mineral weathering reactions

Mineral weathering leads to net CO₂ consumption in all subglacial discharge samples (Fig. 6), and thus the measured CO₂ concentrations at glacial outlets represent only a fraction of the total CO₂ that would have been present in the absence of mineral weathering reactions (CO_2-total; Eq. 10). Net consumption occurs because the CO₂ source from Carb_SA is ubiquitously low compared to sinks from either Carb_CA or Sil_CA (Fig. 6b). The range in Net CO_2-MW is similar between subglacial discharge sites (between 10–150 µM; Fig. 6a), but average values increase from Kiattut Sermiat to Russell to sub-Isunnguata, likely reflecting the relative weatherability of alkaline igneous rocks, granulite facies gneisses, and amphibolite facies gneisses. Kiattut Sermiat is characterized by a relatively high proportion of Carb_CA compared to sub-Isunnguata and Russell sites, which may arise from the presence of trace carbonates in abundant and readily weatherable basaltic intrusions as has been implicated in other studies (Urra et al., 2019). The relatively greater influence of carbonate dissolution compared to silicate dissolution on Total CO_2-MW at Kiattut Sermiat may also relate to more rapid dissolution kinetics of carbonates, which allow carbonate dissolution to have a large influence on major cation and anion load even when carbonates are only present in trace amounts (Deuerling et al., 2019; Tranter, 2005). At sub-Isunnguata and Russell sites, Sil_CA has a greater influence than Carb_CA on Total CO_2-MW, which could result from either a lower abundance of trace carbonates to participate in weathering reactions or relatively longer subglacial residence times that would allow a greater accumulation of silicate weathering products.

Despite the impact of Carb_CA on Total CO_2-MW at Kiattut Sermiat compared to sub-Isunnguata and Russell sites, Carb_SA is notably lower at Kiattut Sermiat than other sites and suggests a limited role for sulfuric acid weathering that may relate to subglacial sulfide oxidation dynamics. Lower abundances of sulfide minerals in the subglacial environment may limit the production of sulfuric acid and could result from differences in lithology between sites, the depletion of sulfide minerals due to prior weathering (Graly et al., 2014), or weathering occurring in anoxic environments that limit the oxidation of sulfide to sulfuric acid (Deuerling et al., 2019). The kinetics of sulfide oxidation may also significantly differ between sites depending on the relative contributions of abiotic compared to microbially mediated sulfide oxidation, as microbially mediated sulfide oxidation is several orders of magnitude faster than abiotic sulfide oxidation (Boyd et al., 2014; Harrold et al., 2016). Rapid microbially mediated sulfide oxidation has been implicated in the development of anaerobic conditions, which could also support subglacial methanogenesis (Sharp et al., 1999; Stibal et al., 2012; Wadham et al., 2010). Observations of higher CH₄ concentrations as well as higher contributions of Carb_SA at sub-Isunnguata and Russell compared to Kiattut Sermiat may therefore be linked to subglacial microbial activity, which is known to vary based on factors such as the presence of organic and fine-grained rock flour to serve as growth substrates, insulation from fluctuations in temperature, and delivery of nutrients and organic matter from supraglacial sources (Sharp et al., 1999). If microbially driven, our results suggest possible linkages between microbial processes and subglacial mineral weathering regimes, with significant impacts to both CH₄ and CO₂ dynamics due to the role of Carb_SA as a CO₂ source (Table 1).

Whether controlled by geochemical, microbial, or mechanical processes, the relationships between CO₂ concentrations and δ¹³C-CO₂ and average daily discharge are similar between sub-Isunnguata and Russell sites (Fig. 8). These similarities suggest that underlying controls in carbonate chemistry may be consistent between sites despite the heterogeneity in measured dissolved CO₂ concentrations. For both sub-Isunnguata and Russell sites, mineral weathering reactions consume CO₂, which implies contributions from in situ CO₂ sources (such as atmospheric CO₂ invasion or OM remineralization) to produce measured CO₂ concentrations. The different CO₂ concentrations observed between sites therefore appear to result from the strength of in situ CO₂ sources relative to CO_2-MW, both of which impart the greatest chemical change at times of low discharge and high subglacial residence time. At Kiattut Sermiat, where measured CO₂ concentrations are lowest, the magnitude of in situ sources is insufficient to maintain atmospheric equilibrium, leading subglacial discharge to be a sink of atmospheric CO₂, while CO_2-total maintains close to atmospheric equilibrium concentrations at the Russell Glacier. At sub-Isunnguata, however, OM remineralization produces more CO₂ than is consumed by mineral weathering and causes meltwater to be a source of CO₂ to the atmosphere. This finding implies that subglacial mineral weathering serves to partially or fully consume CO₂ produced from in situ sources under the Greenland Ice Sheet but does not always serve to directly consume modern atmospheric CO₂.

5 Conclusions

Subglacial reactions impact the concentrations of CO₂ and CH₄ in subglacial discharge of the Greenland Ice Sheet, and differences in the relative magnitudes of microbial and geochemical processes result in a high degree of previously unrecognized heterogeneity between glacial discharge sites of the Greenland Ice Sheet. Our results imply a significant role of OM remineralization in driving this heterogeneity and leading to CO₂ and CH₄ supersaturation at the sub-Isunnguata site and CH₄ supersaturation at the Russell site. Heterogeneity may result in significant uncertainty in total greenhouse gas flux estimates from subglacial systems of the Greenland Ice Sheet, which will be an increasingly important carbon flux as the Arctic warms in the coming decades. While heterogeneous, the uncertainty in greenhouse gas fluxes from Greenland Ice Sheet meltwater may be reduced by a better understanding of the controls and variability of the weathering reactions and microbial processes driving heterogeneous gas concentrations. Such a process-based understanding could also improve estimates of the impact of greenhouse gas variability associated with the growth and retreat of continental ice sheets over glacial–interglacial timescales. Subglacial OM remineralization further implies the existence of links between subglacial OM deposits and export of other biogeochemical solutes from the Greenland Ice Sheet, including nutrients as well as redox-sensitive elements. While the export of nutrients from the Greenland Ice Sheet has been the focus of numerous studies (Bhatia et al., 2013; Hawkings et al., 2016; Lawson et al., 2014), little is currently known regarding the role of OM sources in governing these exports. Given the variability in greenhouse gas concentrations observed in this study, constraining the extent of heterogeneity in outlet glaciers of the Greenland Ice Sheet as well as the biogeochemical, hydrologic, and geologic controls of this heterogeneity will be important for upscaling atmospheric fluxes as well as efforts to predict impacts of ice loss on carbon and nutrient budgets due to current and future melting of the Greenland Ice Sheet.