This study examines current subglacial biogeochemical weathering models using PHREEQCi, a computer-based speciation mass-balance (SMB) model, parameterized using hydrochemical and mineralogical field data from Haut Glacier d'Arolla, Switzerland. The aim is to investigate the utility of SMB models for quantifying subglacial biogeochemical mechanisms, and the mathematical and thermodynamic robustness of current subglacial weathering models in a temporally variable, spatially heterogeneous hydrological and geochemical subglacial environment.
The chemical evolution of meltwaters between their source in the supraglacial environment and their output from the hydroglacial system as bulk meltwaters are modeled, and compared with in situ meltwaters collected from the base of boreholes drilled to the glacier bed. SMB modeling produced a broad range of weathering outcomes, but delivered no unique weathering scenario which could account for the observed changes in water chemistry between input and output waters over the ablation season. Organic carbon oxidation and sulfide oxidation by dissolved O2, coupled to carbonate dissolution and incongruent silicate dissolution, could account for seasonal changes in meltwater chemistry, supporting current subglacial biogeochemical weathering models. Atmospheric CO2 was not required under any weathering scenario, and organic carbon and atmospheric CO2 dissolution was only possible in one weathering scenario, at a mass ratio of 10:1, further suggesting that CO2-driven dissolution is unimportant in subglacial environments. This investigation indicates the utility of SMB models applied to subglacial hydrological systems for determining geochemical and biogeochemical processes and interpreting water quality.
The study of meltwater chemistry in subglacial environments concerns the use of hydrochemistry to determine rates of biogeochemical denudation and carbon cycling in order to assess the long term implications for global biogeochemical cycles (Sharp et al., 1995; Jones et al., 2002), to determine nutrient and metal export and cycling (Mitchell et al., 2001, 2006; Mitchell and Brown, 2007; Hodson et al., 2004; Brown, 2002), and to infer hydrological flow-paths (Brown et al., 1996a; Mitchell et al., 2006). Such studies require reaction processes and solute provenance to be accurately defined (Raiswell, 1984; Sharp et al., 1995; Tranter et al., 1993, 2002, 2005). They have often attempted to de-convolute complex bulk meltwater quality compositions by apportioning each solute to the dissolution of a particular rock type or gas in order to disaggregate solute into crustal, atmospheric, and snowpack sources, and to identify variations in solute sources and weathering mechanisms at seasonal timescales (Raiswell, 1984; Tranter et al., 1993; Fairchild et al., 1994; Sharp et al., 1995; Brown et al., 1996a; Anderson et al., 2000). However, more accurate definitions of solute provenance are necessary if meltwater composition is to be used unequivocally to disaggregate solute fluxes to particular reaction processes. This is particularly pertinent given the recent discovery of flowing liquid water under the Antarctic ice sheet (Wingham et al., 2006; Fricker et al., 2007) and the implications this has for understanding biogeochemical cycles in previously unexplored ecosystems.
Previous and Current Models of Subglacial Biogeochemical Weathering
Studies at Haut Glacier d'Arolla, a warm-based alpine glacier in Switzerland, have been central to our understanding of links between subglacial hydrology, hydrochemistry, and glacier dynamics, due to intensive interdisciplinary field studies over the last two decades (e.g. Tranter et al., 1993; Brown et al., 1996a; Sharp et al., 1995; Brown, 2002; Hubbard et al., 1995, 1998). Until recently, studies at Haut Glacier d'Arolla with the preceding studies of Raiswell (1984) gave rise to a generic chemical weathering model which has been used to explain and deconvolute bulk meltwater chemistries observed at a range of alpine and polar glaciers. This model, which we will call the T-Model (Tranter et al., 1993), is based upon sequentially applied mass balance calculations, representing a budget that describes fluxes of solutes into and out of a glacier, and assigns solutes to specific sources and sinks (Brown et al., 1996a; Sharp et al., 1995; Hodson, et al., 2000) (Table 1, reactions 1 to 5). In the T-Model, subglacial solute acquisition proceeds by acid hydrolysis reactions, fuelled by protons derived from either the dissolution of atmospheric CO2, or the oxidation of sulfide minerals by dissolved oxygen (Tranter et al., 1993).
However, the measurement of in situ subglacial meltwater chemistries derived from the base of boreholes drilled to the base of Haut Glacier d'Arolla have suggested problems with the T-Model which were not apparent from studies using bulk meltwater chemistry alone (Tranter et al., 1997, 2002). First, it was determined that CO2 was not a major driving force for subglacial chemical weathering reactions, but instead there was a progression of reactions from carbonate and silicate hydrolysis through to sulfide oxidation driven by dissolved O2. Once O2 had been depleted and anoxic conditions were achieved, carbonate and silicate weathering was hypothesized to be driven by electrons from the reduction of Fe3+. Central to these new hypotheses was the discovery of viable microbial populations at the beds of alpine and polar glaciers. Microbes are abundant at the water-rock-ice interface beneath Haut Glacier d'Arolla (>106 cells mL−1 in some subglacial waters) (Sharp et al., 1999), as well other alpine and polar glaciers including Tsanfleuron, Switzerland (Sharp et al., 1999), Fox and Franz-Josef Glaciers, New Zealand (Foght et al., 2004), John Evans Glacier, Nunavut, Canada (Skidmore et al., 2000), and Bench Glacier, Alaska, U.S.A. (Skidmore et al., 2005).
Second, organic carbon had been previously overlooked as a source of dissolved inorganic carbon (DIC) (Sharp et al., 1999; Skidmore et al., 2000, 2005). Indeed, subglacial debris contains organic carbon which is a potential source of energy for microbes, including kerogen, bacterial necromass, and overridden organic matter. This suggests that microbially mediated reactions, utilizing the oxidation of organic carbon by O2(aq) [as well as abiotic sulfide oxidation by O2(aq)] would drive the bed in places to anoxia. In these anoxic ‘hot spots’ sulfide oxidation by Fe3+, sulfate reduction (Wadham et al., 2004), and the reduction of organic carbon (methanogenesis) was proposed to occur (Tranter et al., 2002, 2005) (Table 1).
The configuration and dynamics of the subglacial drainage system will likely control such biogeochemical gradients and the spatial dynamics of reaction processes (Tranter et al., 2005). Alpine subglacial hydrological systems broadly comprise either a distributed or a channelized system (Fig. 1). The distributed system comprises linked cavities, permeable sediments, film flow, and broad canals, and it transports delayed-flow waters. This component is fed largely by snowmelt and is the dominant mode of subglacial drainage at the end of the winter, and at low diurnal discharges. Here, access to the atmospheric gases will be limited suggesting that microbes will initially oxidize organic matter and sulfides using dissolved O2, driving the bed towards and to anoxia. Then reactions may proceed anoxically, via Fe3+ reduction and methanogenesis (Tranter et al., 2002; 2005). Conversely, a hydraulically efficient dendritic channelized subglacial system develops during the melt season as surface ablation and subglacial discharges increase (Nienow et al., 1998), evolving to accommodate the increasing volume of meltwater generated during the summer months. These arterial channels expand up-glacier following the retreat of the seasonal supraglacial snowpack (Nienow et al., 1998), and drain waters which are comprised predominantly of icemelt. These channels form the dominant mode of meltwater drainage later in the ablation season, and at maximum diurnal discharges (Fig. 1). Therefore, in contrast to the distributed system, a plentiful supply of oxygenated water will promote microbially mediated organic carbon and sulfide oxidation by dissolved O2. However, it is likely that significant microbial activity will be limited to basal sediments in the distributed system, as well as the channel marginal zones (Tranter et al., 2005). Here the residence time of water is long, and hydraulic shear stress will be low, which will promote the stable attachment of microbes to the sediment or bedrock substrate and will allow prolonged access to both solid phase and dissolved nutrients.
Additionally, these current models only crudely consider the relative importance of silicate and carbonate weathering. Specifically, Ca2+ is the dominant cation in bulk meltwaters in a range of glacierized alpine and polar catchments, and despite being underlain by vastly varied compositions of silicate and carbonate lithology, Ca2+ is thought to derive dominantly from the dissolution of trace carbonates by chemically limited weathering regimes (Brown, 2002; Tranter et al., 1993; Drever and Zobrist, 1992; White et al., 2005). However, the potential for significant silicate-derived Ca2+ has rarely been given any quantitative consideration (cf. Anderson et al., 1997, 2000), despite the implications to solute provenance calculations as a whole. While much of these theories were formulated from studies at Haut Glacier d'Arolla, field studies from a range of warm-based and polythermal glaciers support these assertions (Anderson et al., 1997, 2000; Anderson, 2005; Wadham et al., 2004; Skidmore et al., 2005; Hodson et al., 2004).
Alternative Approaches for Determining Subglacial Geochemical and Biogeochemical Reactions
The development of freely available computer based geochemical reaction models provides a potentially useful tool for improving current models of subglacial biogeochemical weathering in glacierized environments. These models, like the previous T-Model, work on the mass balance approach, but they allow multiple alternative weathering scenarios under different hydrological conditions to be tested simultaneously, unlike previous mass balance calculations which require mass balance equations to be solved sequentially, and ad hoc, in a spreadsheet or in a longhand mass balance matrix (Finley and Drever, 1997). They also easily allow mineral-specific weathering processes to be modeled, the thermodynamic feasibility of the reactions to be assessed, and alternative estimates to be made of the consumption and/or production of important species not easily analyzed [e.g. CH2O(s), CO2(g)], under differing weathering scenarios.
This paper evaluates current subglacial biogeochemical weathering models using speciation mass-balance (SMB) hydrochemical models parameterized with hydrochemical, hydrological, and mineralogical data collected from Haut Glacier d'Arolla during the melt seasons of 1999 and 2000 (Mitchell et al., 2001, 2006; Mitchell and Brown, 2007). The aim is to investigate the utility of SMB models and to provide a better quantitative understanding of these most recent subglacial biogeochemical weathering models. Specifically, we (1) present seasonal estimates on the oxidation of organic carbon, (2) quantify if organic carbon can solely provide all non-CaCO3 DIC, (3) provide a more robust understanding of silicate weathering than has been previously considered in subglacial environments, and (4) test the thermodynamic feasibility that Fe3+(aq) can act as a terminal electron acceptor and oxidize sulfides under anoxic conditions. Haut Glacier d'Arolla provides an ideal case study with which to parameterize the study hydrologically and hydrochemically, and to compare and contrast with the existing studies of subglacial solute provenance and flow routing.
We used the SMB model PHREEQCi (Parkhurst and Appelo, 2000), which calculates the weathering transfers (dissolution/precipitation/de-gassing) of minerals and gases (phases) between two or more hydrologically connected sampling points. In simple terms, this is achieved by defining the input water and output water solution chemistry, and by determining the minerals and gases available to flowing water in the subglacial environment; the model then calculates a number of weathering scenarios, which are specific combinations of minerals and gases that can account for the observed changes in solution chemistry. The model outputs the mass of phase dissolved into solution per liter of water. This assumes the volume of water into the subglacial system from the supraglacial environment equals that exiting the subglacial environment. In practical terms, the composition of the input waters makes no difference to the output models, since the input waters are so dilute relative to the in situ and bulk meltwaters (Mitchell et al., 2001, 2006) that they can be taken to have no dissolved concentration.
The SMB model initially estimates the dissolved chemical form, or speciation, of the measured solutes in the input and output waters, checks that the solutions are charge balanced, and then using a suite of simultaneous mole-balance equations, calculates different weathering scenarios that can account for the changes in solution chemistry between the input and output waters. Such models do not include kinetic parameters. Therefore, while microbial enzymatic activity may increase some reaction rates above the abiotic rate, for example the oxidation of organic carbon (Tranter et al., 2005), the same reaction stoichiometry will apply, and thus the mass balance approach can be used.
This furthers simple mass balance approaches previously used (e.g. Tranter et al., 1993; Finley and Drever, 1997; Sharp et al., 1995; Anderson et al., 2000) by the inclusion of a more complete set of mole-balance equations and the addition of inequality constraints that allow for uncertainties in the analytical data. Mole-balance equations are included for (1) each element or, for a redox-active element, each valence state of the element, (2) alkalinity, (3) electrons, which allows redox processes to be modeled, (4) a charge-balance equation for each aqueous solution, and (5) an equation that relates uncertainty terms for pH, alkalinity, and total dissolved inorganic carbon for each solution. Furthermore, inequalities are used (6) to constrain the size of the uncertainty terms within specified limits (based upon analytical accuracy of the water chemistry data input into the model), and (7) to constrain the sign (+ve = dissolution; −ve = precipitation) of the mole transfer of reactants and products. The mole-balance equations, including the uncertainty terms and redox reactions, for elements and valence states are defined as:Parkhurst and Appelo, 2000).
The saturation index (SI) of the input and output waters with respect to various minerals and gases under equilibrium conditions is also calculated by the model (Equation 2);Plummer et al., 1983; Parkhurst and Appelo, 2000). When IAP = Ksp, then SI = 0, and the aqueous system is at thermodynamic equilibrium with the mineral or gas. When SI > 0, the water is oversaturated with the mineral or gas, suggesting the mineral or gas is likely to precipitate or de-gas to maintain chemical equilibrium. However, when SI < 0 the water is undersaturated with the mineral or gas, suggesting the mineral or gas will dissolve to maintain equilibrium in a given system (Plummer et al., 1983; Parkhurst and Appelo, 2000; Mitchell and Brown, 2007). Thus, mole-balance equations determine the reactions that can mathematically account for the observed changes in water chemistry based upon the speciated water chemistry and the stoichiometry of the proposed reactants and reactions, but this does not take into account the results of saturation calculations directly in the weathering scenarios generated. The model user has to consider saturation data to help constrain the model further (Plummer et al., 1983). Speciation and saturation-state calculations assume thermodynamic equilibrium, which may not be achieved in turbid glacial meltwaters for many hours after initial water-rock interaction (Brown et al., 1996b; Mitchell et al., 2001, 2006). However, few natural hydrological systems are in true thermodynamic equilibrium, but assuming chemical equilibrium provides a good quantitative basis for deconvoluting complex natural water chemistries, and is our best estimate of the likelihood of mineral dissolution or precipitation (Langmuir, 1997).
Haut Glacier d'Arolla is a temperate glacier ∼4.2 km long, with a maximum ice thickness of ∼180 m, and occupies ∼6.3 km2 of a ∼12 km2 catchment, overlying unconsolidated sediments which are between 0.05 and 0.26 m thick (Sharp et al., 1993). The bedrock consists largely of igneous and metamorphic rocks of the Arolla series of the Dente Blanche Nappe (Dal Piaz et al., 1977; Mazurek, 1986) (Fig. 2; Table 2). These major units contain trace amounts of geochemically reactive trace minerals such as calcite (up to 0.58%) and pyrite (up to 0.71%) (Tranter, et al., 1997). Brown et al. (1996b) recorded concentrations of carbonate up to 12% and pyrite up to 5% in the fine fraction (<5 mm diameter) of the moraine.
The major geological units and associated mineralogy at Haut Glacier d'Arolla. Mineralogy based on the composition of fines (<5 mm ø) from the major moraines at Haut Glacier d'Arolla.
The model was parameterized with the chemical composition of waters collected during the 1999 ablation season. Bulk meltwaters were collected twice daily near the glacier snout at 10:00 and 17:00 local time (approximating to minimum and maximum diurnal discharges, respectively) between 16 June (calendar day [CD] 167) and 20 August (CD 232) 1999. Supraglacial meltwaters and snow samples (pre-melt snowpack and fresh snowfall) were also collected. Samples were analyzed for major ions, minor and trace elements, and pH (see Mitchell et al., 2001, 2006 for further details). The pE of meltwaters was initially estimated as 7, representing typical oxygenated surface waters (Langmuir, 1997). This was used instead of redox potential measurements because they are often very inaccurate, reflecting interference of the redox electrode with suspended sediment, and the inability to filter out such sediments rapidly enough before redox potential determination. The model is also parameterized using all the main minerals in the Haut Glacier d'Arolla catchment (Fig. 2 and Table 2). The choice of phases is constrained by idealized mineral stoichiometries present in the MINTEC database used by PHREEQCi. The chosen phases that will be used in the following model development are listed in Table 3. Quartz, feldspars (anorthite; microcline, albite), micas (muscovite, biotite, chlorite), calcite, and pyrite can be defined in the model. Tremolite is chosen to represent amphiboles. Magnetite, epidote, and sphene are ignored in the modeling process due to their low abundance and resistance to chemical weathering (Deer et al., 1992). Clays can also be defined, including kaolinite and montmorillonite, which have been identified in the catchment. We also make organic carbon available as a simple carbohydrate (CH2O), which can be available in the form of bacterial necromass or young organic matter (Wadham et al., 2004; Tranter et al., 2005), as well as CO2(g) and dissolved O2, although this will depend upon the open or closed nature of hydrological environment, as we will demonstrate.
Phases (minerals and gases) used for speciation mass-balance modeling applied to the subglacial hydrological system of Haut Glacier d'Arolla.
First, we investigate the spectrum of feasible geochemical and biogeochemical weathering processes that can account for the observed changes in water quality from supraglacial to bulk meltwaters using the mean meltwater composition from the 1999 ablation season. Such models represent net reactions occurring throughout the subglacial hydrological system. Second, the number of weathering scenarios is reduced by constraining the weathering parameters of the model. This is achieved by allowing individual phases to only dissolve or precipitate. We also test weathering scenarios between supraglacial and in situ subglacial meltwaters. With this model parameterization, any changes in the viable weathering processes during the melt season are investigated by using seasonal variations in supraglacial and bulk meltwaters as the input and output waters (Mitchell et al., 2001, 2006).
Results and Discussion
Discharge and Meltwater Chemistry
Supraglacial waters are very dilute and exhibit very little seasonal variation, but become heavily enriched in solute as they travel though the subglacial hydrological system, as demonstrated by the more concentrated in situ and bulk meltwaters. The median composition of supraglacial meltwater, bulk meltwater, and in situ meltwater was, respectively, for the following ions, Ca2+ = 0.32, 6.6, 15 mg L−1; Mg2+ = 0.014, 0.4, 0.5 mg L−1; Na+ = 0.034, 0.37, 0.58 mg L−1; K+ = 0.036, 0.37, 0.35 mg L−1; SO42− = 0.075, 5.2, 9.7 mg L−1; HCO3− = 1.1, 17, 33 mg L−1; NO3− = 0.025, 0.4, 0.77 mg L−1; Fe = 0.083, 0.39, 0.56 mg L−1; and pH = 7.5, 8.1, 7.2. The bulk discharge magnitude and daily amplitude generally increased as the melt season progressed, reflecting variations in radiation receipts at the glacier surface and the up-glacier recession of the seasonal snowpack (Fig. 3). Four hydrological periods (P1 to P4) have been defined during the 1999 study period, based upon qualitative analysis of bulk discharge records (Fig. 3). P1 and P4 exhibit low discharges relative to P2 and P3. Calcium, Mg2+, K+, Na+, SO42−, and HCO3− exhibit seasonal variations, with declining concentrations during P1 as discharge magnitude and diurnal variability increased, reaching the lowest seasonal concentrations in P2 at the height of the ablation season. There is a general increase in concentrations over P3 and P4 as discharge falls towards the end of the ablation season.
Testing Viable Chemical Weathering Reactions in Subglacial Environments Using Phreeqci
We first investigated all viable reaction combinations allowing the unconstrained dissolution or precipitation of most of the main minerals available in the Haut Glacier d'Arolla catchment (Table 3) using the mean melt season concentrations of Ca2+, Mg2+, Na+, K+, SO42−, Fe, Al, and Si, and pH, for supraglacial and bulk meltwaters (Mitchell et al., 2006). This is required as a starting point to test the model, and represents net geochemical reactions occurring throughout the subglacial hydrological system on average over the melt season. Temporal variations in chemical weathering processes using the seasonal meltwater quality data will be considered subsequently in the Section COMPARISON OF MODELING-DEFINED CHANGES IN SUBGLACIAL WEATHERING PROCESSES OVER SEASONAL TIMESCALES. This produced 136 weathering scenarios, but unsurprisingly no individual weathering scenario that could account for the evolution of meltwater in the heterogeneous subglacial environment. The phase transfers of these weathering scenarios indicate dissolution of calcite, pyrite, tremolite, and O2(aq), and precipitation of ferrihydrite, anorthite, and microcline. While such reactions are feasible based upon reaction stoichiometry, mineral saturation index (SI) calculations indicate that primary silicate minerals present in the Haut Glacier d'Arolla catchment are undersaturated in solution over the range of measured pH and bulk meltwater solutions in the system (Fig. 4). Therefore, primary silicate mineral precipitation is unlikely. More likely is the formation of clay minerals from the incongruent weathering of micas, feldspars, and other aluminosilicate minerals, as has been previously proposed (Fairchild et al., 1994; Tranter et al., 1993; Anderson et al., 1997).
Therefore, kaolinite and montmorillonite were chosen as two representative clays for the model; these clays commonly form from the weathering of micas and feldspars and have been observed in catchment bedrock (Tables 2 and 3). When these clays were set as precipitating phases, 47 weathering scenarios were produced. Individual weathering scenarios (WS) produced by the parameters of this model indicate montmorillonite, or kaolinite and montmorillonite, precipitation is required by the model to remove the incidence of all primary silicate precipitation (Table 4). Kaolinite-only weathering scenarios cannot remove the incidence of primary silicate precipitation. The influence of clays in the mass balance calculations indicates how aluminosilicate dissolution and clay formation are covariant. This demonstrates the model is able to account for the changes in subglacial chemistry via incongruent silicate weathering.
Results of model testing between supraglacial and bulk meltwaters, using the mean composition of supraglacial and bulk meltwaters as input and output waters, respectively. Each line is a single weathering scenario (WS) that can account for the changes in solution chemistry between supraglacial and bulk waters. Weathering scenarios highlighted in light gray indicate the constrained model with oxygen-driven organic carbon oxidation and pyrite oxidation driven carbonate and incongruent silicate weathering. Weathering scenarios highlighted in mid gray are those which include calcite and anorthite weathering. Weathering scenarios highlighted in dark gray are those identified as being the most feasible weathering scenarios, based upon mineral dissolution rates. Units = mass of phase (mmol L−1) transferred between input supraglacial meltwater and bulk meltwaters, with positive values showing dissolution and negative values showing precipitation.
Dissolved oxygen is consumed in all of the weathering scenarios. Oxygen-inclusive weathering scenarios always include pyrite dissolution and ferrihydrite precipitation, indicating pyrite oxidation by dissolved O2. Pyrite is the only likely source of S since gypsum has not been recorded in the Haut Glacier d'Arolla catchment, unlike some other glacierized catchments (Skidmore et al., 2005), and thus the model cannot find any weathering scenarios unless pyrite and dissolved O2 are consumed (Table 4). This is in agreement with current models of solute acquisition in subglacial environments, although the effect of limiting O2 to simulate suboxic conditions will be subsequently tested (Table 1).
Organic carbon may be either consumed or completely absent from the weathering scenarios. CH2O and calcite account for the observed DIC measured in output waters (Table 4). Non-CH2O models always require the precipitation of a primary silicate mineral (Table 4). Conversely, CH2O-inclusive models allow silicates and calcite to dissolve if montmorillonite precipitation, or kaolinite and montmorillonite precipitation, occurs. It is logical to assume that the net dissolution of calcite and primary silicate minerals must occur in subglacial environments to account for solute enrichment between supraglacial and bulk meltwaters. CH2O consumption always requires dissolved O2, indicating the simulated organic carbon oxidation is driven by electron transport to dissolved O2 since this is the only available electron acceptor phase in the model. Organic carbon mineralization will very likely be microbially mediated, as it is an otherwise exceptionally slow abiotic process (Konhauser, 2006). Modeling based on the average seasonal bulk meltwater composition therefore suggests that carbonate and incongruent silicate weathering, driven by organic carbon oxidation and pyrite oxidation by dissolved O2, occurs at some point in the subglacial hydrological system at Haut Glacier d'Arolla. While simple hydrolysis may account for some of the weathering, as is suggested in current subglacial biogeochemical models (Tranter et al., 2002), organic carbon oxidation and sulfide oxidation with dissolved O2 are required at some point in the subglacial environment. Therefore, so far, PHREEQCi modeling is in agreement with current subglacial biogeochemical models (Tranter et al., 2002, 2005).
To constrain the number of weathering scenarios to produce only those that include carbonate and silicate weathering driven by organic carbon oxidation and sulfide oxidation with dissolved O2, in the INVERSE_MODELLING toggle of PHREEQCi, primary silicates, calcite, and CH2O were set to dissolve, and ferrihydrite and clays set to precipitate as before. These constraints significantly reduce the number of weathering scenarios from the original 47 to 18 (Table 4). All these weathering scenarios are comprised of CH2O and sulfide oxidation by dissolved O2 and incongruent aluminosilicate dissolution.
Only five of the main primary silicate phases (albite, biotite, muscovite, microcline, tremolite, anorthite, quartz) may be dissolved in any of the 18 individual weathering scenarios, indicating that under this model configuration the simultaneous weathering of all the main catchment minerals cannot account for the change in water composition from supraglacial to bulk meltwaters. Since the model consistently suggests organic carbon and pyrite oxidation–driven weathering, the feasibility of the 18 weathering scenarios can be given some quantitative assessment with consideration of the relative reactivity of the dissolving carbonate and silicate minerals. Glacial environments exhibit high physical weathering rates, which constantly supply freshly comminuted mineral surfaces. Under these conditions, weathering is largely chemically limited where the most reactive minerals, but not necessarily the most abundant, will be likely to provide the majority of solute (Anderson et al., 1997, 2000; White et al., 2005). Therefore, mineral-specific dissolution rates (after Lasaga et al., 1994; Langmuir, 1997) for minerals present in the Haut Glacier d'Arolla catchment were used to consider the most likely weathering scenarios. The weathering series order for carbonate and silicate minerals present in the Haut Glacier d'Arolla catchment and used in the model is anorthite ≈ calcite > biotite > tremolite > muscovite > quartz > microcline. The dissolution rate data used is derived from laboratory dissolution experiments at pH 5 and 25 °C. Lower water temperatures generally reduce weathering rates, and would therefore be expected in subglacial environments (∼0.1°C), although the effect is less pronounced for biotite (Anderson, 2005; White et al., 1999). However, these data are simply used to indicate the relative dissolution rates under fixed conditions, to provide some estimate as to the most likely weathering scenarios.
Anorthite dissolution is required in all the weathering scenarios, while calcite-free weathering scenarios are generated in 9 of the 18 scenarios (Table 4). Therefore, PHREEQCi suggests that anorthite is required to account for the enrichment of Ca between supraglacial and bulk meltwaters, as has been suggested in current models of subglacial biogeochemical weathering models (Tranter et al., 2002, 2005) but has rarely been quantified (cf. Anderson et al., 1997). Anorthite is highly reactive (2.00 * 10−9 mol m−2 s−1 [Lasaga et al., 1994]), so the model assertions seem likely. However, the absence of calcite weathering is highly unlikely even if in trace quantities in the bedrock, reflecting its high dissolution rates (8.43 * 10−10 mol m−2 s−1 [Lasaga et al., 1994]). Indeed, all current models of chemical weathering in glacierized and mountainous catchments indicate the dominance of carbonates as a solute source even when in trace quantities in the bedrock, owing to the chemically limited weathering regimes which operate in such environments (Anderson et al., 1997, 2000; Jacobson et al., 2002, 2003; Tranter et al., 2002; White et al., 2005). This suggests the 9 calcite and anorthite inclusive weathering scenarios (Table 4) are most feasible for the subglacial environment.
The consistent requirement of albite and anorthite dissolution in each of the weathering scenarios indicates feldspar weathering is essential for the evolution of bulk meltwater chemistry at Haut Glacier d'Arolla. Albite is the only likely major source of Na. Mica (muscovite or biotite) dissolution occurs in 7 of the 9 anorthite and calcite inclusive weathering scenarios, with K-feldspar (microcline) dissolution only occurring in 2 of the 9 weathering scenarios. The apparent dominance of biotite weathering over K-feldspar weathering appears likely as biotite dissolution rates (5.6 * 10−12 mol m−2 s−1 [Lasaga et al., 1994]) are much higher than K-feldspar (3.16 * 10−13 mol m−2 s−1 [Lasaga et al., 1994]). However, muscovite dissolution rates are 3.7 times lower than that of K-feldspar (8.51 * 10−14 mol m−2 s−1). This suggests the 5 biotite inclusive weathering scenarios (with calcite and anorthite dissolution; WS 41, 44, 45, 46, 47; Table 4) are more likely than muscovite weathering scenarios. Indeed, the importance of biotite weathering has been previously demonstrated from glacial moraine chronosequences (Blum and Erel, 1997), and mass balance weathering schemes (Anderson et al., 2000). Tremolite exhibits dissolution rates on the same order as biotite (1.9 *10−12 mol m−2 s−1 [Lasaga et al., 1994]); therefore, tremolite inclusive weathering scenarios, with combined biotite, anorthite, and calcite weathering would seem the most likely weathering scenarios, based upon mineral dissolution rates. This is achieved only in WS41. We conclude that PHREEQCi is therefore able to deconvolute silicate weathering more than can easily be achieved through traditional longhand mass balance calculations (Tranter et al., 1993; Anderson et al., 2000), owing to the ability of the numerical code to solve balance equations simultaneously, rather than sequentially. It is worth noting that the model suggests the formation of secondary phases which have notable ion exchange properties, including montmorillonite and Fe-oxyhydroxies. We cannot account for this in the SMB model, but ion exchange with such minerals may well be expected to have some influence on water chemistry (Souchez and Lorrain, 1975), such that solute ratios measured in the meltwater may not only reflect the stoichiometric release from the mineral lattice, but also ion exchange with these secondary phases (Wadham et al., 1998). However, like all previous mass balance studies, we have to assume that ion exchange is in steady state as very little is known about the ion exchange pool in aqueous systems.
Testing The Combined Influence of Ch2o and CO2
Current models of subglacial biogeochemical weathering suggest that atmospheric CO2-driven weathering is relatively insignificant, but the potential for combined organic carbon oxidation and atmospheric CO2 dissolution as a DIC source has been given little quantitative consideration (cf. Wadham et al., 2004). Some CO2-driven weathering may be expected in turbid, well-mixed waters which are open to the atmosphere and are rapidly transported through englacial and subglacial channels, and also in connected channel marginal locations. In order to test this, CO2(g) was added to defined phases in the model and constrained to dissolve, in conjunction with CH2O. Only two additional weathering scenarios were generated with combined organic carbon oxidation by dissolved O2, and CO2(g) dissolution (WS48 and 49; Table 4). WS48 did not include calcite weathering, and is therefore unlikely to be feasible. However, WS49 did include calcite weathering, along with sulfide oxidation by dissolved O2, and incongruent aluminosilicate dissolution, with albite, biotite, tremolite, and anorthite weathering, and is like WS41 except for the distribution of DIC between CH2O and CO2(g). This weathering scenario suggests that CH2O oxidation is dominant over CO2 dissolution at a mass ratio of ∼10:1. Therefore the results from PHRREQCi support current models of subglacial biogeochemical weathering, and indicate that even if some CO2-driven weathering occurs, it is relatively insignificant to that driven by organic carbon oxidation.
Testing Anoxic Subglacial Weathering
It has been suggested in current subglacial biogeochemical models that bacterial catalysis of organic carbon oxidation by dissolved O2, and sulfide oxidation by dissolved O2, may drive subglacial environments towards anoxia (Tranter et al., 2002, 2005). Once anoxia has developed, methanogenesis has been proposed (Tranter et al., 2005). Indeed, stable isotope evidence from bulk meltwaters from a polythermal High Arctic glacier indicates that reduction of SO42− was occurring and anoxic conditions apparent (Wadham et al., 2004). In order to test the effect of anoxia upon the weathering mechanisms suggested by PHREEQCi, O2(aq) was omitted from the model and CH4 was included as a potential phase to investigate whether methanogenesis could account for observed changes in DIC. The model produced 9 weathering scenarios, all of which included calcite and anorthite dissolution (Table 5). CH2O reduced to CH4, producing H+, which could account for changes in DIC between supraglacial and bulk meltwaters, and was required if no O2(aq) was available.
Results of model testing under anoxic conditions, between supraglacial and bulk meltwaters, and between supraglacial and in situ waters collected from the distributed system. Each line is a single weathering scenario (WS) that can account for the changes in solution chemistry between supraglacial and bulk waters or in situ waters. Units = mass of dissolving phases (mmol L−1) transferred between input supraglacial water and bulk or in situ waters.
The SMB model was also run using supraglacial input waters evolving to in situ subglacial waters collected directly from the base of boreholes, which appeared to be representative of waters collected from the distributed system (Mitchell et al., 2006). In such environments, anoxia and methanogenesis would be most likely to occur due to poor hydraulic connectivity to oxygenated waters (Tranter et al., 2005). These waters exhibit major ion concentrations that are on average 320% the concentration of bulk meltwaters than measured in all bulk meltwaters and have low suspended sediment concentrations, suggesting these waters have had a low water velocity and a long residence time within the distributed system (Tranter et al., 1997; Mitchell et al., 2006). Using the mean composition of supraglacial meltwaters evolving to the mean composition of these waters collected from the distributed system, methanogenesis was able to account for the changes in water chemistry, producing 10 weathering scenarios, providing some quantitative evidence to support the hypothesis that methanogens operate in areas of the bed where there is poor hydraulic connectivity (Table 5) (Tranter et al., 2005). Interestingly, all these weathering scenarios required anorthite weathering, but only 5 of the 10 weathering scenarios included calcite weathering, suggesting Ca-silicate weathering may dominate over calcite as a Ca source in the distributed system. This is likely as longer water residence times would allow Ca-silicate minerals, which are in higher abundance but are generally less reactive, to provide proportionally more solute (Anderson et al., 2000; Drever and Zobrist, 1992).
The modeling therefore demonstrates that methanogenesis can account for the evolution of bulk waters, which integrates both waters from the distributed and channelized system, and also water collected directly from the distributed system. However, methanogenesis is unlikely to occur in subglacial channels or channel marginal locations due to the constant influx of oxygenated waters. Therefore, this demonstrates how viable weathering scenarios from PHREEQCi, while mathematically feasible with regard reaction stoichiometry, have to be considered within the hydrological context of the system.
The suggestion that Fe3+ derived from the reduction of Fe-(oxy) hydroxides or Fe-silicate minerals can act as an oxidizing agent for sulfide minerals under anoxic conditions (Tranter et al., 2002, 2005) was also tested. Recent anoxic incubation of glacial sediments has demonstrated the occurrence of microbially induced Fe reduction by the measurement of Fe2+(aq) (Montross, 2007). However, it remains unclear if aqueous Fe will exist sufficiently in the Fe3+ form under anoxic conditions, since dissolved Fe2+ will evidently predominate under anoxic conditions. We therefore undertook speciation modeling on the mean solution composition for in situ water collected from the distributed system, where total dissolved Fe was measured. By the nature of the model, this estimates abiotic equilibrium speciation, and will thus indicate whether Fe3+ in waters from the subglacial environment can be produced abiotically. The redox potential was varied from the oxic conditions of pE = 7 used in the model thus far, down to an estimated 0 (suboxic) and −5 (anoxic) (Langmuir, 1997) (Fig. 5). Aqueous Fe2+(aq) increasingly dominated Fe speciation below pE 4 through to heavily reduced conditions of pE −5, reaching concentrations <8 µmoles L−1. Fe(III)(aq) species [Fe(OH)2+, Fe(OH)3+, Fe(OH)4−] decreased in concentration from combined concentrations of <8 µmoles L−1 under oxic conditions of pE 7 down to 0 µmoles L−1 at pE 2 and below. This therefore reinforces the suggestion that abiotic processes cannot account for the generation of aqueous Fe3+ (Tranter et al., 2002, 2005), and suggests that microbial oxidation of Fe2+ to Fe3+ is required if Fe3+ is indeed an important oxidizing agent in subglacial environments. This is highly likely given that the reducing power of Fe2+ increases dramatically at pH values higher than 2–3, owing to the formation of ferric hydroxy and oxyhydroxy compounds. The standard redox potential of Fe3+/Fe2+ (E0 = +0.77 V) is relevant only under acidic conditions, whereas at pH 7, typical in in situ and bulk meltwaters, the couples Fe(OH)3/Fe2+ (E′0 = −0.236 V) or Fe(OH)3 + HCO3− FeCO3 (E′0 = +0.200 V) prevail (Widdel et al., 1993). Fe2+ should thus be able to donate electrons readily, facilitating the energy requirements of subglacial Fe oxidizing microbial populations. It is also worth considering that intergranular microenvironments in subglacial sediments will exhibit different aqueous chemistries than measured in the bulk solution in in situ environments. Such differences can often account for differences between biogeochemical mechanisms measured in the laboratory and field (Geesey and Mitchell, 2008). Research into biogeochemical transformations at the micro-scale in subglacial sediments should therefore be pursued.
Comparison of Modeling-Defined Changes in Subglacial Weathering Processes Over Seasonal Timescales
The biogeochemical evolution of bulk meltwaters are the product of reactions occurring both in the distributed and channelized system. The modeling so far demonstrates that while anoxic conditions can account for the evolution of waters collected directly from the distributed system, anoxia is unlikely to predominate over the majority of the bed at the height of the ablation season, but rather in areas with poor hydraulic connectivity. A constant supply of oxygenated waters into the subglacial system at this time suggests there will be a net consumption of dissolved oxygen. Therefore we consider the seasonal dynamics of selected weathering scenarios from the previous modeling phase, where net subglacial weathering reactions from supraglacial to bulk meltwaters are driven by organic carbon and pyrite oxidation by dissolved O2. The seasonal time series bulk meltwater composition from the ablation season (Fig. 3) (Mitchell et al., 2001; 2006) were run through PHREEQCi, using concurrent input supraglacial water compositions with the aforementioned phase dissolution and precipitation constraints.
Based upon mineral dissolution rates, we consider some of the apparently most feasible weathering scenarios, WS41, WS43, and WS49, identified in Sections TESTING VIABLE CHEMICAL WEATHERING REACTIONS IN SUBGLACIAL ENVIRONMENTS USING PHREEQCi and TESTING THE COMBINED INFLUENCE OF CH2O AND CO2; (Table 4). These weathering scenarios are chosen because they exhibit markedly different seasonal weathering dynamics to each other, while WS49 allows us to investigate seasonal changes in organic carbon oxidation versus CO2-driven weathering. This allows us to investigate whether the modeled variations in chemical weathering are consistent with known seasonal changes in the subglacial hydrological system, and changes in meltwater routing within this system at Haut Glacier d'Arolla (Nienow et al., 1998), and other alpine glaciers which undergo the same general seasonal changes in subglacial hydrology.
Carbonate and Silicate and Weathering
The assumptions of current models of subglacial biogeochemical weathering suggest that most Ca2+ derives from CaCO3 dissolution, with the rest derived from Ca-silicates (Tranter et al., 2002, 2005). Other studies in glacierized and alpine catchments (Anderson et al., 1997; 2000) have also suggested that the chemically limited weathering of trace carbonate minerals dominates Ca2+ sources in glacial and alpine environments underlain by igneous and metamorphic bedrock (e.g. Drever and Zobrist, 1992; Brown, 2002). Similarly, WS43 suggests the weathering of trace calcite dominates on seasonal timescales, as shown by the mass of individual dissolving phases transferred between supraglacial and bulk meltwaters (Fig. 6; Table 6). Here the proportion of Ca2+ derived from calcite is comparable to previous studies (e.g. Anderson et al., 2000), ranging between 70 and 100% on average over the ablation season (Table 7).
Median mass of dissolving phases (mg L−1) transferred between input supraglacial meltwater and bulk meltwaters under modeling weathering scenarios (WS) 43, 41, and 49 over the ablation season. Mass transfer of each phase also expressed as percentage of total mass transfer.
Average percentage of species produced or consumed by the dissolution or precipitation of individual phases between input supraglacial meltwater and bulk meltwaters under weathering scenarios (WS) 43, 41, and 49 over the ablation season.
WS43 demonstrates there is very little change in the proportion of carbonate and silicate weathering with changes in bulk meltwater discharge over the ablation season (Fig. 6). This suggests that the proportion of carbonate and silicate weathering may not be significantly affected by the evolution of the drainage system on the seasonal timescale. This seems unlikely as such large seasonal changes in subglacial water routing and water residence times should impart a large impact on the weathering of lithologies with differing geochemical reactivity (Tranter et al., 2005). Indeed, consideration of WS41 and 49 seems more plausible than WS43 with regard to seasonal changes in subglacial flow routing and mineral weathering kinetics. Here, calcite dissolution and incongruent silicate weathering provide an approximately equal mass proportion due to high transfers of anorthite, both of which provide ∼40% of Ca2+, and tremolite, which provides the remaining ∼20% of Ca2+, on average, over the ablation season (Tables 6 and 7). Unlike WS43, calcite becomes proportionally more important while anorthite, tremolite, and total incongruent silicate dissolution becomes proportionally less important as discharge increases, such that calcite and anorthite weathering is approximately equal at the height of the ablation season (Fig. 6), supplying approximately equal amounts of Ca2+ (Fig. 7). Calcite exhaustion is unlikely at higher discharges, and experimental investigations suggest the rate of calcite dissolution is greatest at higher suspended sediment concentrations, as generated at higher discharges (Brown et al., 1996b). Indeed, chemically limited calcite weathering should dominate and be enhanced relative to silicate weathering at higher discharges, where the most reactive but not necessarily most abundant minerals will provide the majority of solute because sediment transport times are rapid relative to mineral weathering rates. Therefore, the observed increase in the proportion of calcite weathering and decrease in the proportion of incongruent silicate dissolution with increasing discharge in WS41 and 49 is most convincing with regard to flow routing and mineral weathering kinetics, and reflects the seasonal expansion of the channelized system.
The modeling suggestion that Ca-silicate weathering may predominate in the distributed system may also be apparent in other subglacial environments. This may be of particular importance to polythermal glaciers, which override bedrock containing trace amounts of calcite, and where the development of channelized systems is limited, and most waters are routed through a distributed system (Wadham et al., 1998; Hodson et al., 2000). Here, silicate minerals, which are more abundant but generally less reactive than carbonates, have sufficient time to dissolve and may ultimately provide more solute than trace reactive minerals (Drever and Zobrist, 1992).
Relative Importance of Sulfide Oxidation, Organic Carbon Oxidation, and CO2 Dissolution as Proton Sources
Pyrite oxidation by dissolved O2 has been shown to preferentially occur in distributed subglacial environments (Tranter et al., 2002). We would similarly expect organic oxidation to proceed preferentially in basal sediments in the distributed system, where microbial activity will mostly be limited to, because the residence time of water is sufficiently long, and hydraulic shear stress will be low, allowing microbial populations to maintain their existence (Tranter et al., 2005). We would therefore expect to see both sulfide oxidation and organic carbon oxidation provide a greater proportion of solute at lower bulk discharges, when water from the distributed system will account for proportionally more of the bulk discharge than at higher discharges. This is the case in WS41 (Figs. 6 and 7). However, in WS43 the proportion of organic carbon oxidation actually increases at higher discharges, which seems unfeasible from a hydrological perspective. This suggests that WS41 is most feasible from a hydrological and biogeochemical perspective, as was also suggested by the seasonal dynamics of carbonate and silicate weathering (Section Carbonate and Silicate and Weathering). Here organic carbon oxidation accounts for 76% of DIC on average over the melt season, with the remaining 24% derived from calcite dissolution. In WS49, which resembled WS41 apart from it being the only weathering scenario to include CO2 as well as CH2O as a DIC source, CO2 provided only 8.5% of DIC, and organic carbon oxidation accounting for 65% DIC, on average over the melt season, with calcite providing 26.5% of DIC (Table 7; Fig. 7). The source of carbon dioxide, required in WS49, may be atmospherically derived from dissolution into turbulent waters in subglacial channels or from gas bubbles trapped in subglacial ice, since microbial sources have already been accounted for in the model by the oxidation of organic carbon.
Both of these estimates of organic carbon oxidation are far higher than the only previous estimate made of subglacial organic carbon mineralization, at Finsterwalderbreen, Spitzbergen, by Wadham et al. (2004), who calculated from δ13C isotope mass balance calculations that ∼10% of the total annual DIC flux was due to microbial organic carbon mineralization. The degree of organic carbon oxidation in different glacierized catchments will be controlled by concentration and availability in the subglacial environment (variations in the extent of organic-rich sediments and soil formation prior to overriding by the glacier, the composition of the bedrock, its kerogen content, the abundance of microbial necromass, wash-in of organic carbon from the glacier surface), the physical conditions in the subglacial hydrological system (water velocity, particle size), and the abundance of chemoheterotrophic microbial populations which may catalyze organic carbon oxidation. At Finsterwalderbreen, the organic carbon content of the bedrock ranged from <0.05 to 0.27 wt% in carbonates, 0.11 to 0.47 wt% in siltstones and sandstones, and 1.2 to 2.3 wt% in shales. At Haut Glacier d'Arolla, the organic carbon content of proglacial sediments is 0.22 wt%, and dissolved organic carbon (DOC) was 0.2 mg L−1 in bulk meltwaters (Montross, 2007). It would seem likely from these figures that the abundance of organic carbon would on average be lower at Haut Glacier d'Arolla than Finsterwalderbreen, which suggests that our estimates of organic carbon oxidation are too large. However, the ingress of organic carbon into the subglacial environment is likely to be much larger at Haut Glacier d'Arolla than Finsterwalderbreen due to the higher degree of surface melt and the greater number of crevasses and moulins which will allow in-washing of extra-glacial soil and other non-glacial organic carbon deposited on the glacier surface. Therefore the higher degree of organic carbon turnover may be expected at Haut Glacier d'Arolla and other alpine glaciers than under polar glaciers, where surface melt and water routing from supraglacial to subglacial environments is more limited. At this time it is hard to determine whether this is the case, and consequently further studies are required to ascertain the controls and kinetics of organic carbon oxidation in a range of subglacial systems, and the sources of organic carbon, and whether new organic carbon is produced in subglacial environments by chemoautotrophy (Priscu et al., in press). Nonetheless, SMB modeling reinforces that chemical weathering in subglacial environments is distinctive from that in most other earth surface environments, where weathering is driven by atmospheric CO2 or biologically derived CO2 in soils (Tranter et al., 2002).
Speciation–mass balance models applied to the subglacial hydrological system at Haut Glacier d'Arolla demonstrate such models provide a potentially useful tool for investigating geochemical and biogeochemical weathering processes in a range of subglacial environments. They allow all stoichiometrically and thermodynamically feasible weathering scenarios to be tested simultaneously quickly and easily, unlike previous mass balance calculations which require mass balance equations to be solved sequentially, ad hoc, or in a time consuming longhand mass balance matrix.
These hypothetical geochemical and biogeochemical weathering scenarios suggest a number of key features characterize subglacial weathering processes at Haut Glacier d'Arolla, which are likely to apply in many subglacial systems owing to the same broad physical, chemical, and biological conditions. Organic carbon oxidation and sulfide oxidation by dissolved O2, coupled to carbonate dissolution and incongruent silicate dissolution could account for seasonal changes in meltwater chemistry, supporting current subglacial biogeochemical weathering models. All carbonate ions in glacial meltwaters cannot realistically be derived from carbonate minerals alone, requiring additional carbon from organic carbon. CO2(g) dissolution was not required in any weathering scenarios reinforcing the apparent unimportance of CO2(g)–water transfers in subglacial environments. Silicate weathering is likely to be dominated by biotite and Ca-feldspar weathering, owing to their relatively high dissolution rates. Indeed, the model also reinforces the suggestion in some glacierized catchments that Ca-feldspar dissolution is an important solute source. Here we quantify under the most likely weathering scenario that Ca-feldspar provides approximately equal amounts of Ca2+ as calcite. We also demonstrate that microbial oxidation of Fe2+ to Fe3+ is required if Fe3+ is indeed an important oxidizing agent in anoxic subglacial environments. This investigation indicates the utility of SMB models applied to subglacial hydrological systems for determining geochemical and biogeochemical processes and interpreting water quality. The modeling reinforces that chemical weathering in subglacial environments is distinctive from that in most other earth surface environments, where weathering is driven by atmospheric CO2 or biologically derived CO2 in soils.
This work was supported by a University of Wales Aberystwyth Studentship (AM) and NERC Grant No. GR3/11216. Many thanks to Mark Skidmore and Scott Montross for TOC/DOC data for Haut Glacier d'Arolla, Ron Fuge for valuable contributions to this investigation at an early stage, Grande Dixence S.A for the provision of discharge data, and to Ian Willis, Doug Mair, Bryn Hubbard, Pete Nienow, Darrel Swift, Urs Fisher, Becky Goodsell and Ian Cochrane. We would especially like to thank two anonymous referees and Suzanne Anderson for their valuable comments which have significantly strengthened this work.