Clay mineral precipitation and low silica in glacier meltwaters explored through reaction-path modelling

Abstract

The subglacial chemical weathering environment is largely controlled by low temperatures and the presence of freshly comminuted minerals with a high surface area. These characteristics are believed to promote dissolution processes that give rise to low silica and high Ca²⁺ fluxes emanating from glacierized basins. We test an alternative hypothesis, that mineral precipitation reactions in the subglacial environment play an equally important role in controlling the water chemistry in glacierized basins. We analyze borehole and proglacial water chemistry from a subarctic polythermal glacier, complemented by mineral XRD analysis of suspended sediment, till and bedrock samples. In conjunction with a thermodynamic analysis of the water and mineral chemistry, we use reaction-path modelling to study the chemical enrichment of water through the glacier system. We find that the high pH of the subglacial environment is conducive to secondary mineral precipitation, and that it is not possible to balance the water chemistry using dissolution reactions alone. We show that low silica can be explained by standard weathering reactions without having to invoke mineral-leaching reactions. Our results suggest that subglacial weathering intensity may be significantly underestimated if the production of secondary minerals is not considered.

Introduction

Glacierized basins have a unique chemical weathering environment, characterized by low temperatures, variable access of meltwater and atmospheric gases to the glacier bed, and high amounts of rock comminution, thereby providing fresh and abundant mineral surfaces. This has motivated much research into how chemical weathering differs between glacierized and non-glacierized basins (e.g. Anderson and others 1997; Lafreniere and Sharp, 2005; Pogge von Strandmann and others, 2006; Wimpenny and others, 2011), and the relevant impacts on atmospheric CO₂ drawdown (e.g. Sharp and others, 1995a; Hodson and others, 2000; Tranter and others, 2002a). More recent work has shifted the focus of subglacial weathering to consider the importance of alternative proton sources for mineral dissolution through microbially mediated pyrite oxidation and organic carbon sources, because atmospheric CO₂ and O₂ supply are often limited at the glacier bed (e.g. Sharp and others, 1999; Tranter and others, 2002b; Montross and others, 2013). Yet there are still a number of unanswered questions pertaining to basic chemical mass-balance problems in the subglacial environment. For example, there is considerable uncertainty in the extent to which weathering is controlled by the leaching of fresh mineral surfaces versus the weathering of mineral surfaces undergoing stoichiometric or incongruent dissolution. Further uncertainty exists in trying to relate the water chemistry at the terminus to chemical reactions in the distributed system, where waters may have residence times that are orders of magnitude longer than in a channelized system.

Inferences of subglacial chemical weathering processes have been largely based on studies of proglacial waters that are significantly diluted by surface melt, and are therefore quite far from saturation for most minerals. These studies often conclude that subglacial weathering is governed by mineral dissolution processes (e.g. mineral leaching reactions, the incongruent dissolution of feldspars to kaolinite and the stoichiometric dissolution of oxides, carbonates and other framework silicates (e.g. Anderson, 2005; Pogge von Strandmann and others, 2006; Wimpenny and others, 2011; Ryu and Jacobson, 2012)). During mineral leaching and incongruent dissolution, the molar proportions of ions in solution do not reflect the molar proportions of elements in the mineral from which the ions were weathered, while the opposite is true for congruent or stoichiometric weathering. As a consequence of emphasizing dissolution reactions, the effect of secondary mineral precipitation and clay formation beyond kaolinite is generally not accounted for in the overall chemical mass balance (e.g. Sharp and others, 1995a; Hodson and others, 2000; Anderson, 2005). As such, chemical anomalies in glacierized environments, namely high Ca²⁺ and low silica fluxes, have been interpreted solely in the context of mineral dissolution reactions.

Previous observations of subglacial water composition, through borehole sampling at Haut Glacier d’Arolla, Switzerland, have hinted at the possibility of secondary mineral precipitation (Tranter and others, 2002b), while results from inverse mass-balance modelling of Haut Glacier d’Arolla waters indicate the need to precipitate smectite to balance the water chemistry (Mitchell and Brown, 2008). Many secondary mineral precipitates have been found on recently deglaciated bedrock as silica (Hallet, 1975), calcite (Hallet, 1976; Sharp and others, 1990) and calcite-cemented siltskins (Carter and others, 2003), while authigenic nanoparticles of ferrihydrite and geothite have been found to coat subglacial sediments (Raiswell and others, 2009; Hawkings and others, 2014). In this study, we aim to test the hypothesis that secondary mineral precipitation plays an important role in controlling the chemical fluxes from glacierized basins. We analyze borehole and proglacial water from a polythermal surge-type glacier underlain solely by granodiorite, and supplement these data with mineralogical analyses of the bedrock, till collected at the glacier terminus and sediments suspended in proglacial and borehole waters. We forward model the water composition using reaction-path modelling, which uses mineral saturation, reaction mechanism, mineral surface area and water temperature to compute rates of mineral dissolution and precipitation. Analysis and modelling of the hydrochemistry is divided into three conceptual stages: (1) initial water/rock reactions; (2) mineral precipitation induced by long water/rock contact times and/or basal freeze-on and (3) post-mixing reactions.

Field Site

The study site is an unnamed glacier in the Donjek Range of the St Elias Mountains, Yukon, Canada (Fig. 1), hereafter referred to as South Glacier. The Donjek Range lies within the central belt of the Kaskawulsh Group in the Alexander terrane (Wheeler, 1963; Campbell and Dodds, 1978; Israel and Cobbett, 2008). Field mapping at the ice margin indicates that the glacier is underlain solely by the Shield pluton granodiorite, dated with K–Ar biotite ages of 133 ± 5 Ma, and hornblende K–Ar ages at 144 ± 6 Ma (Dodds and Campbell, 1988). The northeast ridge of the basin is flanked by a nonconformity, where the pluton imprints a narrow contact aureole on the surrounding graphitic and siliceous argillites and sandstones of late Paleozoic age (Dodds and Campbell, 1988). Small quantities of rock around the contact may be wasting into the glacier basin, but this likely contributes negligibly to the mineralogy observed in the proglacial stream and till. A 2m wide dike of ortho-amphibolite cuts the pluton at the eastern glacier margin, but this dike is not observed on the western margin.

Fig. 1. DigitalGlobe^® image (9 September 2012) of South Glacier field site. The glacier is ∼5.3 km², and is within a basin having a total area up-valley from the terminus of 9.7 km². The glacier spans elevations 1970–2960 m a.s.l. Symbols denote sampling locations for bedrock, till, borehole (BH), supraglacial and proglacial streams.

South Glacier has a total area of of 5.3 km², with a center-line length of ∼5 km and an elevation of 1970–2960 m a.s.l. The accumulation area is temperate throughout the ice depth, while the ablation area is characterized by a surface layer of cold ice that increases in thickness down-glacier (Wilson and others, 2013), with the potential for a cold bed in the lower ablation area (Wilson and others, 2014). In the ablation area, the glacier margins often border debris-covered ice that grades into ice-cored moraines. Annually averaged ice surface velocities of ∼30 m a⁻¹ from 2006–14 show that the upper ablation area is exceeding the expected balance velocity, suggesting that the glacier is currently undergoing a slow surge (Flowers and others, 2011); the last observed surge occurred in 1986/87. In this region of the glacier, sliding speeds account for 50–100% of the total ice surface velocity (De Paoli and Flowers, 2009). Borehole water pressure data from 2008–11 show that faster sliding is observed in conjunction with high basal water pressures, which may occur both in and out of phase with peak diurnal supraglacial melt (Schoof and others, 2014).

Data Collection and Analysis

Sample collection

