The Greenland Ice Sheet as a hot spot of phosphorus weathering and export in the Arctic

The contribution of ice sheets to the global biogeochemical cycle of phosphorus is largely unknown, due to the lack of field data. Here we present the first comprehensive study of phosphorus export from two Greenland Ice Sheet glaciers. Our results indicate that the ice sheet is a hot spot of phosphorus export in the Arctic. Soluble reactive phosphorus (SRP) concentrations, up to 0.35 µM, are similar to those observed in Arctic rivers. Yields of SRP are among the highest in the literature, with denudation rates of 17–27 kg P km−2 yr−1. Particulate phases, as with nonglaciated catchments, dominate phosphorus export (>97% of total phosphorus flux). The labile particulate fraction differs between the two glaciers studied, with significantly higher yields found at the larger glacier (57.3 versus 8.3 kg P km−2 yr−1). Total phosphorus yields are an order of magnitude higher than riverine values reported in the literature. We estimate that the ice sheet contributes ~15% of total bioavailable phosphorus input to the Arctic oceans (~11 Gg yr−1) and dominates total phosphorus input (408 Gg yr−1), which is more than 3 times that estimated from Arctic rivers (126 Gg yr−1). We predict that these fluxes will rise with increasing ice sheet freshwater discharge in the future.


Introduction
Phosphorus (P) is a bioessential macronutrient [Ruttenberg, 2014]. The supply of P is known to influence ecosystem productivity and biogeochemical cycles in marine and terrestrial ecosystems. It is believed to be the ultimate limiting nutrient on geological timescales, in part because it cannot be fixed from the atmosphere [Broecker and Peng, 1982;Elser et al., 2007;Ruttenberg, 2014;Tyrrell, 1999] and may currently be a limiting nutrient in some regions [Wu et al., 2000]. Many aspects of the current global P cycle are poorly understood [Filippelli, 2008], especially terrestrial fluxes of P to the oceans, the principal speciation of P exported and the associated bioavailability of fluxes [Slomp, 2011].
Studies of the transport of P to the ocean have so far focused on riverine runoff, atmospheric deposition, and submarine groundwater discharge. Of these, riverine runoff is the dominant flux [Meybeck, 1982;Moore et al., 2013;Ruttenberg, 2014;Slomp, 2011]. Recent research indicates that glaciers and ice sheets export large quantities of nutrients, such as Fe, N, and dissolved organic carbon, to downstream ecosystems [Bhatia et al., 2013;Death et al., 2014;Hawkings et al., 2015Hawkings et al., , 2014Hood et al., 2015Lawson et al., 2014;Raiswell and Canfield, 2012;Schroth et al., 2014]. Only two studies have suggested a role for glaciers in the global P cycle [Föllmi et al., 2009;Hodson et al., 2004]. These have focused on small glacial systems, which may have lower chemical weathering intensities than larger ice sheet catchments [Wadham et al., 2010]. To date, ice sheets have largely been overlooked as P sources to the oceans.
The Greenland Ice Sheet (GrIS) is the second largest ice sheet on Earth, covering an area of approximately 1.7 × 10 6 km 2 [Knight, 1999]. From 2000 to 2010, it discharged~1000 km 3 of freshwater (~40% as runoff and~60% as solid ice) into the surrounding oceans each year [Bamber et al., 2012;Tedesco et al., 2013]. This is likely to increase with climatic warming .
Two distinct environments likely dominate phosphorus export from the GrIS. The surface is the principal source of snow and ice melt. These are relatively dilute, containing only low concentrations of rock-derived surface and routed to the bed [Bartholomew et al., 2010Chandler et al., 2013;Chu, 2014;Cowton et al., 2013;Sole et al., 2011]. The second source of P from the GrIS is the subglacial environment. A drainage system similar to that found under temperate glaciers is believed to be present [Bartholomew et al., 2010;Chandler et al., 2013;Chu, 2014;Sole et al., 2011], and the environment is characterized by high rock-water ratios and microbially mediated chemical weathering of ground bedrock [Tranter et al., 2002]. Meltwater drainage through the subglacial system adds solutes to waters and is thought to be an important source of nutrients [Hawkings et al., 2015]. Periodic flushing of concentrated stored waters by supraglacial lake drainage events has been observed , and meltwaters emerging from beneath the ice sheet carry large quantities of fine-grained suspended sediment with concentrations >1 g L À1 [Cowton et al., 2012]. These sediments have been demonstrated to carry labile nutrients [Hawkings et al., 2015Lawson et al., 2014;Schroth et al., 2014].
This study investigates the importance of the GrIS in the P cycle. Hydrological and geochemical data are presented from a large (~600 km 2 ) and a small (~36 km 2 ) GrIS catchment of similar geology. We utilize highresolution hydrological data sets alongside seasonal time series of dissolved and sediment-bound P species. We aim to (i) present the seasonal variation of P concentrations, (ii) determine the principal processes influencing the P concentrations, (iii) calculate P denudation rates to assess the intensity of glacial P weathering, and (iv) evaluate the importance of the ice sheet on the Arctic P cycle.

