1. Introduction
The Greenland ice sheet (GrIS) contributes significant quantities of meltwater to the surrounding North Atlantic and Arctic Oceans (Reference SolomonSolomon and others, 2007). During the high-melt year of 2007, for example, the GrIS contributed an estimated 523 km3 a−1 surface runoff (Reference Mernild, Liston, Hiemstra, Steffen, Hanna and ChristensenMernild and others, 2009), equivalent to the combined mean annual discharge from the four large North American pan-Arctic rivers, Yukon, Mackenzie, Peel and Beck (Shiklomanov, 2009). Recent evidence that a large fraction of annual surface meltwater likely drains to the bed of the GrIS (Reference McMillan, Nienow, Shepherd, Benham and SoleMcMillan and others, 2007; Reference DasDas and others, 2008; Reference Krawczynski, Behn, Das and JoughinKrawczynski and others, 2009) suggests that significant portions of the GrIS subglacial hydrological system may undergo a seasonal evolution, akin to those observed beneath alpine glaciers (Reference Shepherd, Hubbard, Nienow, McMillan and JoughinShepherd and others, 2009; Reference Bartholomew, Nienow, Mair, Hubbard, King and SoleBartholomew and others, 2010). This routing suggests there is potential for substantial meltwater interaction with underlying subglacial till and bedrock along seasonally evolving flow paths. Despite recent interest in the GrIS subglacial hydrological system (Reference Box and SkiBox and Ski, 2007; Reference Joughin, Das, King, Smith, Howat and MoonJoughin and others, 2008; Reference Shepherd, Hubbard, Nienow, McMillan and JoughinShepherd and others, 2009; Reference Bartholomew, Nienow, Mair, Hubbard, King and SoleBartholomew and others, 2010; Reference Tsai and RiceTsai and Rice, 2010), understanding of the subglacial drainage seasonal evolution remains limited and poorly constrained.
In alpine glacial systems, subglacial flow paths can vary seasonally between two end-member systems: channelized drainage (quick flow) and distributed drainage (delayed flow) (Reference PatersonPaterson, 1994). Channelized drainage systems are a series of large tunnels incised into the bedrock, the overlying ice or the till, which facilitate localized and rapid water flow (∼1 m s−1) to the glacier front. They are prevalent at the peak of the summer melt season (Reference PatersonPaterson, 1994; Reference Benn and EvansBenn and Evans, 1998) and transport the bulk of surface ice melt that has drained to the bed (Reference RichardsRichards and others, 1996; Reference Nienow, Sharp and WillisNienow and others, 1998). In comparison, distributed drainage systems, envisioned as a series of linked cavities, lie along the ice/bed interface, and may include a constant source of water from basal ice melt and groundwater in contact with glacial till (Reference PatersonPaterson, 1994). Such systems have characteristically slower water transit times (∼0.01 m s−1), higher water pressures and are water-full for most of the year (Reference RichardsRichards and others, 1996; Reference Benn and EvansBenn and Evans, 1998). Additionally, they may transport a significant proportion of early-season snowmelt to the glacier front (Reference Nienow, Sharp and WillisNienow and others, 1998).
The seasonal evolution of subglacial drainage conditions has important geophysical and biogeochemical implications. From a geophysical perspective, the partitioning of meltwater between these two different drainage systems strongly influences basal water pressures and thus sliding velocities (Reference PatersonPaterson, 1994). From a biogeochemical perspective, the degree of water/rock contact dictates the chemical enrichment of discharge waters exported to surrounding marine ecosystems, and may fuel subglacial microbial processes (Reference Sharp, Parkes, Cragg, Fairchild, Lamb and TranterSharp and others, 1999; Reference Skidmore, Foght and SharpSkidmore and others, 2000). Microbial communities could in turn facilitate the release of additional nutrients and metals, amplify chemical weathering reactions and/or utilize previously overridden organic carbon (Reference Tranter, Sharp, Lamb, Brown, Hubbard and WillisTranter and others, 2002; Reference Wadham, Tranter, Tulaczyk and SharpWadham and others, 2008; Reference Bhatia, Das, Longnecker, Charette and KujawinskiBhatia and others, 2010). Thus, the geochemical cycles of major and minor elements in the coastal waters surrounding the GrIS may be strongly influenced by the temporal dynamics of subglacial discharge (as observed in other regions: Reference RaiswellRaiswell and others, 2006; Reference Hood and ScottHood and Scott, 2008; Reference HoodHood and others, 2009), and in particular by the release of water that has been stored at the bed.
The development of a chemical mixing model that can successfully differentiate water source contributions and subglacial flow paths will complement existing geophysical and active-tracer methods used to study seasonally evolving subglacial hydrological systems (Reference Nienow, Sharp and WillisNienow and others, 1998; Reference Bartholomew, Nienow, Mair, Hubbard, King and SoleBartholomew and others, 2010). The interaction of surface meltwater with the glacier bed alters its chemical composition from dilute snow- and ice melt to chemically enriched subglacial discharge waters. Thus, in theory, variations in solute concentrations could be used to infer the evolution of the subglacial drainage network by differentiating water source contributions. This approach is promising because distributed drainage systems produce discharge waters with significantly enriched chemical signatures due to the longer residence time at the bed compared to channelized drainage systems. Initial mixing-model efforts based the separation of discharge components on bulk properties such as electrical conductivity (EC) (Reference CollinsCollins, 1979; Reference Gurnell and FennGurnell and Fenn, 1984). However, a model with EC as its defining chemical parameter is poorly constrained because EC is not conservative in glacial systems, and is subject to temporal variations in chemical signature and to post-mixing solute acquisition (Reference Sharp, Brown, Tranter, Willis and HubbardSharp and others, 1995). Consequently, individual dissolved species (e.g. sulfate, chloride) whose provenance is exclusive to specific discharge components have been increasingly used instead of EC (Reference Tranter and RaiswellTranter and Raiswell, 1991; Reference Brown, Sharp, Tranter, Gurnell and NienowBrown and others, 1994; Reference TranterTranter and others, 1997; Reference Mitchell, Brown and FugeMitchell and others, 2001). Despite this progress, discharge component separations remain challenging because individual solutes may not retain unique signatures over time and may be involved in subglacial biological reactions across glacial catchments (Reference Sharp, Brown, Tranter, Willis and HubbardSharp and others, 1995; Reference BrownBrown, 2002). Thus, the potential of hydrochemical separation methods in glacier systems has yet to be fully realized (Reference Sharp, Brown, Tranter, Willis and HubbardSharp and others, 1995; Reference BrownBrown, 2002).
Here we introduce a new multi-component isotope-mixing model combining the radioactive isotope radon-222 (222Rn) with the stable water isotopes oxygen-18 (18O) and deuterium (D) as passive flow tracers. These natural isotopic tracers have unique end-member signatures for different reservoirs ultimately contributing to the total glacial outflow, and are not subject to additional post-mixing enrichment or alteration. We then apply this model to quantify the relative contributions of different water reservoirs to the bulk meltwater discharge at a small land-terminating outlet glacier along the western margin of the GrIS in three stages: (1) identified conceptual end-member water reservoirs (surface snow, glacial ice and delayed flow), (2) identified unique passive flow isotopic tracers for each water reservoir (δ18O, D and 222Rn), and (3) applied end-member mixing analysis to estimate relative contributions of each conceptual water reservoir. We limit our separation to these three conceptual reservoirs because hydrographic separation of all potential drainage components is not possible without additional isotopes or end-member analyses. Finally, we investigate the potential for using beryllium-7 (7Be), a naturally occurring radioisotope produced in the atmosphere, as a tracer for the transit of snowmelt through the subglacial drainage system. To our knowledge, this is the first time that beryllium has been used in a glacial system as a tracer of hydrological flow. Through the application of radioisotopes, this pilot study represents a new direction in the use of chemical mixing models to delineate subglacial drainage structure, even though full hydrograph separation of all potential drainage components remains elusive.
2. Model Theory
The basic tenets of radon radiochemistry suggest that it has the potential to be an effective tracer of delayed-flow basal waters characteristic of distributed subglacial drainage systems. Radon, an inert noble gas, is a daughter product of radium-226 within the uranium-238 decay series that is naturally present in soil, sediment and rocks. Since the uranium content of solids and the degree of water/rock interaction will determine the amount of radon enrichment in a given water parcel, groundwater will be highly enriched in radon, as should any surface meltwater or basal ice melt that has been stored subglacially. In contrast, surface snowmelt and glacial ice melt that is quickly routed through the subglacial environment has minimal lithogenic sediment contact, and so should be relatively devoid of radon. Radon has been successfully applied as a tracer of submarine groundwater discharge (Reference Cable, Burnett, Chanton and WeatherlyCable and others, 1996; Reference Corbett, Burnett, Cable and ClarkCorbett and others, 1997; Reference Burnett and DulaiovaBurnett and Dulaiova, 2003) but has only recently been applied in a glacial setting (Reference Kies, Hengesch, Tosheva, Jania, Nawrot, Barnet and PacherováKies and others, 2010). Moreover, the high uranium content of the Greenland fractured silicate bedrock suggests that waters in contact with the GrIS subglacial environment will be particularly enriched in radon (Reference Kraemer, Genereux, Kendall and McDonnellKraemer and Genereux, 1998).
Our new model is also based on the hypothesis that the surface snow and glacial ice at a GrIS outlet glacier have unique δ18O and δD signatures. Stable water isotopes have been used extensively in glacial settings for a variety of applications from ice-core paleotemperature reconstructions (e.g. Reference JouzelJouzel and others, 1997) to delineation of drainage basins on the GrIS (Reference Reeh and ThomsenReeh and Thomsen, 1993). The GrIS margin, in contrast to alpine catchments (Reference TheakstoneTheakstone, 2003), possesses large isotopic differences between snow and ice components (e.g. Reference Reeh and ThomsenReeh and Thomsen, 1993). This is because seasonal surface snow carries the enriched signature of present-day precipitation at relatively high temperatures and low elevation across the ablation zone, whereas underlying marginal glacial ice will have comparatively more depleted values, reflecting its depositional and flow history from higher elevations and/or colder times (Dansgaard and others, 1971).
3. Field Site and Sampling Overview
Sampling for this study was conducted at two primary field locations along the GrIS southwestern margin. Samples were collected in 2008 in the vicinity of a small (∼5 km2; see Section 3.3 for more details) land-terminating outlet glacier 125 km south of Jakobshavn Isbræ (‘N’ glacier; 68°02′ 34″ N, 50°16′08″ W), with a surface elevation range from ∼100 to ∼500 m (Fig. 1a and b). Supraglacial samples were collected from surface snow and a meltwater pond on the surface of ‘N’ glacier (Fig. 1c), and glacial ice samples at the margin of the glacier; subglacial samples were collected from the outflow stream (Fig. 1d) exiting the base of the glacier at two locations: the mouth where the outflow first emerged, and a downstream site ∼0.15 km from the outflow mouth. Groundwater samples were taken from the bank of the ‘N’ glacier proglacial stream and the flood plain of an adjacent larger glacier (‘M’ glacier) (Fig. 1e). Most samples were acquired while field personnel were on-site in late spring (16 May–1 June) and at the height of the summer melt season (10–17 July); automated instrumentation was used to measure and/or sample selected parameters between 2 June and 15 July. Constraints on end-member values were supplemented with supraglacial snow, ice and meltwater samples collected in July 2007 and 2008 within the ablation zone (980 m elevation) 70 km north of the primary study site (68°34′16″ N, 49°21′29″ W).
(a) Landsat image of the 2008 ice-margin field site. Sampling was concentrated at ‘N’ glacier and on the ‘M’ flood plain. (b) Air photo of the land-terminating outlet glacier (‘N’ glacier) focused upon in this study, with specific sample sites labeled. (c) Photo of the supraglacial pond sampled on the surface of ‘N’ glacier. (d) Photo of the mouth of the ‘N’ glacier outflow stream. (e) Photo of the ‘M’ flood-plain sample site.
3.1. Meteorological measurements
Local meteorological conditions were obtained using a HOBO U30-NRC weather station equipped with a tipping-bucket rain gauge and data logger installed at 100 m elevation and ∼1.5 m above the ground in the proglacial area in the vicinity of ‘N’ glacier. Shielded air temperatures and precipitation (0.2 mm resolution) were recorded every 5 min for the duration of the 2008 field season (16 May–17 July) with some gaps (73 min on day 138, 7 min on day 146, and 72.45 hours from day 151 to day 154). Hourly and daily moving averages were calculated for temperature, and daily total values were summed for precipitation. All times are reported in Greenland local time (GMT −3 h).
3.2. Discharge measurements
Stream discharge was measured at the ‘N’ glacier outflow stream using the velocity–area method and pressure transducers (HOBO U20 Water Level Logger and InSitu Level TROLL 300 Logger) at a location ∼0.15 km downstream of the glacier mouth. Stream velocities were measured using a manual flowmeter (General Oceanics Mechanical Flow-meter, model 2030R). A horizontal transect was established across the stream, and triplicate velocity measurements were taken at evenly spaced subsections (verticals) along the stream transect. The triplicate measurements were averaged to produce a single velocity at each vertical. Discharges for each stream subsection were calculated as the product of the subsection velocity and area (Reference DingmanDingman, 2002) and summed to obtain a total stream discharge (m3 s−1). Pressure transducers were used to continuously measure stream stage (depth). Water pressure was sensed in 10 min intervals from 31 May (19:00) to 16 July, and was converted to stream stage after correction for atmospheric pressure using a record sensed by an InSitu BaroTROLL Logger. A stage–discharge rating curve (r 2 = 0.76, p < 0.01) was developed using 12 discharge measurements. The rating curve was used to produce a continuous discharge record for the period with continuous stage measurements with point discharges only available from 19 to 27 May. Total meltwater discharge measured from ‘N’ glacier between 31 May and 16 July is 6.4 × 106 m3. The error associated with the discharge is estimated to be ±17% (following Reference DingmanDingman, 2002).
3.3. Catchment delineation
Lacking adequately resolved ice thickness and surface and basal topography for this region, we rely on interferometric synthetic aperture radar (InSAR)-derived ice velocity (personal communication from I. Joughin, 2011) to delineate the ‘N’ glacier surface area. We defined the catchment area at the downstream end to be bounded by the margins of the outlet glacier at the ice-sheet edge. We delimited the top of the catchment area where the background ice-sheet flow diverged from ‘N’ glacier. We were further limited in this approach by the coarse resolution (500 m) of the InSAR, so we defined divergence as separation of streamlines by one or more grid spacings. Using this method, we estimate the glacier length to be 5 km and the surface area to be 5 km2, which we defined as the catchment area for the purposes of evaluating surface meltwater input. To evaluate the reasonableness of this estimate, we calculated the total surface melt over this area required to match the cumulative discharge measured at the front of ‘N’ glacier. This calculation yielded a mean melt rate of ∼0.03 m d−1 (i.e. 1.28 m/46 days), a value well within the range of summer melt rates previously reported for the western margin of the GrIS (Reference BoxBox and others, 2006), thus providing an independent assessment of our catchment area.
4. End-Member Water Reservoirs and Isotope Measurements
Three conceptual end-member water reservoirs contributing to the bulk subglacial discharge from ‘N’ glacier were defined: (1) surface snow, (2) glacier ice and (3) delayed flow, where surface snow and glacier ice represent water sources, and delayed flow is a hydrological flow path. Samples were collected from the bulk discharge as well as each end-member (where possible) across the catchment. Samples were then analyzed for δ18O, δD, 222Rn and 7Be. In this study, delayed flow is operationally defined as water stored at the base on a timescale of days to weeks. These waters could consist of supraglacial waters stored at the base, basal-ice melt and groundwater. The timescale is dictated by the time required for radon (τ1/2 = 3.8 days) to approach secular equilibrium activity with its parent radium-226. At this activity, the production rate of radon is equal to its decay rate, and the radon content of waters stored at the bed is constant.
4.1. Stable water isotopes
Both seasonal-snow and glacial-ice δ18O and δD values vary with surface elevation, so samples were collected from this study to constrain these values. Samples include measurements from surface snow (∼300 m elevation), glacial ice (∼100, ∼300 and ∼1000 m elevation), groundwater (∼100 m elevation) and basal ice (∼100 m elevation). We then estimated the range of expected δ18O values over the ‘N’ glacier catchment (100–500 m surface elevation). Snow δ18O values were estimated using a 0.5‰ depletion per 100 m rise in elevation to account for the effect of altitude on δ18O (Reference Clark and FritzClark and Fritz, 1997). Glacial ice δ18O values were estimated using an empirical relationship defined by Reference Reeh and ThomsenReeh and Thomsen (1993) at a nearby GrIS location. The average difference (+2.5‰) between the measured ice δ18O at both 100 and 300 m and that calculated using the relationship defined by Reference Reeh and ThomsenReeh and Thomsen (1993) was used to correct the calculated δ18O ice value at 500 m. Bulk discharge samples for δ18O and δD measurements from the ‘N’ glacier outflow stream mouth and downstream sites were also collected. The discharge samples were collected at least daily from 18 May to 1 June and 10 to 16 July. Between 2 June and 9 July an ISCO 3700 autosampler (Teledyne Isco Inc.) was used to collect samples in intervals ranging from 1.5 to 4.5 days.
All samples were collected in acid-cleaned and sample-rinsed 250 or 1000 mL polypropylene bottles from which two 10 mL aliquots were taken for δ18O and δD analysis. Samples were frozen upon return to the laboratory until analysis. Thawed water samples were analyzed for δ18O and δD at the University of California Davis Stable Isotope Facility on a Laser Water Isotope Analyzer V2 (Los Gatos Research, Inc., Mountain View, CA, USA) with precisions of ≤0.3‰ for δ18O and ≤0.8‰ for δD.
4.2. Radon-222
Water samples were collected for end-member 222Rn activity from a supraglacial meltwater pond near the edge of ‘N’ glacier (22 and 31 May), from groundwater along the stream bank 0.15 km from the mouth of ‘N’ glacier (25 May) and from groundwater in the flood plain of ‘M’ glacier (28 May). Groundwater samples were taken at ∼0.4 m depth, using a stainless-steel drive point piezometer. Daily ‘N’ glacier outflow samples were collected on 18, 21–23, 27 and 29–31 May and 10–16 July from the mouth of ‘N’ glacier. Higher-resolution time-series samples (6 hour intervals) were taken on 31 May and 12 July. EC was also measured on-site using a Russell RF060C meter (Thermo Electron). Radon-222 samples were collected without head-space in glass 250 mL bottles, and were quantified using a RAD-7 continuous radon monitor (Durridge Inc.) (Reference Burnett and DulaiovaBurnett and Dulaiova, 2003). Typical RAD-7 uncertainties were 14%, with a range of 8–37% for the lowest measured 222Rn activities in this study. All samples were analyzed within 24 hours of collection. Results were corrected for radioactive decay between the time of collection and analysis and reported as an activity in disintegrations per minute per liter (dpm L−1). A model II (geometric mean) regression was used to compare the radon and EC data since both are measured (dependent) parameters with different units (Reference RickerRicker, 1973; Sofal and Rohlf, 1995).
In order to determine the maximum potential radon activities in saturated subglacial sediments, laboratory equilibration experiments were conducted using sediment collected from the proglacial area at the mouth of the ‘N’ glacier outflow (n = 1) and downstream (n = 2). Approximately 100 g of wet sediment were incubated with ∼0.5 L radium-free water in sealed 1 L high-density polyethelyene bottles, following methods described by Reference Dulaiova, Gonneea, Henderson and CharetteDulaiova and others (2008). The sediment was incubated for at least 3 weeks and the radon was subsequently quantified via an alpha scintillation technique. Each sediment equilibration sample was analyzed twice. The measured radon activities (dpm g−1) in the wet sediment were converted to pore-water radon activities (dpm L−1) using a wet bulk density of 2.3 g cm−3 and a porosity of 0.2 (Reference Dulaiova, Gonneea, Henderson and CharetteDulaiova and others, 2008).
4.3. Beryllium-7
Water samples (∼190 L) for measuring 7Be activity were collected from a supraglacial meltwater pond near the edge of ‘N’ glacier (22 May) and from a supraglacial meltwater stream at the inland ice-sheet site (20 July). Two ‘N’ glacier discharge samples were collected for 7Be (21 May and 11 July). Water was collected in a large plastic container and processed on-site. The supraglacial samples were not filtered prior to collection due to the lack of particles in these waters; however, the ‘N’ glacier outflow stream sample was filtered through a 10 mm Hytrex II cartridge prior to collection. A 1 mL aliquot of stable 9Be (10 000 ppm) was added as a yield monitor, and iron oxide (Fe(OH)3) fibers were used to pre-concentrate both the 7Be and 9Be from the water sample (Reference Andrews, Hartin and BuesselerAndrews and others, 2008). Periodic aliquots (20 mL) of the fiber column filtrate were taken and subsequently analyzed on an inductively coupled plasma mass spectrometer (ICP-MS) at the Woods Hole Oceanographic Institution (WHOI) to determine the collection efficiencies for 7Be (Reference Andrews, Hartin and BuesselerAndrews and others, 2008). The fibers were combusted at 820°C for 16 hours and the ash was analyzed for 7Be via gamma spectroscopy (Reference Andrews, Hartin and BuesselerAndrews and others, 2008). Samples were counted for 2 days, and corrected for decay since the time of collection. Beryllium recovery on the fibers averaged 75%. Results are presented as activities in units of dpm L−1.
5. Model Descriptions
5.1. Isotope-mixing model
A multi-component isotope-mixing model using the stable water isotope values and radon measurements described above was constructed to quantify the relative fraction of flow contributed by each end-member water source to the total discharge exiting ‘N’ glacier using
where f is the fraction of flow and subscripts indicate the following: 1 is the snow end-member, 2 is the glacial ice end-member, 3 is the delayed-flow end-member and 4 is the ‘N’ glacier outflow. The end-member water source and ‘N’ glacier outflow stream isotope values (Table 1) were used to solve the system of equations using singular value decomposition. The radon content in the snow and glacial ice end-members (222Rn1 and 222Rn2) was set to zero (negligible in situ source of 222Rn). To initially solve the model, the highest radon activity measured in the ‘N’ outflow stream (210 dpm L−1) was operationally defined as 222Rn3, effectively normalizing the entire dataset to this maximum concentration. In this study the sediment radon flux to the delayed-flow reservoir is assumed to be a continuous, steady-state process, so a loss term due to radioactive decay is not included in the radon end-member mixing equation.
Isotope tracer values used to initially solve the end-member mixing model equations (Equations (1–4)). For the delayed-flow waters a basal ice sample collected at 100 m elevation was used for the δ18O and δD ratios, while the maximum ‘N’ glacier outflow radon activity was used as the 222 Rn end-member. The radon activity of the surface snow and glacial ice reservoirs was set to zero
5.2. Mixing-model sensitivity analysis
A sensitivity analysis employing a range of end-member values (Table 2) was conducted to put envelopes of uncertainty on the model fraction results (f 1, f 2, f 3) by simultaneously varying each of the end-member water reservoir values across a range of reasonably determined limits. The range of surface snow δ18O values utilized was chosen to incorporate a maximum isotopic depletion from the original snow during metamorphism and melting (Reference Taylor, Feng, Kirchner, Osterhuber, Klaue and RenshawTaylor and others, 2001). Though the mean snowmelt will become progressively enriched throughout the summer melt season due to the early removal of isotopically light water (Reference Cooper, Kendall and McDonnellCooper, 1998; Reference Taylor, Feng, Kirchner, Osterhuber, Klaue and RenshawTaylor and others, 2001), this isotopic enrichment is difficult to predict without additional samples of the total snowpack oxygen isotope signature. However, the model results will be most affected by the potential for depleted snow end-member values, which encroach on the glacial ice δ18O values. The glacial ice δ18O values were chosen to reflect the range of ice isotopic values across the elevation range of the ‘N’ glacier catchment that could potentially contribute to the outflow waters. The range of delayed-flow δ18O values employed in the sensitivity analysis encapsulated the δ18O content of the groundwater and basal ice at ‘N’ glacier. Corresponding dD ranges for each end-member water source were calculated directly from these δ18O values using the global meteoric waterline (δD = 8δ18O + 10). The range of potential 222Rn activities in the delayed-flow end-member was defined using the maximum activity measured in the ‘N’ outflow stream as the lower bound and the lowest groundwater radon activity measured in this study as an upper bound. A better estimate of the 222Rn activity in the delayed-flow waters could be attained by sampling outflow waters during the early season before they are diluted by any surface input. This could be achieved with standard automated continuous radon monitors that measure the 222Rn activity of the outflow stream (e.g. Reference Dulaiova, Peterson, Burnett and Lane-SmithDulaiova and others, 2005; Reference Schmidt, Schlueter, Melles and SchubertSchmidt and others, 2008).
Range of isotope tracer values used in end-member mixing model sensitivity analysis
5.3. Transit time model
The surface snow and glacial ice fraction results (f 1 and f 2) from the isotope model were used to estimate a transit time for snowmelt from the surface to its exit at the glacier front using 7Be. 7Be is continuously produced in the atmosphere and is deposited to the surface environment (e.g. the ice-sheet surface) via wet (precipitation events) and dry (aerosol) deposition (Reference Nimz, Kendall and McDonnellNimz, 1998). Its unique atmospheric source in combination with its short half-life (τ 1/2 = 53.3 days) suggests that 7Be should be present in the surface snow end-member, while occurring below detection levels in glacial ice and delayed-flow waters that are older than ∼300 days. Thus, assuming constant production on the surface and no 7Be in the delayed-flow reservoir (f 3), the 7Be activity in the outflow stream can be described using a steady-state model:
where t is time (days), A is the decay constant for 7Be (0.013 d−1), 7Be1, t0 is the 7Be activity of the surface snow end-member water source, 7Be2 is the 7Be activity of the glacial ice, and 7Be4eff, t is the effective 7Be activity in the ‘N’ glacier outflow stream. 7Be4eff, t accounts for any scavenging of 7Be onto subglacial particles and is calculated as the sum of the 7Be in the ‘N’ glacier outflow stream (7Bedissolved) and the 7Be scavenged by particles in the subglacial environment (7Beparticulate). The former was measured in the ‘N’ outflow stream on 21 May, and the latter was estimated using
where K d is particle-water coefficient ((7Be/mass of particles)/(7Be/mass of water)) and C p is suspended sediment concentration (mg mL−1) measured in May. A K d value of 5000 was assumed, which is within the range of previously published 7Be K d values (103–104) (Reference Olsen, Larsen, Lowry, Cutshall and NicholsOlsen and others, 1986). Rearranging Equation (5) to solve for time yields
Here time is defined as the transit time from the glacier surface to the front since 7Be is only produced on the surface, and we have accounted for any sinks (scavenging) in the subglacial environment.
6. Results and Discussion
6.1. Climatology and discharge
‘N’ glacier discharge was highly sensitive to air-temperature fluctuations, with the two clearly co-varying throughout the melt season (Fig. 2a and b). We did not capture the melt onset, as there was already meltwater discharge on our arrival at the field site on 19 May. Nonetheless, a cold period towards the beginning of our record (days 142–144), characterized by subfreezing air temperatures, reduced discharge almost to zero. Following this cold period, air temperature and discharge increased from mid-May to early June as the melt season progressed. Daily-average temperatures at the ice edge then generally remained above 6°C for the remainder of the study period. Daily-average discharge also stabilized around 1.6–1.8 m3 s−1 until late July, when a cold period precipitated a drop in discharge back towards early-season values. We also observed a strong diurnal cycle in subglacial stream discharge that was highly responsive to, although offset from, the insolation-driven diurnal temperature cycle. Peak discharge lagged peak air temperature by an average of 2.4 hours (range 1–5.6 hours), while the average lag in minimum discharge was 2.7 hours (range 0.35– 4.7 hours) from the minimum temperature. The mean diurnal amplitude in the temperature and discharge records was 7.3°C and 0.39 m3 s−1, with maximum and minimum daily discharges occurring between 15:00 and 21:00 and 4:00 and 8:00, respectively.
(a) Plot of 3 hour (thin black line) and daily (thick black line) average air temperatures, with total daily rainfall in the gray bars. (b) Three-hour (thin black line) and daily (thick black line) average discharge at‘N’ glacier. (c) Radon activities. (d) δ18O content (black dots) and five-point moving average (black line) in the ‘N’ glacier outflow stream. The discharge record is confined to point measurements from 19 to 31 May. Discontinuous lines in the temperature and discharge records reflect gaps in the data.
6.2. Radon as a tracer for delayed-flow waters
We did not detect any radon in a sample from a supraglacial meltwater pond on the surface of ‘N’ glacier, thus confirming our assumption that surface snow and glacial ice would be devoid of radon due to negligible sediment inventories on the ice-sheet surface. Conversely, as expected, radon activities in the groundwater samples were very high (1626 ± 48.8 and 2750 ± 63.3 dpm L−1 ). These activities were consistent with the laboratory-derived pore-water radon activities (range 1285 ± 43.3 to 3045 ± 132 dpm L−1), thus indicating that the groundwater samples collected in this study represent saturated flow. We observed seasonal differences in the amount of radon detected in the ‘N’ glacier outflow stream waters. The mean activity of the May samples was much higher (75.6 dpm L−1) than that of the July samples (25.4 dpm L−1). Furthermore, we observed a greater range in the radon activities of the outflow stream in May (16.9–210 dpm L−1) compared to July (10.4–35.7 dpm L−1) (Fig. 2c). Thus, during at least some periods in the late spring, the‘N’stream outflow waters had high radon activities, whereas at the height of the summer melt season the outflow waters had universally low radon activities. For comparison, the open ocean has an average radon activity of ∼0.01 dpm L−1 (Reference Broecker and PengBroecker and Peng, 1982). In groundwater, radon activities can vary greatly, but, in general, activities range from hundreds to thousands of dpm L−1 (Reference Charette, Moore, Burnett, Krishnaswami and Kirk CochranCharette and others, 2008 and references therein). We were not able to resolve a diurnal cycle in the ‘N’ outflow radon activity during the 12 July high-resolution time series due to the low activities measured throughout this day (27.3–31.8 dpm L−1). Earlier in the season, however, we found that the maximum daily radon activity on 31 May (75.5 dpm L−1 at 06:45) occurs within the period when daily discharge was at a minimum (although we lack continuous discharge measurements during that time). Moreover, we observed the lowest radon activity (47.5 dpm L−1) at 18:25, when daily discharge was at a maximum. These preliminary data indicate that radon may be useful in resolving the diurnal contribution of delayed-flow waters to total outflow in the early season, but this application requires more frequent sampling. Thus, we limit further discussion of radon to seasonal trends.
Radon in water that is physically decoupled from its sediment or rock source is subject to decay on a timescale determined by its half-life (Reference Kraemer, Genereux, Kendall and McDonnellKraemer and Genereux, 1998). Thus, subglacial outflow radon activities in line with published groundwater values require substantial steady-state sediment–water interaction. To quantify the potential for suspended sediment to explain the observed radon values, we collected replicate unfiltered samples from the mouth of the ‘N’ outflow; one sample was analyzed immediately, while the other was measured after one radon half-life. The decay-corrected activities of the samples were the same, indicating that 226Ra in the stored sample sediment was not a significant source of radon. We therefore conclude that the presence of radon in bulk outflow waters necessitates some delayed-flow component that has had substantial interaction with the bed. As further evidence of this idea, the regression (model II, geometric mean) of radon and EC was significant (r 2 = 0.87, p < 0.01) (Fig. 3). However, unlike EC, radon is not subject to the dissolution chemistries of a wide range of solutes and thus can potentially be utilized for quantitative hydrograph separation.
Model II regression (geometric mean) of EC and radon activity in the ‘N’ glacier outflow stream waters.
6.3. Radon evasion in subglacial channels
Although radon is water-soluble, the radon will partition into the air phase in an air–water system (Reference Kraemer, Genereux, Kendall and McDonnellKraemer and Genereux, 1998). Loss due to evasion is a function of temperature and the amount of radon present in the air, so evasion could be problematic late in the melt season when subglacial channels may not be entirely water-full. Moreover, water flow in the subglacial channels is often faster and more turbulent at the peak of the summer melt season, which may enhance gas exchange loss (Reference Kraemer, Genereux, Kendall and McDonnellKraemer and Genereux, 1998). Thus, if evasion were the dominant process influencing radon in the ‘N’ glacier outflow stream we would expect to find lower radon activities toward the end of a falling discharge limb. Instead we observed increases in radon during times when discharge was decreasing and subglacial channels were less full (e.g. days 192–195), and the largest radon activities occurred at the discharge minimum (Table 3). Thus, dilution of the delayed-flow waters with radon-free surface input appears to have had the greatest effect on radon values in the ‘N’ outflow stream.
222Rn activities in the ‘N’ glacier outflow stream on falling and rising discharge limbs from 10 to 16 July 2008
6.4. Water isotopes as a tracer in a GrIS outlet glacier
We found distinct δ18O values for the surface-snow (−11 to −13‰) and glacial-ice (−26 to −30‰) end-members measured across our study region (Table 4). This difference means that the snow and glacial ice reservoir δ18O values do not overlap and thus are useful as passive flow tracers. Conversely, our‘basal-ice’ (ice collected at 100 m elevation at the glacier margin) sample values (−29.6‰) were not sufficiently isotopically distinct from glacial ice values estimated across the surface of the catchment (−25.5 to −28.2‰) to separate the delayed-flow water source from the glacier ice reservoir. Thus it was necessary to delineate these reservoirs from each other with the radon end-member mixing equation. An additional groundwater sample near the front of ‘N’ glacier also possessed a depleted δ18O signature (−27.8‰), suggesting that groundwater in this region is derived primarily from glacial ice melt. Although we limit our discussion to the δ18O values, the trends observed in the oxygen isotope values are also applicable to the deuterium data, since δ18O and δD co-vary on a global scale (Reference CraigCraig, 1961).
Measured and estimated(*) δ18O ratios of surface-snow, glacial-ice and groundwater samples
6.5. Isotope-mixing model
The end-member mixing equations used in this study assume a simplified drainage system limited to three conceptual end-member water sources: (1) surface snow, (2) glacial ice and (3) delayed-flow waters. This was necessary to determine the applicability of radon as a hydrological tracer in a glacial setting, and to produce an initial chemical mixing model, although we recognize that we have oversimplified the subglacial drainage system by categorizing glacial flow components into three broad water source reservoirs (Reference Sharp, Brown, Tranter, Willis and HubbardSharp and others, 1995).
Model results showed that the snowmelt and delayed-flow waters comprised a greater fraction of the total outflow in May (Fig. 4a) than in July (Fig. 4b). In May, delayed flow dominated the discharge (mean 41%), followed by nearly equal contributions from surface snowmelt (mean 23%) and glacial ice melt (mean 26%). In July, however, the mean fractional contributions from the surface snowmelt and delayed-flow reservoirs decreased to 6% and 12%, respectively, while the mean glacial ice contribution rose to 82%. This finding was likely due to the removal of seasonal snow from the glacier surface by this time, and dilution of delayed-flow reservoirs with increased glacial ice melt. Scaling the model results with the measured discharge allowed us to compare the discharge contribution of surface snow, glacial ice and delayed flow to the total ‘N’ stream discharge from 18 May to 1 June (Fig. 4c) and 11 to 17 July (Fig. 4d). The average snow component (n = 11) of total discharge decreased by more than half from May (mean 0.17 ± 0.04 m3 s−1; 1 standard error) to July (mean 0.07 ± 0.01 m3 s−1), whereas the average glacial ice component (n = 11) more than doubled from May (mean 0.43 ± 0.10 m3 s−1) to July (mean 1.1 ± 0.09 m3 s−1). By comparison, the average delayed-flow component (n = 11) remained a relatively constant contribution between May (0.19 ± 0.05 m3 s−1) and July (0.16 ± 0.01 m3 s−1).
Stacked bar plots of isotope-mixing model solutions for the fractions of surface snow, glacial ice and delayed-flow waters comprising the ‘N’ glacier outflow stream waters from (a) 18 May to 1 June and (b) 11 to 17 July, and scaled contributions from each reservoir in (c) May and (d) July.
The results of our isotope-mixing model were a direct consequence of the shifts we observed in the 222Rn activities and δ18O and δD ratios of the ‘N’ outflow stream composition from May to July. The highest radon activities were found during times of lowest discharge (Fig. 5a). During the 3 day subfreezing period in May when discharge dropped to near zero, radon activities in the ‘N’ outflow stream peaked at >100 dpm L−1 (Fig. 5a). In July, even though the radon activities were overall much lower than in May, elevated radon (35.7 dpm L−1) coincides with a prominent drop in discharge on day 195 (Fig. 5b). This behavior can best be explained by varying levels of dilution of the delayed-flow waters with a supraglacial water source devoid of radon. This reasoning is consistent with a seasonal evolution of the subglacial drainage structure from a distributed system characterized by chemically enriched outflow waters to a channelized system that facilitates rapid transit of dilute glacial ice melt.
(a) Radon activities (± one standard error) in the ‘N’ glacier outflow stream plotted against daily average discharge, and (b) 3 hour average discharge and measured radon activities in the ‘N’ glacier outflow stream from 11 to 17 July. Discontinuous lines in the discharge record reflect gaps in the data.
The difference between the stable-isotope signatures of the snow and ice reservoirs at ‘N’ glacier is sufficiently large that a change in δ18O runoff composition can likely be attributed to a water source change. In late spring when dischargeis low, and snowmelt feeds a predominantly distributed subglacial drainage system, we measured enriched δ18O values in the ‘N’ outflow stream, compared to later in the season (Fig. 2d). Indeed, the most enriched values occurred over the 3 day span that discharge dropped to near zero. Subsequently there is a decrease in the stable-isotope signature of the ‘N’ glacier discharge as the melt season progresses from late spring to the summer. This trend is consistent with a seasonal shift in water source reservoirs from a snow and ice contribution in late spring to a purely glacial ice contribution at the height of the summer melt season.
In addition to the overall seasonal decline, our δ18O record exhibited higher-frequency variability, suggesting that changes in meltwater source contributions and/or drainage system evolution may also have occurred during synoptic-scale events. Although we lacked the temporal resolution required to explore this variability in full, one such event late in the melt season is reasonably well resolved. On day 195 there was a notable spike in the ‘N’ glacier subglacial stream δ18O values above the late-season mean coincident with a prominent drop in air temperature and discharge, and an increase in radon activity (Fig. 2). One possible explanation for this event is cooling air temperatures across the glacier surface leading to a decrease in δ18O-depleted glacier ice melt. This decrease in total surface meltwater input to the subglacial resulted in a relatively higher base flow contribution (characterized by residual stored δ18O-enriched snowmelt) to the bulk runoff during this event. Another possible explanation for this isotopic excursion is a rainfall event (0.2 mm) recorded that day (Fig. 2a) which would also yield an enriched δ18O signature. Larger rainfall events in our record (e.g. 8.4 mm on day 166), however, did not correspond to enriched δ18O runoff values. Furthermore, rainfall events should increase stream discharge, whereas we observed a decrease in stream discharge during this event.
6.6. Mixing-model sensitivity analysis
The sensitivity analysis revealed that for the entire dataset (May and July) the contribution from the surface snow reservoir to the total outflow varied from a mean maximum of 26% to a mean minimum of 0.9%. Similarly, the delayed-flow fraction varied from a mean maximum of 26% to a slightly higher mean minimum of 3.4%. Not surprisingly, the glacial ice fraction exhibited the highest mean maximum and minimum values, varying from 97% to 49%. Results of the sensitivity analysis are also displayed as the maximum and minimum flow contributions (m3 s−1) from the surface snow, glacial ice and delayed-flow reservoirs from May (Fig. 6a) and July (Fig. 6b). In order to identify meaningful estimates, we constrained the sensitivity analysis so that no flow contribution was permitted to fluctuate below zero. When discharge was very low, we were not able to differentiate the flow contribution (m3 s−1) from the different component reservoirs accurately (Fig. 6a). Additionally, though we were able to drive the snow contribution to zero, there was always a delayed-flow component that is diluted by an increasing ice component throughout the season. This analysis illustrates that we are currently able to determine the flow contribution from each of the defined water source reservoirs within an absolute maximum and minimum value. Further improvements to the flow estimates would benefit most from better characterization of the delayed-flow radon end-member.
Sensitivity analysis illustrating maximum (solid lines) and minimum (dashed lines) flow contributions from the surface snow, glacial ice and delayed-flow water sources from (a) 18 May to 1 June and (b) 11 to 17 July.
7. Transit Time Estimates
We observed detectable 7Be activity (7.7 dpm L−1) in the ‘N’ supraglacial pond in May, which we use as an analogue for the potential 7Be activity of surface snow in this study. This activity is within the range of previously reported 7Be activities in fresh snow at Summit, Greenland (2.67–76.5 dpm L−1), but below the reported median (15.3 dpm L−1) (Reference DibbDibb, 1990). The wide variability in reported fresh-snow 7Be activities likely reflects atmospheric inventory depletion and wet deposition-related dilution effects of precipitation-event frequency and duration (Reference Nimz, Kendall and McDonnellNimz, 1998). For example, high 7Be activities may result from a short snowfall event following a period of minimal precipitation. Conversely, lower snow activities may be explained by a relatively large snowfall event following a series of recent precipitation events. Our sample may have had lower 7Be activity compared to the fresh snow collected at Summit because it represents a composite of fresh and older snow on the surface of ‘N’ glacier. Comparatively, meltwater derived from recent glacial ice melt measured at the inland site in July had an extremely low 7Be activity (0.04 dpm L−1), an indication that its original 7Be inventory had been lost via decay. Thus, 7Be in outflow water can only be derived from a young supraglacial source that has originated at the surface <1 year before. On 21 May the 7Be activity of the ‘N’ glacier outflow stream was 1.05 dpm L−1, and on 11 July the 7Be activity was 0.03 dpm L−1. The low July 7Be signal was similar to the recent glacial ice-melt 7Be signal, and most likely represented a switch in end-member contribution from snowmelt to ice melt. The May value, however, was consistent with a hydrological connection between surface melt and subsurface discharge at this point in the season. However, we cannot rule out the possibility that the May 7Be signal in the ‘N’ glacier outflow stream resulted from the release of supraglacial waters that had been stored at the bed for <300 days or in basal crevasses (Reference Harper, Bradford, Humphrey and MeierbachtolHarper and others, 2010).
We used the fractions from the isotope-mixing model for the surface snow and glacial ice contributions on 21 May to solve for a transit time (Equation (7)). Since 7Be is particle-reactive (Reference Hawley, Robbins and EadieHawley and others, 1986; Reference Olsen, Larsen, Lowry, Cutshall and NicholsOlsen and others, 1986) and thus could have adsorbed to solids during travel through the subglacial environment, we also included a correction that describes scavenging of 7Be onto particles in the subglacial environment (Equation (6)). Given these assumptions, we estimated that supraglacial waters took ∼7.5 days to travel from the surface to the glacier mouth. For comparison, distributed drainage-system transit times (estimated using velocities from dye-tracing experiments in Reference Nienow, Sharp and WillisNienow and others, 1998) at Haut Glacier d’Arolla, Switzerland, which has a similar catchment size to ‘N’ glacier, range from ∼6 to <1 day(s) for a 5 km flow path. However, we should note that our calculated time estimate depends on the accuracy of our assumed partition coefficient, K d, and the surface 7Be activity. Nonetheless, 7Be may hold promise as a tracer for snowmelt in early-season distributed drainage systems with similar or longer transit times.
8. Synthesis
Results from our multi-component isotope-mixing model provide broad separation of water reservoir contributions, thus providing a potential new direction in the application of chemical mixing models to study glacier hydrology. Based on the results from this study, we furthermore suggest that these methods could be successfully scaled up to investigate the subglacial hydrology of the much larger outlet glaciers that drain the bulk of the GrIS. Some practical and technological challenges remain to be solved in regularly sampling discharge and radon at large land-terminating glaciers (e.g. in flood plains, large channels or braided river environments) and in large marine-terminating glaciers (tidewater environments). This effort would also require comprehensive sampling of the glacial ice end-member water isotope values across each of these larger catchments.
Focusing our study on a small land-terminating glacier on the southwestern margin of the GrIS, we show that there is a relatively constant and chemically enriched delayed (basal) flow component present throughout the melt season. These delayed-flow waters comprise a greater fraction of the total discharge in May compared to July, and are diluted first by snowmelt and then by increasing amounts of rapidly flowing ice melt as the season progresses. In alpine glaciers, chemically enriched delayed-flow waters (e.g. snowmelt, basal melt, groundwater) are characteristic of distributed drainage systems, which transmit meltwater slowly through the glacier via a hydraulically inefficient network. As the snowline retreats and surface meltwater input to the bed increases, the subglacial drainage system structure evolves to a channelized drainage system, which can more efficiently export the surface glacial ice melt. Though such seasonal subglacial drainage evolution is well documented in alpine systems (Reference Hubbard and NienowHubbard and Nienow, 1997; Reference Nienow, Sharp and WillisNienow and others, 1998; Reference Cuffey and PatersonCuffey and Paterson, 2010), its existence under land-terminating sectors of the GrIS has only recently been hypothesized (Reference Shepherd, Hubbard, Nienow, McMillan and JoughinShepherd and others, 2009; Reference Bartholomew, Nienow, Mair, Hubbard, King and SoleBartholomew and others, 2010) and has not been directly observed. The findings from this study offer a hydrochemical line of evidence for this hypothesis, albeit at a comparatively much smaller outlet glacier, and bolster the need to incorporate these dynamically significant subglacial processes into GrIS modeling efforts.
From a biogeochemical perspective, knowledge of the seasonal controls on end-member water-source contributors to bulk discharge provides greater understanding of the potential for high temporal variability of carbon, nutrient and metal export from subglacial environments to downstream marine ecosystems. Previous studies have suggested that water draining a distributed drainage system contains much greater concentrations of these biogeochemically important species, compared to the waters draining a channelized system (Reference Tranter, Skidmore and WadhamTranter and others, 2005). Thus, total annual flux calculations of carbon, nutrient and metal export require knowledge of the base flow (m3 a−1) exiting a glacier. Our isotope-mixing model shows promise at being able to provide reasonable quantitative estimates of snow, ice and delayed-flow components comprising bulk meltwater discharge from a land-terminating GrIS glacier. These flow estimates can be used as a first-order approximation of base flow emanating from similar catchments around the GrIS throughout the year.
Acknowledgements
This research was supported by the WHOI Clark Arctic Research Initiative (E.B.K., S.B.D., M.A.C.), the WHOI Ocean Ventures Fund (M.P.B.), the US National Science Foundation ARC-05200077 (S.B.D.), NASA (S.B.D.), the Natural Sciences and Engineering Research Council of Canada (M.P.B.) and the WHOI Climate Change Institute (M.P.B.). We acknowledge I. Joughin for kindly providing InSAR velocity data. We are very grateful to M. Behn and D. Glover for assistance with data analysis, to M. Sharp and K. Longnecker for comments and advice that improved the manuscript, and to B. Gready, M. Behn, I. Joughin, M. Evans, A. Criscitello and R. Harris for assistance in the field.