Hydrochemical and suspended sediment samples were collected from the proglacial stream ∼100 m downstream from the glacier terminus, and from boreholes that were drilled to the base of the glacier in the central ablation area (Fig. 1). The sampling schedule is shown in Figure 2. Three water samples were collected from 28 to 30 May 2013, but the bulk of the sampling was carried out in July 2013, with a full 24 hour sampling period on 23/24 July. Suspended sediments were collected on 28 May and throughout July. Subglacial water was sampled in boreholes that were drilled using a hot-water drill (as described by Schoof and others, 2014), and the samples were retrieved using a Niskin sampler, following the design of Blake and Clarke (1991). Borehole samples were measured for pH, alkalinity, dissolved oxygen (DO), electrical conductivity (EC) and cation and anion composition on 26 July, and for suspended sediment mineralogy on 7 and 26 July.

Fig. 2. Water discharge measurements (red squares) and interpolation (blue curve) for July 2013. Sampling schedule for borehole and proglacial water and suspended sediment (ss) shown by vertical bars at bottom. The gray bars indicate the timing and extent of the transitions between rating curves (RC).

All samples were filtered within 30 min of collection, through a 0.45 μm cellulose acetate membrane, using a Nalgene vacuum filtration unit. Highly turbid samples were first filtered through a 0.80 μm mixed cellulose ester membrane. For each sample, the filtration unit was rinsed with distilled water and then with filtrate. The initial water volume was measured in the upper chamber using a bubble level. To estimate the suspended sediment concentration (SSC), used filters were placed in pre-weighed glass vials with tweezers, and later dried at 60°C for a minimum of 24 hours, at which point they were weighed every hour to ensure no further weight loss. Upon filtration, cation and anion samples were stored in 125 mL HDPE Nalgene bottles at 0–4°C for up to 4 months before analysis. Cation samples were acidified with one drop of 1M nitric acid. Suspended sediment samples from the proglacial stream were obtained by collecting ∼15 L of water, which was then sealed and left to settle. After 3 days, ∼14.5 L were siphoned off and the sediment and remaining water were stored in 500 mL HDPE bottles with no remaining air head. As an exception to this procedure, sample J23 #42b was collected by preserving the water siphoned from sample J23 #42 after 3 days of settling time.

For every sample, the in situ DO was measured using a Hach LDO sensor connected to a Hach HQ40d multi-parameter meter, and the in situ temperature-corrected EC was measured using an Orion conductivity probe connected to an Orion 4-star meter. The pH was measured on filtered samples. Measurements were taken with an Orion Ag/AgCl combination pH electrode connected to an Orion 4-star meter, which was temperature corrected by placing the conductivity probe in a separate vial of the same temperature. The temperature was kept constant during this procedure by placing the samples in an ice bath. At the beginning of each sampling day, the pH probe was calibrated with pH 4, 7 and 10 Orion buffer solutions. Total alkalinity was measured using a 0.01 mol L⁻¹ HCL solution dripped from a Hach digital titrometer to an endpoint titration of pH 4.5 using a few drops of a methyl-orange indicator. Alkalinity was averaged over three titrations per sample, giving an average error of ±8 µmol L⁻¹ (∼1.4% for July proglacial samples).

Till and bedrock samples were collected from the locations shown in Figure 1. The bedrock samples were taken by hammering out blocks from the intact bedrock just above the ice surface. Proglacial till sample 1 was collected under the terminus of the glacier by access through an open channel, while sample 2 was collected from a 2 m thick column of subaerially exposed till ∼50 m downstream from the terminus.

Sample analysis

The dissolved metals were analyzed using a Horiba Jobin-Yvon Ultima II inductively coupled plasma optical emission spectrometer (ICP-OES), and the anions were analyzed through ion chromatography, using a Dionex ICS-3000 system. The analytical errors for all species are ±3%, with the error growing exponentially near the detection limits given in Table 1. Procedural blanks were run through the ICP during the middle and end of sample analysis, with all cations below detection limits.

Table 1. Average South Glacier water composition by sample type with one standard deviation. Cations shown in lower rows with anions and remaining parameters in upper rows. All detection limits and sample concentrations are reported in µmol L⁻¹. SMF denotes the sulphate mass fraction, Alk_T denotes the total alkalinity and CBE denotes the charge-balance error

The suspended sediment samples were wet-sieved through a 75 μm mesh with tap water, followed by a flush with distilled water and subsequent drying in an oven at 90°C until no further weight loss was observed. Samples were then ground down with ethanol to <10 μm using a vibratory McCrone micronizing mill. The till and bedrock samples were first crushed in a jaw crusher, followed by a swing mill, and were later crushed down to 10 μm. All mineral samples were analyzed by powdered X-ray diffraction (XRD), and the normalized weight percent mineralogies were determined from the X-ray diffractograms using the International Center for Diffraction Database PDF-4 and Search Match software by Bruker, with further refinement using the Reitveld program Topas 4.2. All samples were reported as being ‘quantitative’, with the exception of sample J23 #42b, which only had enough clay for a ‘semi-quantitative’ analysis. The mineralogy of the bedrock as given by XRD was confirmed through optical microscopy. A whole-rock chemical digest was also run on the bedrock samples, which allowed us to estimate the An content of the plagioclase by setting the product of the weight percent plagioclase and the estimated molar fraction of Na equal to the Na content given by the chemical digest. Total rock digest was performed with lithium metaborate/lithium tetraborate fusion followed by 5% nitric acid digestions and analysis with ICP-OES.

Suspended sediment grain size was measured on separate samples using a Malvern Mastersizer 2000 particle size analyzer. Samples were wet-sieved through a 90 μm mesh, soaked in a 0.05% Na-metaphosphate dispersant solution for >24 hours, then placed in a sonic bath prior to injection. The grain-size distribution is likely finer than measured given the presence of inorganic colloids (Chanudet and Fillela, 2006) that did not completely deflocculate, and the presence of platy particles that align perpendicular to the Mastersizer laser.

Proglacial water discharge

Dilution gauging was performed according to the schedule in Figure 2. For each dilution, 200–600 g of table salt was premixed with river water and dumped in as slug injections. The conductivity of the pulse was recorded using a non-commercial half-bridge conductivity sensor (Stone and others, 1993) connected to a Campbell Scientific CR1000 data logger set to record at 1 s intervals. Experimental results of Day (1977) suggest that a mixing length of 25 times the average stream width is appropriate for high-gradient turbulent streams, so we used a mixing length of 50 m given the ∼2m stream width. The relationship between the salt concentration and the resulting conductivity was determined for each injection by measuring the increase in conductivity associated with adding known salt concentrations to a known volume of stream water. Conductivity required temperature correction, as the stream temperature at the terminus varied by up to 1°C on daily cycles. Temperatures were recorded with an Onset Hobo Tidbit V2 temperature logger.

Water stage was measured at 10 s intervals using a Campbell Scientific SR50 acoustic depth sounder, and stage data were filtered using a fast Fourier transform with a Butterworth filter centered at 2.9 × 10⁻⁴ Hz. Rating curves were generated between the instantaneous discharge and water level for each day of dilution gauging. High-runoff and mass-wasting events caused significant changes in stream channel geometry over the month of July, producing significant variation in the rating curves. The error function (erf) was used to create a smooth transition in slope and intercept between sequential rating curves; the onset and duration of the transition were based on field observations and the requirement that discharge remain positive (see Fig. 2). An estimated uncertainty of ±0.65 m³ s⁻¹ was applied to the whole record, with an additional uncertainty ranging linearly from ±0–50% over the course of each transition period (Fig. 2).

Geochemical Modelling Methods

Glacier water compositions are often modelled through an inverse mass-balance approach by stepwise calculations (Sharp and others, 1995a; Anderson and others, 2000; Hodson and others, 2000; Yde and others, 2005) and simultaneous solutions (Mitchell and Brown, 2008), with the goal of partitioning all ions in solution into their appropriate mineral phases and atmospheric inputs (e.g. Garrels and Mackenzie, 1967). Mass-balance calculations are also performed by comparing the ratios of ions in solution with the relative stoichiometries of plausible reactions (e.g. Raiswell, 1984; Tranter and others, 2002b). This approach puts an emphasis on sulphate and bicarbonate concentrations as a means of quantifying the proton availability for coupled dissolution reactions. We start with this approach by performing a series of inverse mass-balance calculations to highlight the silica deficiency in the South Glacier waters.