Study Areas
Leverett Glacier (LG; 67.06°N, 50.17°W; Table 1 and Figure 1a) is a large polythermal-based, land-terminating outlet of the GrIS. The glacier stretches~80 km from the ice sheet margin and has a hydraulically active catchment area of >600 km 2 [Cowton et al., 2012;Palmer et al., 2011]. Mean annual summer discharge (2009)(2010)(2011)(2012) ranged from 110 to 200 m 3 s À1 and annual specific discharge of 1.8-3.7 m. Runoff from LG feeds a large glacial river system (Watson River), which discharges into the Davis Strait via Søndre Strømfjord. The catchment bedrock is composed predominantly of ancient crystalline gneiss/granite, which is the dominant geology underlying the GrIS [Henriksen et al., 2009]. Hydrology is similar to temperate glaciers, with supraglacial meltwater input driving subglacial drainage evolution over the melt season, from a slow, inefficient distributed system to an efficient channelized system . LG was sampled during the record GrIS ablation season of 2012 (11 May to 21 July). Kiattuut Sermiat (KS; 61.2°N, 45.3°W; Table 1 and Figure 1b) is a smaller coastal glacier on the southern tip of the GrIS. It is also believed to have a polythermal basal regime. KS covers an estimated 36 km 2 , is around 16 km in length, and terminates in a proglacial lake. Mean ice velocity is <60 m yr À1 and at least the final few hundred meters of ice is believed to be static ("dead ice") [Rignot and Mouginot, 2012]. Mean annual summer discharge was~25 m 3 s À1 , and specific annual discharge was 6.1 m in 2013 (Table 1). Regional geology is crystalline granite, broadly similar in composition to LG. Smaller complexes of diorite and pyroxenebiotite monozonite and basaltic intrusions are also present [Henriksen et al., 2009]. KS was sampled during the 2013 ablation season (21 April to 11 August).

Hydrological Monitoring
LG and KS meltwater rivers were monitored for discharge, suspended sediment concentrations, and electrical conductivity (EC) at stable bedrock sections (Figure 1) using methods previously described by Bartholomew et al. [2011] (Text S1 in the supporting information). pH was recorded manually when samples were taken using a handheld Hanna® HI9124 meter with standard glass electrode probe at LG and continuously using an ISFET (ion-selective field-effect transistor) Honeywell Durafet® pH sensor at KS. Handheld meters were calibrated daily using National Institute of Standards and Technology-certified pH 7 and 10 buffers (Fisherbrand TM ), and the ISFET sensor was calibrated monthly with low ionic strength pH 4.01 and pH 6.96 buffers (Reagecon). Errors are estimated at ±0.05 pH units.

Sample Collection, Processing, and Storage
Bulk runoff samples were collected at least once daily during the 2012 (LG) and 2013 (KS) ablation seasons ( Figure 1). Collection time was between 10.00 and 12.00 h, and on occasion between 18.00 and 19.00 h. Meltwater was collected in 2 L high-density polyethylene (HDPE) Nalgene® bottles. Samples for dissolved P Global Biogeochemical Cycles 10.1002/2015GB005237 analysis were immediately filtered through a 47 mm 0.45 μm cellulose nitrate filter (Whatman®), mounted onto a polyethersulfone filter unit (Nalgene®) at LG, or a 0.45 μm Whatman® GD/XP PES syringe filter using a PP/PE syringe at KS (Cyanmar®). Filtration was carried out in a designated "clean" sampling tent (LG), or immediately on site (KS; weather permitting). Three replicate samples collected at LG using the filter unit and syringe filter methods showed no significant difference in final measured concentration (±3.6%). Supraglacial meltwater samples at KS were collected in a similar manner. Dissolved P samples were stored in clean 30 mL HDPE bottles (Nalgene®), three times rinsed with filtered sample water, and stored frozen and in the dark at À10°C, using portable field freezers, before transport back to Bristol, where they were kept at À25°C prior to analysis. Freezing is an effective and commonly used preservation method for storing low-level nutrient samples for short periods [Dore et al., 1996]. Samples for suspended sediment-bound P fractions were collected from the same bulk meltwater samples. Around 200-1000 mL of meltwater was filtered through a 47 mm 0.45 μm cellulose nitrate filter (Whatman®), with the sediment retained and stored frozen until analysis. Measurements were based on the molybdenum blue method. Precision was ±3.2% (LG; n = 5) and ±2.9% (KS; n = 7), and accuracy was +1.5% (LG) and À5.7% (KS), as determined from comparison with gravimetrically diluted 1000 mg L À1 PO 4 certified stock standards (Sigma TraceCERT®). Field blanks were below the limit of detection (LoD), 0.01 μM (n = 9).

Particulate Phosphorus
A four-step sequential extraction methodology, adapted from Stibal et al. [2008], was used. This method removes increasingly unreactive, operationally defined, P phases. Sediment was removed carefully by gentle scraping from air-dried filter membranes. Approximately 50 mg of sediment was added to a 2 mL PP microcentrifuge tube. Extraction 1 removes highly labile, "loosely adsorbed" P (MgCl 2 -P). A 1 M MgCl 2 solution (1.5 mL) was added to the microcentrifuge tube, to give a final solution:sediment ratio similar to that used by others [Hodson et al., 2004]. The tubes were placed on a reciprocating shaker for 16 h, centrifuged at 3000 rpm for~10 min, and the supernatant extracted and filtered into a clean 2 mL centrifuge tube using a 1 mL syringe (PP/PE) and a 0.45 μm Whatman® Puradisc syringe filter (PP membrane). The sediment was subsequently shaken with 1.5 mL of Milli-Q water for 2 h, the supernatant extracted, filtered as above, and added to the initial extract. Extraction 2 removes Fe-and Al-bound P, which is believed to be moderately labile (NaOH-P). As this phase is calibrated to the growth of algae cultures, it is often referred to as "algal available" [Dorich et al., 1985]. 1.5 mL of 0.1 M NaOH was added to the retained sediment, and the same process as per extraction 1 was followed. Extraction 3 removes Ca-and Mg-bound P, which is more refractory and less reactive (HCl-P). Samples were treated as above, but with a 1 M HCl solution added to the retained sediment. The final extraction, extraction 4, was used to determine "residual" and organic-bound P (Res-P). The remaining sediment was transferred to a 10 mL glass digestion tube, to which 5 mL of a sulfuric acid/potassium persulfate solution was added [Jeffries et al., 1979]. Samples were autoclaved for 30 min at 120°C and 100,000 Pa, allowed to cool, and then filtered into clean PP microcentrifuge tubes as above. Samples were kept refrigerated at <4°C for no longer than 24 h prior to analysis. Solutions were analyzed using the SRP method described in section 2.4.1 with matrix-matched standards. Precision was found to be ±0.7%, ±0.9%, ±0.2%, and ±0.2%, and accuracy was +0.2%, À0.4%, +2.3%, and À2.3% based on six replicate standards for extractions 1, 2, 3, and 4, respectively. All extractions were diluted with Milli-Q water before analysis to bring them within the instrumental determination range. Extraction 1 was diluted 1:3, extraction 2 1:10, extraction 3 1:100, and extraction 4 1:15. Dry weight corrections were made based on the mean dry weight corrections of six sediment samples from LG (±3.9%) and twenty one samples from KS (±3.6%).

Principle Components Analysis
Principle components analysis (PCA) was performed using SPSS® Statistics software (IBM®). PCA is a data transformation technique used to reveal underlying relationships in complex data sets. PCA was conducted using Direct Oblimin rotation with Kaiser Normalization.
2.6. Catchment and Ice Sheet P Fluxes 2.6.1. Dissolved Fluxes of P Catchment meltwater fluxes for the sampled years (2012 and 2013) were calculated from the seasonal discharge record following the method of Cowton et al. [2012]. Total catchment runoff was derived by summing discharge from each 5 or 10 min data point (frequency at which discharge was logged) over the monitoring period (LG = days 134-251; KS = days 112-223). Dissolved P fluxes were calculated from the product of the water flux and the discharge-weighted mean concentrations of SRP and DOP (see section 2.4), following previous work [e.g., Meybeck, 1982;Hawkings et al., 2014]. A P flux range was also calculated from total water flux and minimum and maximum recorded P concentrations.
Dissolved P fluxes from the GrIS were estimated from the minimum, discharge-weighted mean and maximum concentrations from our LG data set, multiplied by the mean 2000-2012 GrIS water flux taken from Tedesco et al. [2013], estimated as 437 (± 97) km 3 yr À1 . A crude estimate of future fluxes was calculated from Tedesco et al. [2013] discharge data from 2012, a record melt year, when the meltwater flux was estimated to be 665 km 3 . We use P data from LG, rather than KS, since this is more likely to represent the hydrology of other large outlet catchments that dominate meltwater export from the GrIS (see Text S1) Chu, 2014;Hawkings et al., 2014].

Particulate Fluxes of P (PP)
LG and KS suspended sediment loads were calculated from suspended sediment concentrations (section 2.2) multiplied by discharge at each logged time point [Cowton et al., 2012;Hawkings et al., 2014]. A minimum, discharge-weighted mean, and maximum concentration for each PP phase (μg g À1 dry sediment; section 2.4.2) was combined with sediment flux data to derive fluxes. Global Biogeochemical Cycles

10.1002/2015GB005237
Sediment fluxes from the GrIS are poorly constrained but are likely globally significant [Cowton et al., 2012;Hay, 1998;Hudson et al., 2014;Jeandel and Oelkers, 2015;Lisitzin, 2002]. We use minimum (0.643 g L À1 ), dischargeweighted mean (1.109 g L À1 ) and maximum (3.876 g L À1 ) suspended sediment concentrations measured at LG to derive estimates of sediment fluxes, as per Hawkings et al. [2014]. GrIS sediment yields are calculated by multiplying these values with the meltwater flux from Tedesco et al. [2013]. This method generates sediment fluxes of 0.28, 0.48, and 1.69 Gt yr À1 . We combine these sediment fluxes with the discharge-weighted mean PP phases from LG to derive a range of estimates (see section 2.6.1) [Hudson et al., 2014].

Glacier Hydrology
The hydrology of LG is well studied Chandler et al., 2013;Cowton et al., 2013]. Seasonal evolution of efficient drainage allows progressively more remote regions of the bed to be accessed. A large, well-documented, melt event was observed at LG during 2012, when over 90% of the GrIS surface was at melting point (day 194; Figure 2a) [Nghiem et al., 2012], corresponding with a record discharge of >800 m 3 s À1 . The melt season is interspersed with meltwater pulse events (E1-E6; Figure 2a) where discharge, suspended sediment, and EC spikes occur in response to surface lake drainage events, forcing stored waters from the glacier bed . Suspended sediment concentrations are high throughout, with a mean concentration of >1.1 g L À1 (Table 1) [Cowton et al., 2012].
KS discharge increases with supraglacial snow and ice melt throughout the ablation season, similar to LG and other temperate glaciers studied [Collins, 1979;Swift et al., 2005]. A spike in suspended sediment concentration, coupled with a sharp decline in EC and major ion concentration, accompanies a large increase in discharge and pH around day 157 ( Figure 2b). This indicates the evolution of an efficient drainage system able to route water to the margin more rapidly [Mair et al., 2003], as the waters exiting the glacier are representative of those with a lower residence time (less concentrated, higher pH). Supraglacial lakes in Greenland do not form at the elevation range on which KS lies (~100-650 m above sea level); hence, there are no associated meltwater pulse events, and seasonal trends largely reflect incident supraglacial melt rates. Suspended sediment concentrations are generally lower than other glaciers, suggesting that settling in the proglacial lake may be a sink of sediments [Bagshaw et al., 2014;Liermann et al., 2012] and/or that KS glacier physical erosion rates are low due to the slow velocity (Table 1).