An alternative approach is to forward model the chemical composition with constraints on the mineral reaction rates, surface areas and stoichiometries. Under the assumptions of a chemical steady state and a dilute input source, Anderson (2005) implemented a forward modelling approach by conceptualizing the flow-through reactor as

(1)

where Q _w is the discharge (m³ s⁻¹), and for mineral j, r_j is the weathering rate (mol m⁻² s⁻¹), β_ij is the molar stoichiometry for species i and S_j (m²) is the surface area. In Eqn (1) the output concentration of the ith species, c_i (mol m⁻³), is corrected for the input concentration, , of the supraglacial water composition using the ion to Cl⁻ ratios of supraglacial water. We briefly consider this approach to get a second estimate of the silica deficiency in South Glacier waters.

To avoid assuming a chemical steady state and minerals far from saturation, we focus on forward modelling the hydrochemical evolution along a reaction path, by incorporating the effect of the saturation state and secondary mineral production beyond kaolinite (e.g. Bethke, 2007). Although our implementation of reaction-path modelling is poorly constrained, we use it to illustrate the qualitative behavior of the system for a range of plausible parameters. We attempt to replicate the low observed silica concentration through ordinary geochemical reactions influenced by physical processes in the subglacial environment (e.g. long residence times, basal freeze-on and mixing), and use reaction-path modelling to argue that it is not necessary to invoke mineral leaching to explain low silica ratios in glacial meltwaters.

Theoretical basis for reaction-path modelling

When considering the rate mechanism and the effect of saturation, while neglecting any transport, the overall reaction rate for a given mineral can be described by transition state theory as (Lasaga, 1998)

(2)

where n is the number of moles of mineral j dissolved, k_l is a rate constant that depends on the reaction mechanism, l, is the activity of hydrogen and is raised to the power m, as defined by the reaction mechanism l, Q is the activity product and K is the equilibrium constant for the reaction at temperature T (K). The reaction mechanism, l, is either acidic, neutral or basic. The quantity log(Q/K) is known as the saturation index (SI), and the equilibrium constant is set by the law of mass action as governed by thermodynamics. We use the Lawrence Livermore National Laboratory database ‘thermo.dat’ for the thermodynamic data of the minerals (Delany and Lundeen, 1990). Since thermodynamic data are not given for the plagioclase solution series, the K value for the plagioclase composition is calculated using the Margules parameters to account for the excess free energy available from solid solution (Holland and Powell, 2003). The intrinsic reaction rate, k, is a function of temperature, and defined by the Arrhenius relationship as

(3)

where k ₀ is the rate constant at the reference temperature (T ₀ = 298.15 K), E is the activation energy and R is the gas constant. Values for k ₀, m and E come from numerous experimental data compiled by Palandri and Kharaka (2004). No experimental data are present in the database for biotite above a neutral pH, so the basic pH mechanism is not included in the biotite weathering rates.

Input parameters

We use the Geochemist’s Workbench™ (GWB) React module to model the reaction path and plot the results on a series of log activity diagrams with the observed water chemistry. The GWB software speciates the sampled and modelled water compositions, and allows us to define the mineral reaction rate parameters, initial water composition, water volume, specific surface area for each mineral, mass of minerals to be reacted, reaction time and temperature. Most minerals are set to dissolve at rates governed by their kinetics, but the minerals gibbsite and kaolinite are assumed to be in equilibrium with the water at each time step, and are gradually added or removed over the course of the reaction. As a result, these minerals cannot become supersaturated, and their precipitation rates are governed by the rates of the dissolving minerals.

In the subglacial environment, the concentrations of gases in solution are controlled by access to the atmosphere, the composition and pressure of air bubbles in the basal ice (Hallet, 1976; Killawee and others, 1998), water pressure, microbial activity in the presence of organic matter (Sharp and others, 1999) and the overall mineralogy. As a result, the constraints on gas composition are poorly defined, so we often use the gas composition as a model parameter by defining equilibrium with an initial or final gas composition. The mineral surface areas, amount of mixing and simulation time are also used to constrain the reaction path.

Conceptual model

One of the largest uncertainties in modelling subglacial geochemistry is the amount of mixing between the channelized and distributed drainage systems. The distributed drainage system can have a wide range of water and sediment residence times and water/rock ratios, depending on whether water flows through unconsolidated sediment (e.g. Boulton and Jones, 1979), an interconnected system of linked cavities (e.g. Lliboutry, 1968; Kamb, 1987), or some form of water sheet or film (e.g. Weertman, 1957; Creyts and Schoof, 2009). Despite the complexity of the drainage system, models sometimes assume binary mixing by calculating changes in electrical conductivity of the proglacial stream, where the highest-conductivity waters reflect the slow/distributed endmember, and the low-conductivity waters result from mixing with the fast/channelized system (e.g. Collins, 1979; Raymond and others, 1995). However, mixing may not be a simple binary process between channelized and distributed waters (Tranter and others, 2002b), and conductivity variations in proglacial waters do not reflect conservative mixing, due to post-mixing reactions (Tranter and others, 1993; Sharp and others, 1995b; Brown and others, 1996). These complexities have led to the use of chemical markers, such as strontium isotopes, to calculate mixing (Hindshaw and others, 2011). Based on the above, and observations of hydrochemistry at South Glacier, we implement a three-stage conceptual model of the reaction path:

• Stage 1: reaction of dilute water with sediment.

• Stage 2: precipitation reactions driven by saturation, as a result of long water residence times and/or basal freeze-on.

• Stage 3: mixing of dilute supraglacial water with the subglacial water and sediment.

Reaction-path modelling is the focus of stages 1 and 3, while stage 2 is mostly guided by observations of the mineral saturation state and the presence of mineral precipitates within the system.

Observations

Proglacial and supraglacial water chemistry

The South Glacier proglacial stream displayed a typical glacier meltwater composition (Table 1). Although the basin is presumed to be composed solely of granodiorite, Ca²⁺ concentrations dominated the water composition, while Si concentrations were low. Over the July 2013 sampling period, the conductivity of the proglacial stream ranged between 50 and 142 μS cm⁻¹. The mean molar concentration of the proglacial water is used to guide stage 3 of the reaction-path modelling, and follows the order of

where subscript T denotes total, and the subscript number in parentheses indicates the average enrichment relative to Na⁺. The average charge-balance error of the proglacial water was 3.5%. The concentrations of major cations and anions were significantly lower in the supraglacial stream sample, with the exception of Cl⁻. The sulphate mass fraction, defined by Tranter and others (2002b) as SMF = [SO₄ ²⁻]/([SO₄ ²⁻] + [HCO₃ ⁻]), where concentrations are in equivalents, was ≪0.5 for all samples, indicating that sulphide oxidation was not the prevailing acid source within the water types we observed.

The concentration of all species decreased linearly with discharge, again with the exception of Cl⁻. For all species, the flux increased with discharge. There was an increase in pH from 6.8 to 8.2 over the course of the July sampling period, but this was the only parameter that varied systematically with time. At the terminus elevation of 1970 m, DO varied from 11.6 to 12.2 mg L⁻¹, and the stream temperature varied from 0 to 1 ± 0.2°C. Unlike the findings of Brown and others (1994), we were not able to correlate the DO content with the sulphate content or the discharge, and DO concentrations at the terminus were measured as being supersaturated. The average SSC from 11 to 27 June was 2.3 ± 1.0 g L⁻¹ (n = 31).

Water draining from the eastern tributary glacier (see Fig. 1) flows into an ice-marginal channel, and was very similar in composition to the water from the main proglacial stream. Other perennial streams draining the slopes of the basin had EC values of ∼20 μS cm⁻¹. The discharge was not measured in any of these streams, but appeared to be low relative to the proglacial stream. We assume that these streams have a negligible influence on the proglacial water composition.

Borehole water composition

The borehole water composition was substantially different from that of the proglacial water. Of the six samples collected, all had similar ionic ratios but dissimilar concentrations, as indicated by a range in conductivity of 1.5–106.7 μS cm⁻¹. Water for the borehole drill came from the supraglacial stream where the supraglacial sample was collected. The large range in borehole water composition likely results from variable amounts of mixing with dilute water, either from natural processes in the subglacial environment (discussed below), or from mixing with the drill water. We take caution in interpreting borehole samples because the boreholes were not instrumented with conductivity or DO sensors, as suggested by Gordon and others (1998, 2001). In addition, we do not know the chemistry of the basal ice, and although we speculate that basal freeze-on is occurring, we have no way to estimate the thickness of a freeze-on layer.

We justify using the borehole samples as indicators of the subglacial hydrochemical environment with the following arguments. Borehole camera footage at the base of a sampled borehole showed that the hole intersected the glacier bed, revealing fine-grained sediment and boulder-sized clasts. The water collected from all other boreholes had a similar composition, suggesting an environmental similarity between the boreholes. When focusing on the

and Ca²⁺ concentrations, the borehole water collected at South Glacier was similar to the ‘type C’ borehole water collected at Haut Glacier d’Arolla, which is thought to be basally connected (Tranter and others, 2002b). Type C water is characterized by proportionally high Si, Na⁺ and K⁺, and low

, Ca²⁺, Mg²⁺ and

. The absolute concentrations of Na⁺, K⁺, Mg²⁺, Si and

were slightly different between South Glacier borehole waters and Haut Glacier d’Arolla type C waters, which we attribute to differences in bedrock mineralogy between the two sites. The average molar composition of South Glacier borehole water follows the order

Unlike the proglacial input, we use individual sample compositions to guide stages 1 and 3 of the modelling, and not the average composition. The South Glacier borehole water samples were basic, with an average pH of 8.3 and a maximum of 9.29. The DO values of the only measured samples were 6.3 and 10.0 mg L⁻¹, which yields saturation values of 65% and 97%, respectively, for water at 0°C and atmospheric pressure at the ice surface elevation of drilling. These values are lower than the minimum proglacial and supraglacial concentrations. The DO was much lower for the sample that had a high EC, which likely reflects a lesser degree of mixing between the oxygen-depleted borehole waters and the oxygen-saturated supraglacial input. Given the low sulphate concentrations, the SMF was, on average, 0.05 in the borehole waters, which is much lower than the 0.25 SMF of the proglacial waters collected in July. Borehole waters were also much higher in Ba and B than the average proglacial water composition (Table 1). Lastly, while the grain size in the suspended load from the proglacial stream had a mean of 21.11 ± 5.65 μm (n = 30), the borehole grain size was on average finer, at ∼16.55 ± 10.75 μm (n = 27). The difference in grain-size distribution between borehole and proglacial waters likely reflects a grain-size fractionation from the sampling procedure.

Suspended sediment and bedrock mineralogy

The mineralogy of the suspended sediment in the proglacial and borehole waters was substantially different from the bedrock mineralogy (Table 2). Calcite, illite/muscovite and laumontite were below XRD detection limits in the bedrock samples, yet these minerals were present in all suspended sediments in both borehole and proglacial waters at low concentrations. In the proglacial waters, laumontite was in the range 4.9–6.8 wt%, illite/muscovite in the range 6–10 wt% and calcite in the range 0.2–1.1 wt%. Ankerite/ dolomite only appeared in sample PG M30 (0.9 wt%), smectite was only in samples PG M30 (6.8 wt%) and PG J23 #42b (7.0 wt.%), and gypsum was only in sample PG J23 #42b (1.3 wt.%). For the primary silicates, the weight percents of plagioclase and quartz were significantly lower in the suspended sediments (∼38 wt.%) than in the bedrock (∼50 wt.%), while the opposite was true for K-feldspar, actinolite and clinochlore. Although kaolinite and gibbsite were not reported in the XRD analyses, they may have been present as either coatings on grain boundaries, in an amorphous phase within the solution, or below detection limits.

Table 2. Estimated mineral compositions (%) with one standard deviation for bedrock, till and suspended sediment (SS) samples collected in July (J) and May (M) in the proglacial (PG) stream and from boreholes (BH). All samples are reported as quantitative, with the exception of the semi-quantitative sample J23 #42b, which was collected at the same time as sample SS J23 #42, but represents the mineral composition of the suspended sediments after 3 days of settling time

The chemical composition of the suspended sediment was inferred from the product of the estimated mineral stoichiometries and the calculated weight percent mineralogy. Our estimates indicate that the suspended sediments were enriched in Al, Mg and K, and depleted in Na and Si in relation to the bedrock. Similarly, Hasholt and Hagedorn (2000) found suspended sediments enriched in Al in relation to riverbed sediments, due to the mechanical separation of aluminosilicate minerals, while Hindshaw and others (2014) sampled river sediments that were also depleted in Na⁺ relative to the estimated bedrock composition. Both studies were conducted in Greenland. Mineralogy of the South Glacier till indicates a depletion in Mg, resulting from less clinochlore and the absence of micas within the till. Within the range of observed proglacial discharge, the suspended sediment concentration and flux both increased linearly with the discharge, with no apparent change in the mineralogy.

Modelling Results

We first quantify the silica deficiency of South Glacier meltwaters by highlighting results from mass-balance and flow-through reactor calculations. We then consider the chemical evolution of the meltwater through reaction-path modelling guided by an analysis of the mineral saturation and mineralogical data.

Chemical mass balance and the flow-through reactor

To approximate the silica deficiency in the proglacial water, we apply a series of stepwise inverse mass-balance calculations to individual samples (not the average compositions shown in Table 1) by first partitioning all of the Na⁺ in solution to andesine, based on the incongruent dissolution of andesine to kaolinite as

(4)

From Eqn (4) we see that the dissolution of andesine should add Na⁺ and Si to the water at a ratio of 1 to 2. However, the average Na : Si is 1.66 for the South Glacier proglacial water and 1.30 for borehole water, creating Si deficiencies by factors of 3.3 and 2.6, respectively. After this initial calculation, there is no additional Si to partition into the remaining silicate minerals. Alternatively, we could assume congruent dissolution of andesine as

(5)

which would lead to a Na : Si in solution of 0.25, but this would double the Si deficiencies of proglacial and borehole waters.

Results from the flow-through reactor model indicate that the modelled Si flux is at least 4 and up to 35 times the observed flux when considering the incongruent dissolution of the feldspars to kaolinite, and the congruent dissolution of actinolite, biotite and calcite. In brief, we compute the theoretical silica flux by calculating the surface area of each mineral that must dissolve to produce the observed Ca²⁺, Mg²⁺, K⁺ and Na⁺ fluxes. Because the system is underdetermined, surface area ratios are further constrained by thin section analysis and modal weight percent mineralogy given by XRD from bedrock samples. We obtain a best-fit estimate of the Si flux by tuning the aspect ratio and stoichiometry for biotite. The variation in excess Si values from 4 to 35 arises from testing a range of mineral weathering rates that are dependent on a pH from 5 to 10 (Palandri and Kharaka, 2004). The modelled silica excess further increases when considering a stoichiometric release of ions from all minerals, with a minimum excess ∼7 and maximum 59 times the observed flux. The above values are calculated under the assumption of steady-state weathering. We now move to stage 1 of the reaction-path modelling to explore the chemical evolution and silica deficiency through transient weathering simulations.

Stage 1: reaction of dilute water with sediment

Stage 1 begins by simulating the reaction of sediment in water at equilibrium with atmospheric CO₂. We use the supraglacial water composition, and a mineral content based on the average composition of the bedrock mineralogy given by XRD. The surface areas available for reaction are calculated for each mineral as a function of the weight percent values given by XRD, mineral densities, aspect ratios (height : width) and average grain sizes. The product of these variables allows us to compute a mineral surface exposure rate per liter of water (cm² h⁻¹ L⁻¹) (Table 3). Based on thin-section analysis from bedrock samples at South Glacier, and observed subglacial sediment shapes of primary minerals within the literature (e.g. Pandey and others, 2002), we assume that all minerals are rectangular, and assume an aspect ratio of 50 for biotite, 0.5 for actinolite and 1 for all other minerals. We use the grain size as a fitting parameter to produce a realistic reaction path.