Dissolved Phosphorus Concentrations in Meltwaters
Elevated SRP concentrations are closely aligned with the flushing of stored waters from the subglacial system by supraglacial lake drainage events and subsequent drainage of subglacial waters from newly connected parts of the ice sheet bed (E2, E3, E5, and E6; Figure 3a). Maximum concentrations of SRP (up to 0.35 μM) were recorded toward the end of the monitoring period. These are associated with the E5 and E6 events, the 2012 record melt event, and waters with a pH >9 (Figures 2a, 3a, and 4a). DOP concentrations at LG are significantly lower than SRP (Table 2 and Figure 3). DOP contributes <10% of the dissolved flux, and there is little temporal trend, with values Global Biogeochemical Cycles 10.1002/2015GB005237 close to those of the supraglacial meltwaters we measured (Table 2), indicating a supraglacial source [Stibal et al., 2008].
KS discharge-weighted mean SRP concentrations are less than LG (0.09 versus 0.23 μM; Table 2 and Figure 3b), despite the generally more concentrated waters (EC of 26.6 versus 8.5 μS cm À1 ; Table 1). There is little temporal trend, unlike major ion concentrations (Figure 2b), and no significant loading of SRP with other variables (Figure 4). Bulk meltwater DOP concentrations are similar to supraglacial DOP concentrations, indicating a supraglacial source, as with LG.

Particulate Phosphorus Concentrations in Meltwaters
The mean PP content of LG suspended sediments (850 μg g À1 ) is similar to that measured in shield rocks (Table 3) [Hodson et al., 2004;Porder and Ramachandran, 2013]. The majority (>90%) of the PP is HCl extractable and is believed to be predominantly fluorapatite. Labile PP (MgCl 2 + NaOH extractable) accounts for~2% of total PP, with the majority extracted as NaOH-P (P bound to amorphous Fe; Table 3). Around 4% of extractable P is present as residual P. KS suspended sediments contain more total phosphorus than LG (~1200 μg g À1 ; Table 3). The higher P content is likely due to the catchment geology, which although similar to LG, includes basaltic intrusions. Basic rocks such as basalt generally have higher P concentrations (~1000 ppm) than shield rocks (~700 ppm) [Porder and Ramachandran, 2013]. The majority (>97%) of particulate phosphorus is in the HCl-P fraction, as with LG. Labile and residual phases contributed less than at LG, <1% and~2%, respectively.

Catchment Phosphorus Yields
We used catchment-normalized P yields [Föllmi et al., 2009] to compare the efficiency of P weathering in GrIS glacial catchments compared to nonglacial rivers. SRP yields from LG and KS were high and estimated at 26.8 (8.5-43.3) and 16.7 (5.6-33.3) kg km À2 yr À1 , respectively (Table 4). There is significant variation in catchment labile PP yields (Table 4). These differences (57.3 versus 8.3 kg km À2 yr À1 for LG and KS, respectively) are driven mainly by contrasting suspended sediment loads (Table 1), as extractable concentrations are of similar magnitude. Similarly, the LG TP yield (3250 kg km À2 yr À1 ) is nearly 4 times higher than the KS TP yield (859 kg km À2 yr À1 ; Table 4).

Flux Estimates
Fluxes of bioavailable P at LG change throughout the ablation season ( Figure 5). Periods of enhanced P export coincide with outburst events or the extreme melt event of 2012 (days 192-197), with the highest flux period occurring during the latter. This demonstrates the importance of increased meltwater input into the subglacial system as a mechanism of enhancing the P flux [Hawkings et al., 2015]. KS P fluxes are fairly consistent after the distinct hydrological change on day 157 ( Figure 5). Fluxes, ranging from 0.1 to 0.3 g s À1 , are nearly 2 orders of magnitude lower than LG, with dissolved fluxes more important than particulate fluxes.
We use our data to estimate annual P fluxes from LG, KS, and the GrIS ( Table 5). The mean annual dissolved flux from LG totaled 18.3 Mg yr À1 , with SRP the dominant form (~90%). Labile PP accounted for 34.4 (12.5-84.7) Mg yr À1 , approximately twice the SRP flux. The total P flux is estimated to be 1986 Mg yr À1 . KS fluxes are <10% those of LG, with estimated dissolved P export of 0.8 (0.2-2.5) Mg yr À1 , largely due to lower annual Global Biogeochemical Cycles 10.1002/2015GB005237 discharge. The KS labile PP flux is lower than the dissolved P (DP) flux, at 0.3 (0.1-0.5) Mg yr À1 (Table 5). Only a small fraction of total P export from KS is "labile," with most exported as more refractory phosphorus in sediments. Our flux estimates indicate that PP dominates the P flux from LG and KS (>97%). PP therefore dominates export from the GrIS; we estimate that 408 (238-1436) Gg P yr À1 is discharged into the surrounding coastal regions annually. The most highly reactive labile fractions are estimated to account for 10.9 (5.3-34.0) Gg P yr À1 , with dissolved fractions comprising~33%.