Table 3. Model parameters for all simulations. Units are given in cm² h⁻¹ L⁻¹ to indicate the sediment surface area with which 1 L of water comes in contact per hour. Gibbsite and kaolinite are allowed to precipitate and dissolve in all simulations. We use phlogopite for the composition of biotite based on the availability of experimental data for phlogopite within the literature

The stage 1 simulations are designed to mimic the initial water-to-rock contact, and we are therefore interested in the results of the first few hours of reaction time. Experimental work relevant to glacier hydrochemistry has been carried out on these timescales by Brown and others (1996), who mixed distilled water with sediment at a concentration of 4 gL⁻¹ and a mean grain size of 65 µm. Using the conditions described by Brown and others (1996), and our initial mineral surface areas as calculated from the XRD results, we find that it is not possible to simulate the solute increase observed by Brown and others (1996) over the same 3 hour reaction period. This indicates that the simulated surface areas and/or dissolution rate constants are too low.

To constrain the total surface area available for reaction, we react enough calcite and sulphide to match the final [Ca²⁺] and of Brown and others’ (1996) experimental data. We then increase the sediment concentration so that the Mg²⁺, K⁺ and Na⁺ concentrations are roughly comparable to those observed over a 3 hour period. The experiment of Brown and others (1996) was open to the atmosphere, yet they found that the p(CO₂) still dropped to 10⁻⁵ atm (1 Pa) after a 3 hour sampling period, suggesting that CO₂ was being consumed at a much faster rate than it was able to diffuse into solution (Raiswell and Thomas, 1984). For this reason, we show two model runs: one run that has no endpoint constraint on CO₂ (run 1, Fig. 3a), and a second where we force the p(CO₂) to drop linearly to an endpoint of 10⁻⁵ atm (1 Pa) (run 2, Fig. 3a). Because the availability of oxygen is poorly constrained we do not model the pyrite oxidation rate using kinetics. Instead, we select an endpoint , and linearly add H₂SO₄ over the reaction path to simulate the same stoichiometric release of H⁺ and as pyrite oxidation.

Fig. 3. Reaction rate model Stage 1. (a) Reaction path with the supraglacial water composition as the initial condition. a is activity. (b) Change in species concentration versus normalized reaction time for run 2. See Table 3 for details.

The increase in pH is comparable between our simulation and the experiment of Brown and others (1996); however, the experimental results show a much more rapid initial rise in pH, which is attributed to the high initial reaction rates in the presence of freshly comminuted mineral surfaces (Petrovich, 1981). Although Brown and others (1996) do not show any Si data for comparison, we note that the silica concentration in the simulation is much higher than the concentrations of K⁺, Mg²⁺ and Na⁺. Figure 3a shows that the reaction path for run 2 intersects the proglacial water samples, but Figure 3b, which shows the evolution of the species concentration with time for run 2, indicates that the Si concentration is substantially higher than K⁺, Mg²⁺ and Na⁺. Again, this is not the case for a typical meltwater composition. This model produces a similar Si excess to the inverse mass-balance or flow-through reactor approach discussed in the previous section.

Run 1 in Figure 3 shows that gibbsite will react to form muscovite, while in run 2 gibbsite will react to form kaolinite, which later weathers to muscovite in the presence of K⁺. As the reaction path moves up the gibbsite/kaolinite boundary, the conversion of gibbsite to kaolinite can be coupled to the K-feldspar/kaolinite reaction, and as a result no excess silica is produced (e.g. Faure, 1998). This decreases silica, but not enough to produce the low Si ratios observed in solution, because too much silica is released before the path intersects the gibbsite/kaolinite boundary. The paths in Figure 3 are just two of many possible reaction paths. For example, run 1 can be reconstructed in a similar way, by adding more calcite, or reducing the amount of sulphuric acid. Again, these alternative paths do not create low-Si waters. In these simulations we use supraglacial water initially in equilibrium with the atmosphere, but one could also consider a basal meltwater input whereby the initial gas composition and gibbsite abundance are substantially different.

Stage 2: precipitation reactions and mineral saturation

If stage 1 waters become incorporated (or remain) in the slow drainage system, they are subject to chemical alteration through an increased water-to-rock contact time, increased solute concentration through basal freeze-on, or changes in the partial pressure of gases that result from changes in water pressure. When these processes occur faster than mineral precipitation can remove ions from solution, the water can become supersaturated with respect to various minerals. Figure 4 shows that some borehole and proglacial waters are at or above saturation with respect to calcite, laumontite, illite, smectite (beidellite-Mg) and clinochlore at elevated pH and cation activity levels. We focus on these minerals because they are found in the suspended sediment samples. With the exception of the samples collected in May, almost all samples are supersaturated with respect to muscovite (Fig. 4), gibbsite and kaolinite, among other iron oxides and oxyhydroxides (not shown). Given that only six borehole samples were collected, and all from a small area of the glacier, there are likely other subglacial waters with higher pH values and thus higher mineral saturation states. The borehole waters that were collected may be even higher in concentration, as refreezing (discussed below) concentrates solutes at the ice/ water interface in calm environments (Killawee and others, 1998). Waters at this interface cannot be measured because our samples are taken from a bulk sample spanning 30 cm above the bed. Similarly, solutes may be concentrated within etch pits at the grain boundaries of basal sediments (e.g. White and Brantley, 2003), but again, the composition of this thin film would not be reflected in the bulk water composition.

Fig. 4. Ratio of activity product, Q, to the equilibrium constant, K, versus ion activity or pH for various secondary minerals. Saturation is achieved where Q/K = 1. PG indicates proglacial.

Approaching saturation

To investigate the effect of long residence times, we ran the stage 1, run 1 simulation for an extended period of time, but at the end of the stage 1, run 1 time period, we decreased the input rate of sulphuric acid to reflect the low sulphate content of the subglacial waters. To investigate the effect of basal freeze-on, we linearly removed water from the system at the endpoint of stage 1, run 1. These two processes lead to similar results in terms of the secondary minerals that reach or approach saturation (Fig. 5). Calcite and clinochlore are among the first minerals to reach saturation, followed by laumontite. Here we use calcic beidellite (a proxy for calcic smectite), which slowly approaches saturation along with illite and dolomite. By the end of the simulation, K-feldspar has Q/K = 0.1 while andesine only reaches Q/K = 10⁻⁴ (not shown). Since the reaction rate for calcite rapidly approaches zero while the reaction rate for andesine is relatively unaltered, the simulated waters experience an increased Na⁺ to Ca²⁺ ratio, as described elsewhere in the literature (e.g. Anderson, 2005; Wadham and others, 2010).

Fig. 5. Secondary mineral saturation versus time obtained by running stage 1, run 1 for an extended period of time with a decrease in sulphate input at the end of the stage 1 simulation. Similar results are obtained by linearly removing water at the end of the simulation to represent basal freeze-on. Beidellite-Ca is a proxy for the calcic smectite, and clinochlore is represented by the 14 Angstrom (14A) variety.

In these simulations, the water becomes greatly enriched in solutes, yet the path does not lead to the observed composition of the borehole waters. Instead, gibbsite is consumed and the reaction trace moves right across the kaolinite/muscovite boundary (not shown). One way to move up from this boundary is to add enough gibbsite to the initial solution that the coupled gibbsite/kaolinite/muscovite reactions conserve Si. The impact of adding gibbsite is considered in the discussion.

The influence of secondary silicate precipitation on water composition

Although the modelling results are qualitatively useful, there are many limitations to what the stage 2 modelling can accomplish. Due to the high degree of substitution within the clay minerals, their thermodynamic properties are difficult to establish (Faure, 1998), and the XRD data do not provide an exact chemical composition for clays with variable composition. The saturation state is therefore estimated from an idealized endmember composition. For example, illite and muscovite may have an undetermined amount of K⁺ substitution between the two endmembers. A further limitation in using GWB to model secondary mineral formation is in constraining the precipitation rates. The reaction path is also dependent on initial conditions, such as mineral surface area and reaction rates, both of which are poorly constrained. For example, the reaction path can move up from the kaolinite/muscovite boundary if kaolinite is exhausted, but the amount of kaolinite depends on the initial conditions. To overcome these limitations, we shift our focus to consider the influence of mineral precipitation reactions on the water chemistry, by investigating the stoichiometry of plausible precipitation reactions.

The secondary minerals that are stable in South Glacier waters (Fig. 4) and found within the suspended sediments (Table 2) can be precipitated through the coupled dissolution of the primary minerals, gibbsite and kaolinite, along with other minerals as in the equations below:

(6)

Clinochlore occurs in much greater weight percent quantities in the proglacial and borehole water sediments than in the bedrock, and is supersaturated in some borehole waters. The precipitation of clinochlore is a good means of consuming Mg²⁺ and Si, while increasing K⁺. Furthermore, this reaction produces a substantial quantity of acid, which may help drive mineral dissolution as CO₂ and O₂ become depleted. Other possible examples of precipitation reactions in the subglacial environment of South Glacier are

(7)

and

(8)

Equations (7) and (8) also supply additional protons while consuming large amounts of silica. Note however that Eqns (6–8) are written as idealized compositions, and the silica consumption will vary greatly with substitution. The production of protons during these precipitation reactions may act to put an upper limit on the pH of subglacial waters.

Alternatively, it may be possible for the reaction path to move up from the kaolinite/muscovite boundary and into the muscovite stability field through the preferential release of K⁺ from interlayered mica sheets (e.g. Blum, 1997; Hasholt and Hagedorn, 2000; Anderson, 2005). The interlayer depletion of K⁺ is often accompanied by the addition of Mg²⁺ as biotite weathers to chlorite (Eggleton, 1986), and this could also help explain the depletion of Mg²⁺ in the borehole water. The importance of biotite dissolution as a source of K⁺ has been shown using Sr isotopes (Anderson and others, 1997; Sharp and others, 2002; Hagedorn and Hasholt, 2004; Hindshaw and others, 2011).

Basal freeze-on and mixing inferred from chloride in borehole and supraglacial samples

The motivation to explore basal freeze-on as a means of solute concentration arises from the elevated Cl⁻ concentrations in subglacial waters (see Fig. 6), and the observation of debris-rich ice near the terminus of the glacier. Chlorine is a conservative element that is released in negligible quantities when plutonic rocks are weathered (Holland, 1984), and for this reason, may be a good indicator of mixing between supraglacial and subglacial water (Sharp and others, 1995b). Cl⁻ can occur in elevated concentrations in snow and ice, as early-season elution preferentially causes a retention of Cl⁻ in the snowpack (Tranter and others, 1987; Hodgkins and others, 1997). Our analysis is limited by the absence of snow and ice chemistry data, and only one supraglacial sample with a charge-balance error of ∼50%. As a result, we cannot constrain input variations in Cl⁻. The elevated Cl⁻ could result from either basal freeze-on or from a concentrated input source. Either way, the key observation is that Cl⁻ is high in the subglacial water, thereby making it possible to estimate mixing fractions.

Fig. 6. Na⁺ versus Cl⁻ for proglacial and borehole water samples. Highly chlorinated samples likely result in enrichment from basal freeze-on. Symbol size is a function of discharge in the proglacial stream, which ranges from 0.8 to 4.2 m³ s⁻¹. The shaded squares used to fit line 2 represent the chemistry from 18 July 2014.

In analyzing the relationship between Na⁺ and Cl⁻ across all water types, we identify several linear correlations (Fig. 6). An apparent linear correlation for the borehole waters is shown by line 1 (Fig. 6) with a fit of y = x + 42.3. In an attempt to explain this correlation we explore two potential pathways in which the borehole waters could be subject to mixing and freezing. As a first possibility, the supraglacial water could undergo an initial mineral dissolution stage, releasing Na⁺ into solution before any freezing occurs (green line, Fig. 6). Where the green line intersects line 1 (point A, Fig. 6) freezing may be initiated, thereby increasing the Cl⁻ concentration (blue line, Fig. 6). In theory, a freezing line would have an intercept of zero, and for the line to pass through point A, the slope would be 24. Since we observe a slope of 1 for line 1, we would expect that 23 moles of Na⁺ are removed for every mole increase in Cl⁻ (thin black line, Fig. 6), indicating that Na⁺ is being removed by secondary mineral precipitation or cation exchange processes (e.g. Lorrain and Souchez, 1972; Tranter and others, 2002b). Plotting Mg²⁺ against Cl⁻ allows for the same analysis, and we find that Mg²⁺ is removed at a rate of 0.18 for every mole increase in Cl⁻. If smectite is the only secondary phase to remove Na⁺, and if we assume a Na smectite stoichiometry of Mg: Na as 2 : 0.3, then the amount of Na⁺ removed should be a maximum of 0.027 mole per mole Cl⁻, yet Na⁺ is removed at ∼850 times this rate. This could be a gross underestimate, since we do not account for Mg²⁺ that is released through cation exchange; however, the amount of Na⁺ required to compensate for the surplus through cation exchange would still be unreasonably large. The more plausible explanation for line 1 is that the borehole endmember composition at A’ (Fig. 6) results from freezing of supraglacial water (red line, Fig. 6). Line 1 can therefore be a mixing line between endmembers A and A’, where endmember A is again explained by solute enrichment prior to refreezing. Although we suggest that line 1 results from mixing, and the red line results from freezing, the slopes of the lines are poorly constrained, and the above-mentioned processes could feasibly be accompanied by some amount of mineral precipitation or cation exchange.

Line 2 is fitted to the 18 July proglacial samples, as shown by the filled squares in Figure 6. We omit sample PG J18 #28 from the regression because it shows anomalously high Si concentrations for proglacial water, and we exclude sample PG J18 #39, which is anomalously high in Ba and B. From the remaining samples collected on 18 July, we observe an average flux 484% higher in Cl⁻, 93% higher in Si, 81% higher in K⁺ and 71% higher in B than the average of the remaining samples for the July sampling period. The discharge is higher for this period, so the flux is expected to increase, yet there are only modest increases in Mg²⁺ at 37% and sulphate at 31%. This indicates that over the 18 July sampling period, the proglacial waters tend toward the borehole water composition. We have no further information to understand the cause of this event.