Catchment Comparison
There are clear differences in the P cycle at LG and KS. Catchment hydrology likely exerts a first-order influence on phosphorus concentrations since LG and KS bedrock is similarly composed of igneous crystalline rock geology. Concentration of bedrock P is elevated at KS (Table 3); thus, if catchment hydrology were similar, one would expect elevated concentrations at KS. Catchment size, and hence meltwater residence time, therefore appears to be an important factor in dictating SRP concentration where bedrock geology is relatively uniform [Wadham et al., 2010]. More remote regions farther from efficient drainage pathways are likely at the bed of the larger LG, allowing longer-term storage of subglacial waters. In comparison, the bed of the smaller KS is likely to be extensively flushed each year due to the closer proximity of waters to efficient drainage pathways. We therefore hypothesize that the higher concentrations of SRP at LG are the result of longer residence time waters, giving rise to enhanced biogeochemical weathering of bedrock. This hypothesis is likely to hold true for PP species as LG sediments have higher concentrations of labile phases.
Mechanisms of release also differ between LG and KS (section 3.1). Elevated P fluxes at LG are associated with subglacial outburst events (E1-E6) or extreme melt periods ( Figure 5) and are often sustained following such events. These hydrological events play a critical role as they enhance fluxes in the short term and drive changes in subglacial drainage by promoting subglacial water flow from higher up catchment. This likely allows drainage of increasingly isolated regions of the subglacial system, where waters high in P are likely to exist. Hence, this drainage system evolution promotes elevated P fluxes several days after outburst events.

10.1002/2015GB005237
It is probable that seasonal P trends will change at LG with climatic warming [Hawkings et al., 2015]. For example, the region of supraglacial lake formation and drainage on the GrIS is predicted to increase in a warmer climate [Leeson et al., 2015], potentially leading to a higher incidence of outburst events in the catchment. Furthermore, high P fluxes may initiate earlier, with elevated P fluxes continuing until later, as the melt season grows in length and intensity. Conversely, over time more regular flushing may lower bioavailable P concentrations, as the residence time of stored waters falls. There is large uncertainty in the response of suspended sediment flux to enhanced warming, as little data are currently available [Cowton et al., 2012;Hawkings et al., 2015;Hudson et al., 2014]. Changes to suspended sediment dynamics will have a large impact on fluxes of PP downstream. Evidence from past deglaciation events suggests that sediment fluxes may increase in the future [Jeandel and Oelkers, 2015].
An efficient, channelized drainage system is established fairly early in the season at KS. Little change is evident after approximately day 160, as reflected in KS temporal flux records ( Figure 5). Fewer isolated regions of subglacial drainage are likely to exist due to the smaller catchment size, and relatively rapid onset of efficient drainage in the ablation season, as found at Alpine glaciers [Swift et al., 2005]. Supraglacial lake-induced outburst events do not occur at KS.

Controls on Phosphorus Concentrations in GrIS Meltwaters
Most dissolved P in glacial meltwaters is derived from the dissolution of P-containing rock. Both LG and KS are located on crystalline bedrock with fluorapatite by far the most dominant mineral form of P [Porder and Ramachandran, 2013]. Physical erosion of bedrock by glacier action exposes finely ground apatite to biogeochemical weathering, liberating SRP. Hence, elevated concentrations of SRP are likely to arise from glacial comminution of bedrock.
Fe appears important in the glacial P cycle at LG (Figures 6 and 7), as has been observed in nonglacial river waters [Gunnars et al., 2002;Ruttenberg, 2014;Slomp, 2011]. High colloidal/nanoparticulate Fe concentrations (filterable through a <0.45 μm membrane; CNFe) were recorded at LG during this monitoring period . A significant correlation between SRP and CNFe up to CNFe concentration of 2 μM (R 2 = 0.6; P < 0.01; Figures 4 and 6) indicates sorption of SRP onto highly reactive Fe (oxy)hydroxide nanoparticles or colloids <0.45 μm in size [Froelich, 1988]. Changes in SRP and CNFe are also closely correlated to pH with concentrations highest at pH >9, where the solubility of iron (oxy)hydroxides increases (Figures 2a and 4) [Drever, 1997]. The iron (oxy)hydroxide point of zero charge is also at pH~9. When water reaches this pH the desorption of negatively charged SRP from the iron (oxy)hydroxide surface into solution is favored. This may also explain the weaker correction between SRP and CNFe when CNFe concentrations are >2 μM (Figure 6).
Buffering of SRP by particulate Fe fractions also appears to be important (Figures 4 and 7). There is a significant correlation between SRP and NaOH-P (Fe-bound P; R 2 = 0.6, P < 0.01; Figures 4 and 7). The MgCl 2 -P phase is likely a transition between dissolved and Fe-bound P (Figure 7) [Froelich, 1988;Ruttenberg, 2014].  Global Biogeochemical Cycles