To estimate the amount of mixing between the subglacial and supraglacial water, and the amount of refreezing in the subglacial environment, we investigate the relationship between [Cl⁻] and proglacial discharge, and the ratio of subglacial to supraglacial Cl⁻. The increase in concentration, c, of a conservative element is inversely proportional to the amount of solution remaining, x, upon refreezing as 1/x = c/c _i, where c _i is the concentration of the solution prior to freezing. To increase the supraglacial [Cl⁻] from 1.8 μmol L⁻¹ to an endpoint of 70.1 μmol L⁻¹ (as indicated by the most chlorinated proglacial sample, #25, in Fig. 6) would require 97% of the water to freeze. To explore a more concentrated endmember, we could assume a final concentration that is double that of sample #25, but the amount of freezing only increases to 98.6%. Given the freezing point depression that arises from increased solute concentration, we verified that it is possible to freeze the South Glacier water samples to 98.6% at 273.15 K using the program Frezchem (Marion and Grant, 1994).

The concentration of a conservative element in a mixed solution can be represented as c _p = f _d c _d + (1 − f _d)c _c, where the subscripts p, c and d represent the proglacial, channelized and distributed components, f _d represents the mixing fraction, and we approximate c _c using the supraglacial [Cl⁻]. By rearranging this equation as f _d = (c _p − c _c)/(c _d − c _c), we can estimate the contribution of the distributed component to the proglacial stream as a function of the discharge (Fig. 7). Assuming that the endmember Cl⁻ composition, c _d, is represented by proglacial sample #25 (where f _d = 1), the mean mixing fraction from the delayed system over the July sampling period is calculated to be 0.14 ± 0.23. If we assume that c _d is twice as concentrated as sample #25, the resulting mean mixing fraction becomes 0.072 ± 0.11.

Fig. 7. Fractional contribution of distributed system, f _d, versus the proglacial water discharge, Q _w, estimated from [Cl⁻] in proglacial and supraglacial waters. Symbol types indicate sampling day. Mixing from distributed system increases as proglacial discharge increases, but the relationship appears to break down above ∼2.3 m³ s⁻¹.

At low discharge (⪝2 m³ s⁻¹), there appears to be a linear relationship between mixing fraction and discharge. Samples collected on the same day occur both above and below ∼2m³ s⁻¹, so diurnal changes in the f _d to Q relationship may reflect diurnal changes in the drainage system itself.

Stage 3: mixing of dilute supraglacial water with subglacial water and sediment

We now investigate post-mixing reactions by considering a simple mixing of supraglacial water with subglacial water, and consider additional cases where gibbsite and kaolinite are allowed to precipitate or dissolve both with and without the dissolution of primary silicates. To represent the mixing of supraglacial and subglacial water, we start with the simple mixing of the supraglacial sample with borehole samples BH 27_1 and BH 27_2. These point calculations are done at mixing ratios of 1 : 2, 2 : 1 and 10 : 1 (inverted triangular symbols in Fig. 8a). We find that the modelled compositions for the 2 : 1 and 10 : 1 mixtures have much lower K⁺/H⁺ and silica concentrations, yet from the previous section we calculated a mixing ratio of roughly 10 : 1. Either we are underestimating the component of distributed flow water and the subglacial K⁺ and Si concentrations, or there are additional reactions decreasing the acidity or increasing the ion concentrations during the post-mixing reactions.

Fig. 8. (a) Mixing of supraglacial and subglacial water (post-mixing) shown through simple mixing (inverted triangles), dilution of subglacial water with supraglacial water, kaolinite and gibbsite with no other mineral reactions (run 1), with mineral reactions and progressive dilution (run 2) and instantaneous mixing (run 3). a is activity. (b) Evolution of ion concentrations throughout the simulations of run 2 and run 3, which are shown in (a). See Table 3 for model parameters.

For a comparison to simple mixing, run 1 in Figure 8a (dashed curve) shows the reaction-path calculations for progressively diluting BH 27_2 with supraglacial water (see Table 3 for inputs). In this model, we assume an endpoint composition such that the water is in equilibrium with atmospheric CO₂. The reaction path is modified by the gibbsite/muscovite/kaolinite reactions, which distinguishes it from simple mixing. However, the simple mixing and reaction-path mixing are similar, in that the reaction path intersects the observed proglacial composition with a fraction from the distributed system that is much higher than estimated by the mixing calculations.

We now investigate the effects of adding sediment to the post-mixing reactions using two different approaches. In run 2 we simulate the instantaneous addition of sediment to the borehole water sample BH 27_1, then progressively dilute it with supraglacial water and add sediment over the duration of the run (Table 3). For run 3, we start with an initial water composition that reflects the simple mixing of supraglacial water with BH 27_1 at a ratio of 10 : 1, then react the water with sediment. In both runs 2 and 3, the final sediment concentration is 4 g L⁻¹. Reaction rates are initially higher in run 3 than run 2 because the solution is more dilute in run 3, and minerals are further from saturation. Run 2 is likely a more realistic simulation, but the end results for both runs are similar (Fig. 8b). In both runs, the silica concentration increases rapidly during the initial stages of the simulation, but as the reaction path reaches the gibbsite/ kaolinite boundary, Si is removed from the water creating a plateau in the silica concentration. Further dissolution causes K⁺, Mg²⁺ and Na⁺ to approach or eventually surpass the silica concentration by the end of the simulation. As a result, the model produces comparable ion ratios and concentrations to the observed proglacial stream chemistry.

Discussion

Glacier hydrochemistry in the framework of dissolution reactions

Without considering mineral precipitation, high Ca²⁺ to Si ratios are usually explained by rapid calcite dissolution rates (Raiswell, 1984; White and others, 1999a), while a higher than expected K⁺ to Si ratio is often attributed to the preferential release of interlayered K⁺ during the abrasion of biotite (Drever and Hurcomb, 1986; Hagedorn and Hasholt, 2004; Anderson, 2005). To explain high Na⁺ and Mg²⁺ to Si ratios, researchers often reference dissolution experiments of fresh mineral surfaces in acidic solutions, whereby a non-stoichiometric dissolution leaves behind a leached surface layer enriched in Al and Si (e.g. Wollast, 1967; Paces, 1973). For example, the carbonation of fresh feldspar surfaces is described as (e.g. Tranter and others, 1993, 2002b; Yde and others, 2005; Mitchell and Brown, 2008; Graly and others, 2014)

(9)

This proton-exchange reaction decreases with increasing pH, and at sufficiently high pH the reaction can be reversed (Wollast and Chou, 1985). In the glacial environment, borehole waters are often basic, proglacial waters are slightly basic (Tranter and others, 2002b) and supraglacial waters generally have a pH ∼6 (Wimpenny and others, 2010; Ryu and Jacobson, 2012); none of these conditions leads to enhanced leaching. If leaching were the dominant reaction, then it would follow that weathering is controlled predominantly by fresh mineral surfaces. These assumptions have yet to be verified, as the mineral surface age distribution under glaciers remains unknown, and the difference in weathering rates between fresh and weathered surfaces is not well established in the literature. It has also been hypothesized that the relative silica fluxes may be low in glacierized basins, due to a decrease in the weathering rates of primary silicates at low temperatures (White and others, 1999b; Anderson, 2005), but the observed high cation to Si ratios in proglacial waters cannot be explained by this temperature effect alone (Hindshaw and others, 2011).

Dissolution is not the only mechanism that is cited for controlling the fraction of silica in glacier meltwaters. It has been stated that the silica concentration is modified by the adsorption of cations onto mineral surfaces (Lorrain and Souchez, 1972; Tranter and others, 2002b), and the adsorption of silica onto the surface of clay particles (Siever and Woodford, 1973; Tranter and others, 2002b). The latter of the two processes is amplified by increasing the pH (Siever and Woodford, 1973). The measured silica content may also be altered through polymerization, which may prevent the silica from passing through the filter. At pH <9, over ∼95% of the aqueous silica exists as the H₄SiO₄ monomer, but higher pH values can lead to elevated and concentrations, thereby increasing the likelihood of polymerization in natural waters (Marshall, 1964; Sjöberg, 1996). Lastly, basal freeze-on can selectively segregate solutes into the basal ice (Killawee and others, 1998; Rea and others, 2004); however, the fraction of basal solutes is likely negligible until the final stages of freezing, when little solution remains (Hallet, 1976).

As an alternative means to explain the proglacial water quality, we have focused our modelling on weathering rates that depend on the saturation state, and mineral precipitation reactions. Our results do not contradict the idea that the water chemistry can be influenced by dissolution processes (e.g. release of interlayered K⁺ through biotite abrasion, mineral leaching or microbially mediated dissolution reactions). However, without having to rely solely on these dissolution processes, our results are consistent with the observations of a high pH, a low silica flux and the presence of secondary minerals in both proglacial and subglacial waters.

The influence of gibbsite on subglacial water composition

In the log (a _K+ /a _H+ ) − log a SiO₂(aq) phase space, the reaction path follows the gibbsite/kaolinite boundary, implying that any precipitated gibbsite is later being consumed by kaolinite, thus keeping silica low. Alternatively, if gibbsite is not present, we can assume that silica is still being consumed as kaolinite precipitates the excess Al³⁺ sourced from mineral dissolution. Although gibbsite is not found in the South Glacier waters, all samples are supersaturated with respect to gibbsite, with the exception of one sample collected in May. The stage 1 run 1 simulation produced 14 μmol L⁻¹ of gibbsite. This value is below the detection limit for XRD by a factor of ∼10, thereby providing a possible explanation for the absence of gibbsite in the XRD data. It is also possible that the gibbsite occurs in an amorphous phase, and is therefore not detectable by XRD. Alternatively, since gibbsite forms the octahedral layer for many of the two- and three-layer clays (Weaver and Pollard, 1973; Faure, 1998), it may be consumed by other secondary clays through further precipitation reactions.

Sulphate content and basal freeze-on

The observed South Glacier waters are low in sulphate, but it is possible that our boreholes did not intersect a high-sulphate water type, as found by Tranter and others (2002b). The low sulphate ratios likely result from a lack of oxygen at the bed. Since the sulphate concentration in the South Glacier water never exceeds the allowable limit imposed by the maximum oxygen saturation, we have no evidence for microbially mediated sulphide oxidation from the inorganic chemistry alone. Where regelation or basal freeze-on occur, premelted films around grain boundaries may act to concentrate sulphate and other ions to the extreme point of gypsum saturation. In this case, the enriched sulphate content in the proglacial stream may reflect gypsum dissolution, in addition to pyrite oxidation. Basal freeze-on may help to bring about mineral saturation in a shorter period of time and may be necessary to cause gypsum saturation, but it is not clear whether basal freeze-on is required for the saturation of secondary clays. The extent to which basal freeze-on controls the water chemistry could be better constrained by studies of temperate glaciers, where basal freeze-on is limited (e.g. to adverse bed slopes; Creyts and Clarke, 2010).

Total dissolution, CO₂ consumption and secondary mineral production

We compare the proglacial stream concentrations of K⁺ and Si in the secondary mineral phase to the K⁺ and Si concentrations in the aqueous phase (Fig. 9, with K⁺ shown by squares and Si by circles). For simplicity, we only consider the flux of K⁺ and Si for the secondary mineral phase within laumontite and illite/muscovite, which is estimated by multiplying the average laumontite and illite/muscovite weight percents over the July sampling period by the SSCs measured during dissolved ion sample collection. The analysis shows that the concentration of Si in the precipitated phase is ∼100 times higher than in the aqueous phase, while the K⁺ concentration is only ∼10 times higher than in the aqueous phase. If subglacial mineral precipitation is responsible for the observed secondary mineral fluxes, then we conclude that there is substantially more dissolution in the subglacial environment than would be predicted by measuring the dissolved ion flux alone. If this statement is true, then CO₂ consumption in glacierized basins is likely being underestimated by orders of magnitude when the chemical flux of the precipitated phase is not considered.

Fig. 9. Si and K⁺ in suspended laumontite and illite/muscovite versus aqueous Si and K⁺ in the proglacial stream. [Si] and [K⁺] in the solid phase are calculated by multiplying the mole percent of the mineral by the total suspended sediment concentration measured for each dissolved ion sample. The mole percent represents an average over five samples collected in July 2013, where less than a 5% variation in molar percent was observed for both illite/ muscovite and laumontite (Table 2). The apparent inverse relationship between elements in the aqueous and mineral phase results from the influence of water discharge.

When comparing the saturation state of the proglacial water to the mineral phases of the suspended sediment, it is clear that the mineralogy of the suspended load does not reflect the immediate chemistry of the water, and it would be unreasonable to couple the solid phase elemental flux to the dissolved ion flux. Given that water discharge and suspended sediment concentrations decrease in the winter months (Vivian and Zumstein, 1973; Benn and Evans, 1998), it is likely that the secondary mineral flux decreases substantially in winter, while solute fluxes remain relatively high. Longer records than ours would be required to compare the elemental flux from the solid and dissolved loads over an entire season.

On average, laumontite and illite/muscovite accounted for 15.2 ± 7.2% by weight (n = 9) of the suspended sediments in the proglacial stream over the July sampling period. In considering all suspended sediments potentially derived from secondary mineral precipitation (laumontite, illite/muscovite, calcite, smectite, chlorite, gypsum), the secondary mineral flux is calculated to be 33.9 ± 6.0% (n = 9) of the total suspended sediment flux. Normally, the suspended sediment flux is attributed entirely to physical erosion, but here we suggest that a significant portion of the flux is attributed to chemical erosion. However, partitioning erosion rates between physical and chemical may be difficult, given that chemical weathering in the subglacial environment is driven largely by high physical erosion rates (Anderson, 2005).

The mineralogy of the till indicates that secondary minerals (e.g. calcite and laumontite) can be stored within the till, while other secondary phases (e.g. illite) along with other primary minerals (e.g. biotite) are transported out of the system through mechanical separation. Therefore, till composition may not be a good indicator of authigenic minerals forming subglacially, nor a reliable representation of bedrock mineralogy, even in a geologically homogeneous basin. As subglacial calcite precipitates have been found to coat bedrock surfaces, and we find that calcite is at saturation for several of the borehole samples, we are also led to question the source of calcite in the suspended sediments. In addition to disseminated calcite within the bedrock, we infer that calcite could be sourced from mechanical abrasion of subglacial calcite deposits, calcite coatings on grain boundaries or precipitation from solution.

Research implications

As discussed above, this work provides motivation for reanalyzing the magnitude of chemical weathering, and thus the amount of CO₂ drawdown in glacierized basins. On a global scale, such a reanalysis may lead to a better understanding of the temporal role of subglacial chemical weathering on climate during glacial periods (e.g. Gibbs and Kump, 1994; Tranter and others, 2002a). In focusing on the mineralogy of sediment leaving the glacier system, this research may prove useful for understanding the genesis and source of clay in offshore sediments (e.g. Ehrmann and others, 1992) and other present and past proglacial environments where till geochemistry may be of importance (e.g. Anderson and others, 2000; McMartin and McClenaghan, 2001). Moreover, a better understanding of chemically driven sediment production rates in the subglacial environment may provide insight into the production and mechanics of till, which plays an important role in controlling basal water pressure and thus the dynamics of basal sliding (e.g. Iverson, 2010). Lastly, we suggest that reaction-path modelling may be of use in future glacier geochemistry studies, as we acquire better knowledge of dissolution and precipitation rates as a function of sediment age (e.g. White and Brantley, 2003; Anderson, 2005) and microbial activity (e.g. Montross and others, 2013; Kenward, 2014).

Conclusion

Hydrochemical and mineralogical sampling from a polythermal glacier in Yukon has allowed us to test the hypothesis that the hydrochemical signature of this basin is largely controlled by the combined influence of dissolution and secondary mineral precipitation reactions. Our approach to modelling was driven by observed water compositions depleted in silica and saturated with respect to secondary minerals that were present in suspended sediments within borehole and proglacial waters, but not within the bedrock. We found that mineral dissolution in the absence of abundant CO₂ and O₂ creates a basic environment that is conducive to secondary mineral precipitation.

We explored the silica deficiency through inverse mass-balance and forward flow-through reactor models, under the assumption of steady-state dissolution, but were not able to balance the water chemistry. The reaction-path modelling also highlighted silica deficiencies within the initial stages of water/rock reaction, and was able to simulate waters that were saturated with respect to the appropriate mineral phases. Where the modelling was limited by mineral precipitation rates, the stoichiometry of thermodynamically favorable precipitation reactions helped to explain the silica deficiency in solution, shed light on an additional proton source and provide insight into the low Mg²⁺ and relatively high K⁺ concentrations.

Simulations of long water/rock contact times were able to produce waters that reached secondary mineral saturation, but given the elevated Cl⁻ concentrations in supraglacial and subglacial waters, we speculate that basal freeze-on could have also played a role in concentrating solutes and driving secondary mineral saturation. We used chloride to estimate the contribution of distributed system waters to the proglacial stream, and thus to understand the extent to which the distributed system water composition controls the observed proglacial hydrochemistry. Reaction-path modelling indicated that the low-silica signature of subglacial water can be preserved, even into the proglacial environment, by post-mixing reactions which evolve the composition of subglacial water along the gibbsite/kaolinite boundary. Finally, the large measured flux of secondary minerals exiting the South Glacier system suggests a substantially higher net dissolution than would be predicted from the dissolved ion flux alone, and hence a much higher net consumption of CO₂.