10.1002/2015GB005237
High pH waters in the latter part of the sampling period (approximately day 170, Figure 2a), coupled to increasing SRP and MgCl 2 -P (Figure 7), suggest P is derived from NaOH-P. This may be linked to the higher solubility and point of zero charge of P-bearing Fe (oxy)hydroxides at pH >9 [Drever, 1997].
Conversely, we infer that at KS little cycling between phases is occurring. There is no significant relationship between Fe and P, and SRP and PP. This is likely to be due to the lower concentrations of highly reactive suspended sediments for weathering, and the comparatively short water residence times of small versus large glacier systems [Wadham et al., 2010]. We postulate that water residence times (in comparison with LG) are the critical control. This is consistent with observations of Fe-P relationships at other small glaciers [Föllmi et al., 2009;Hodson et al., 2004], attributed to the relatively short residence times of waters exiting the glaciers and limited secondary Fe mineral formation.
To our knowledge, this is the first time a temporal trend in more labile PP fractions has been shown (at LG; Figure 8a). NaOH-P is high in early season sediments and generally declines toward the end of the sampling period, while MgCl 2 -P does the opposite (Figure 8a). PCA for LG (Figure 4a) shows that MgCl 2 -P positively loads with SRP and day of year on component 2, while NaOH-P and Res-P are negatively loaded. NaOH-P and Res-P are also positively loaded on component 1, alongside major ions. There is a subtler trend at KS (Figure 4b), with significant values (taken as > 0.2 or < À0.2) on component 1 obtained for MgCl 2 -P (À0.379) and NaOH-P (0.478). These P species load with similar variables to LG. PP fractional concentrations are therefore closely tied to the geochemical conditions in glacial meltwaters. These temporal trends demonstrate the capacity for PP concentrations to change over relatively short timescales (days to weeks) in natural waters, such as seasonally evolving glacial meltwater runoff. a Estimates of atmospheric phosphorus fluxes and riverine phosphorus fluxes are also given. b Water and sediment flux are presented as km 3 yr À1 and Tg yr À1 , respectively. c Derived from minimum, discharge-weighted mean, and maximum concentrations as respective glaciers multiplied by glacier meltwater flux. d Derived from minimum, discharge-weighted mean, and maximum labile PP concentrations (MgCl 2 -P + NaOH-P) at respective glaciers multiplied by glacier suspended sediment flux. e Derived from minimum, discharge-weighted mean, and maximum concentrations at LG multiplied by GrIS meltwater flux . ) suspended sediment concentrations, and the discharge-weighted mean labile PP concentrations (MgCl 2 -P + NaOH-P) at LG (see Table 3; as in Hawkings et al. [2014] Holmes et al. [2012], assuming~50% of measured DP is SRP and~50% is DOP. Minimum and maximum DOP values are assumed to be~50% of DP [Ruttenberg, 2014]. j Using total P data from Guo et al. [2004] or Ruttenberg [2014] as indicated by "source." Due to lack of data, we assume that minimum, mean, and maximum labile P concentrations are 10% [Meng et al., 2015], 25%, and 40% [Ruttenberg, 2014] of total PP, respectively. k Due to lack of data it is assumed that~50% of DP is DOP [Ruttenberg, 2014].  Leverett Glacier meltwater SRP concentrations are similar to the global weighted mean nonglacial river concentration of 0.32 μM (Table 2) [Meybeck, 1982]. Concentrations of SRP in Arctic rivers range from 0.03 to 0.35 μM, with most <0.2 μM (Table 2). DOP can account for >50% of the total DP load in nonglacial riverine waters (Table 2) [Meybeck, 1982;Ruttenberg, 2014]. DOP also likely dominates dissolved fluxes from Arctic  watercourses. For example, mean concentrations of~0.9 μM and~2.3 μM reported by Meybeck [1982] and Emmerton et al. [2008] in Arctic rivers are much higher than glacial meltwaters. Higher DOP is consistent with elevated dissolved organic carbon concentrations observed in (sub)arctic rivers versus glacial meltwaters [Lawson et al., 2014;Lobbes et al., 2000]. However, the lability of DOP is currently not well understood, and significant enzymatic processing may be required before utilization [Cotner and Wetzel, 1992;Paytan and McLaughlin, 2007]. Hodson et al. [2004] presented NaOH-P (+ MgCl 2 -P) data from the suspended and proglacial sediments of five valley glaciers. Values were typically much lower than those found in this study, with the exception of Storglaciären (Table 3). Föllmi et al. [2009] published extraction data from temperate valley glacier sediments, but the method used was different to this study, and since the phases are operationally defined, both data sets are not directly comparable. As a result, the more aggressive extraction they used for Fe-/Al-bound P may explain elevated concentrations of 20-50 μg g À1 versus 10-15 μg g À1 in this study (Table 3).

Particulate Species
PP in glacial meltwaters appears to be significantly different in composition to PP from nonglacial rivers, although surprisingly little data are available for comparison. Riverine sediments potentially contain more labile P fractions. Specifically, organic PP comprises a much larger proportion of exported material (20-40%) [Meybeck, 1982;Ruttenberg, 2014]. Fe-/Al-bound P and mineral P (such as apatite) are generally responsible for the majority of inorganic PP (>40% of total PP) [Berner and Rao, 1994;Paytan and McLaughlin, 2007].
GrIS concentrations for Fe-bound PP phases tend to be significantly lower than those from riverine studies (Table 3). This is likely partly due to the extraction technique used. The NaOH extraction in this study is more selective for labile P (algal available) bound to the surface of Fe (oxy)hydroxides [Eijsink et al., 1997;Levy and Schlesinger, 1999], while the CBD (citrate-bicarbonate-dithionite solution) method [Ruttenberg, 1992] used in a number of riverine studies [e.g., Berner and Rao, 1994;Meng et al., 2015] typically gives higher estimates, because it removes more refractory Fe-P phases [Eijsink et al., 1997;Golterman, 1996]. We performed an additional extraction on fresh LG suspended sediments to test this. A well-characterized ascorbate extraction for highly reactive amorphous ferrihydrite was used [Raiswell et al., 2010], and 26.7 (± 8.9) μg g À1 (n = 32) of Fe-P was extracted versus 14.6 (± 7.5) μg g À1 for the NaOH-P extraction. This indicates that the NaOH-P extraction is selective for only the most labile Fe-P component. Our values are therefore conservative estimates compared to riverine Fe-P using the CBD method.
The organic PP content of glacial suspended sediments (Res-P) was low compared to riverine waters (Table 3). This is likely due to the low organic carbon content of LG sediments and the relative immaturity of freshly weathered subglacial debris compared to highly weathered riverine particulate matter [Lawson et al., 2014].
4.4. The Greenland Ice Sheet as a "Hot Spot" of Phosphorus Denudation 4.4.1. SRP Denudation Our results suggest that the GrIS is among the most efficient regions of P weathering and export in the world. Dissolved P yields are among the highest observed in literature (Table 2). These rates are elevated compared to those reported in other glacial catchments, apart from Hood and Scott [2008] (Table 4), where SRP yields >30 kg km À2 yr À1 were recorded at two North American glaciated watersheds.
SRP yields from LG and KS are larger than those of nonglacial riverine catchments. They are nearly 2 orders of magnitude larger than the modeled global mean SRP denudation rate of Harrison et al. [2010], an order of magnitude greater than the global mean estimate of Meybeck [1982], and are similar to those reported in some of the world's largest watersheds such as the Amazon, the Mississippi, and the Changjiang (Table 4) [Berner and Rao, 1994;Harrison et al., 2010;Sutula et al., 2004]. Our SRP yields are substantially higher than those from Arctic river catchments, which comprise the majority of freshwater delivery to the Arctic oceans (we estimate a discharge-weighted mean of 3.7 kg km À2 yr À1 ; Table 4) [Harrison et al., 2010]. Even if dissolved P yields were twice these estimates (to include organic P), GrIS catchment yields would still be significantly higher. These high observed glacial SRP denudation rates are likely due to bedrock comminution exposing finely ground P minerals for weathering, and the flushing of long-term stored waters from the ice sheet bed.

Particulate Phosphorus Denudation
From other studies on glacial catchments only Storglaciären has comparative labile PP yields to LG (Table 4) [Hodson et al., 2004]. The low yield estimated for cold-based Austre Broggerbreen (Table 4) [Hodson et al., 2004] Global Biogeochemical Cycles

10.1002/2015GB005237
is probably due to the absence of an active subglacial drainage system, and hence lower physical and biogeochemical weathering, demonstrating the importance of active subglacial hydrology. Differences in extraction technique mean that data from glaciers in the Swiss Alps are not directly comparable [Föllmi et al., 2009], although their values lie between those of LG and KS (Table 4).
Riverine labile PP yields are comparable to GrIS catchments (Table 4), despite labile particulate phosphorus concentrations in riverine particulates tending to be much higher (Table 3) [Ruttenberg, 2014]. This is due to the lower suspended sediment loads (and by association erosional force) of rivers compared to glacial meltwaters (<0.1 g L À1 compared to >1 g L À1 , respectively) [Meybeck, 1982]. We could only find measurements of PP concentrations in major Arctic rivers for the Yukon [Guo et al., 2004], and for a smaller Arctic river [Cai et al., 2008]. We combine Arctic river sediment yield data from Hay [1998] with discharge-weighted total PP data from Guo et al. [2004] (571 μg g À1 ), and mean labile PP phase data from the Amazon and Mississippi Rivers, where by far the most extensive PP fractionation records are available for riverine environments (~40% of total PP; Table 3) [Berner and Rao, 1994;Sutula et al., 2004], to derive estimates of Arctic river basin yields for comparison. We acknowledge that these estimates are crude, because they are based on values reported from temperate/tropical rivers and use a different Fe-P extraction to this study. Despite this, yields are similar to GrIS glacial catchments (Table 4). We estimate labile P yields of 1.4-21.3 kg km À2 yr À1 for the Arctic rivers, which are <37% of LG yields and similar to KS estimates.
A particularly important finding is that total P yields from LG and KS are substantially higher than all nonglacierized catchments in the literature (Table 4).
LG total P yields are approximately 50 times higher than Arctic river yields and the global average of 69.3 kg km À2 yr À1 [Beusen et al., 2005] and are at least an order of magnitude greater than estimates from rivers with elevated physical weathering and particulate loads, such as the Ganges, the Changjiang, and the Irrawaddy (Table 4).
LG total P denudation rates also exceed estimates from other glacial catchments by~60-1700%, while KS total P yields are similar. Thus, high rates of physical erosion and biogeochemical weathering of sediments in the subglacial environment, coupled with large annual freshwater discharge, indicate that the GrIS is a hot spot of P weathering and export in the Arctic.

Importance of Greenland Ice Sheet P Fluxes to the Arctic Oceans
We generated a labile P Arctic river flux estimate of 70.4 (20.4-153.6) Gg P yr À1 , using the methods detailed in Table 5. Atmospheric P deposition for the Arctic is taken from Mahowald et al. [2008]. They estimate annual Arctic Ocean atmospheric TP deposition of 6.9 Gg yr À1 , with 2.2 Gg solubilized to SRP (the labile fraction; Table 5). These combined P fluxes total 72.6 Gg P yr À1 ; thus, we estimate meltwater runoff from the GrIS constitutes approximately 15% (10.9 Gg yr À1 ) of the labile P flux to the Arctic oceans. More significantly, we estimate that TP fluxes from the ice sheet (~408 Gg yr À1 ) are more than 3 times that of the flux from Arctic rivers (~125 Gg yr À1 ) and approximately 60 times larger than atmospheric TP deposition (~7 Gg yr À1 ). GrIS fluxes are similar to those of the Chianjiang, Amazon, and Mississippi Rivers (Table 5).
It has generally been observed that a fraction of sediment P is rapidly released to solution in artificial estuarine waters [Froelich, 1988]. This explains the characteristic nonconservative behavior of SRP through salinity gradients. The mineral fraction (believed to be apatite) that dominates TP export from the GrIS is not usually considered labile. However, some of the less refractory PP is likely to be solubilized in the proglacial region, and in suspended and deposited sediments in near coastal regions [Oelkers et al., 2008], as glacial suspended sediments have been observed to be extremely fine and highly reactive [Brown et al., 1996]. Lab experimentation suggests finely ground apatite dissolves in less than a year in seawater at 4°C [Chairat et al., 2007;Jeandel and Oelkers, 2015]. Thus, we believe that the high fluxes from the GrIS could have a large effect on P availability in downstream marine ecosystems after processing.
A large majority of sediment P does not reach the open ocean in riverine systems and therefore does not play an active role in short-term biogeochemical processes. Nixon et al. [1996] estimate that 75% of TP from large rivers would be deposited and buried in deltas, assuming 25% is solubilized [Froelich, 1988]. It is unclear if this holds true for the GrIS, where large deltas do not exist but large fjord systems do, as there is currently no data. We acknowledge that fluxes will also depend on other physical oceanographic factors around the GrIS, which might not favor significant off-shelf export [Hopwood et al., 2015]. This could significantly lower the flux of P into the open ocean, as with nonglacial regions.

10.1002/2015GB005237
Our estimates do not include the P flux from icebergs. Icebergs constitute~60% of the GrIS freshwater flux [Bamber et al., 2012], rafting fresh glacial ice melt, and entrained debris off-shelf into the open ocean [Smith et al., 2013]. If we assume that iceberg-entrained debris has a similar labile PP extractable content as meltwater runoff (~10 μg g À1 ) and that the mean sediment content of icebergs is~0.5 g L À1 [Raiswell and Canfield, 2012;Smith et al., 2013], we calculate a crude flux of 3.1 Gg yr À1 . This value is approximately the same as the atmospheric flux to the Arctic oceans and~5% of the estimated labile flux from Arctic rivers.
We use discharge data from 2012 as a crude indicator of how fluxes may change in a warmer climate (Table 5). These estimates indicate labile fluxes of P were enhanced by~50% compared to the 2000-2012 mean, assuming nutrient concentrations and suspended sediment loads are roughly scalable to meltwater flux, as has been indicated in other studies [Hawkings et al., 2015;Hudson et al., 2014], although large uncertainties exist. The predicted increase in glacial meltwater discharge of >100% by 2100  will likely outstrip the rate of increase in riverine discharge to the Arctic oceans, which is estimated to be 18-70% during the same time period [Peterson et al., 2002]. We therefore postulate that GrIS meltwater runoff may become more important to the Arctic P cycle in the future.

Importance of Findings in a Wider Context
During recent glacial maxima ice sheets have been present over much of North America and northern Europe. These ice masses discharged significant quantities of glacial meltwater and associated sediments into the oceans during the glacial period, and during initial deglaciation [Clark et al., 2009]. We also hypothesize that this delivered large quantities of bioavailable phosphorus to downstream environments [Föllmi et al., 2009]. The snowball Earth glaciations during the Cryogenian (720 to 635 Ma) are believed to be important in elevating phosphorus concentrations and productivity in the Precambrian oceans [Planavsky et al., 2010]. This may have subsequently driven the Neoproterozoic oxygenation, and evolution of metazoan life [Planavsky et al., 2010;Sahoo et al., 2012]. Our results indicate that widespread glacial weathering could have supplied large quantities of bioavailable P to the ocean during these periods.

Summary and Conclusions
The Arctic Ocean exports large quantities of nutrients into the Pacific and Atlantic Oceans, but the sources of phosphorus are currently unclear [Torres-Valdés et al., 2013]. We argue that glaciers efficiently transport large quantities of phosphorus to the surrounding coastal regions. High rates of physical erosion, combined with an abundance of meltwater at the bed of the ice sheet, drive biogeochemical weathering of phosphorus minerals. We find dissolved phosphorus yields are at least equal to some of the world's largest rivers, such as the Mississippi and the Amazon. Labile fluxes from GrIS meltwaters may constitute as much as 15% of the total Arctic riverine flux (~11 Gg P yr À1 ), which has previously been considered the main source of nutrients for the surrounding ocean [Torres-Valdés et al., 2013], and are significantly larger than the flux from atmospheric sources [Mahowald et al., 2008]. Leverett Glacier, a large catchment of the GrIS, was found to have total phosphorus yields of over 3000 kg P km À2 yr À1 , which is over an order of magnitude higher than yields from riverine catchments. The ice sheet total phosphorus flux of~400 Gg yr À1 rivals that of the largest river systems in world. Although there are large uncertainties, we propose that fluxes are annually dynamic [Hawkings et al., 2015], increasing with meltwater export under warming climate scenarios. The significance of GrIS P fluxes depends on the deposition and burial of sediments in the fjord and coastal regions. This is a large unknown as no data currently exist and should be a priority of future research. Our findings have implications for our understanding of nutrient cycling during glacial transition periods and during past snowball Earth glaciations. ship to J.H. A.T. was funded by a NERC studentship and MOSS scholarship. P.N. was supported by grants from the Carnegie Trust for University of Scotland and the University of Edinburgh Development Trust. Additional support was provided by the Leverhulme Trust, via a Leverhulme research fellowship to J.L.W. We thank all of those who assisted with fieldwork at LG and KS, as well as technical staff in LOWTEX laboratories at the University of Bristol. The data used in this article are available from the corresponding author (jon.hawkings@bristol.ac.uk) upon request. We are grateful to our anonymous reviewers for their constructive comments on the manuscript.