An improved estimate of microbially mediated carbon fluxes from the Greenland ice sheet

Abstract

Microbially mediated carbon fluxes on the surface of the Greenland ice sheet (GrIS) were recently quantified by Hodson and others (2010) using measurements of the surface coverage of debris (cryoconite) and rates of biological production associated with debris near the ice-sheet margin. We present updated models that do not assume the same spatial uniformity in key parameters employed by Hodson and others (2010) because they make use of biomass distribution and biological production data from a 79 km transect of the GrIS. Further, the models presented here also include for the first time biomass associated with both cryoconite holes and surficial algae. The predicted annual carbon flux for a small (1600 km²) section of ice surrounding the field transect is about four times that estimated using spatially uniform biomass and production in this area. When surficial algae are included, the model predicts about 11 times more carbon fixation via photosynthesis per year than the cryoconite-only models. We therefore suggest that supraglacial carbon fluxes from the GrIS have previously been underestimated by more than an order of magnitude and that the hitherto overlooked surficial algal ecosystem can be the most crucial contributor. The GrIS is shown to be in a relatively stable state of net autotrophy according to our model and so a strong link between bare-ice area and total carbon fluxes is evident. The implication is a biomass feedback to surface albedo and enhanced ablation as a result. Climate predictions for the year 2100 show that greater carbon fixation could also result from climate warming.

Introduction

Biologically active sediment (cryoconite) in cryoconite holes on glacier and ice-sheet surfaces has been recognized as a significant contributor to regional carbon cycling (Anesio and others, 2009;Hodson and others, 2010). Microbial communities in cryoconite holes both respire and photo- synthesize at rates comparable to microbial ecosystems in much warmer nutrient-rich environments (Hodson and others, 2007, 2010;Anesio and others, 2009, 2010; Telling and others, 2011). Anesio and others (2009) predicted the carbon flux for all glaciated areas outside Antarctica and found that net carbon fixation of ~64Gga^-1 is likely. Their analysis assumed constant rates of biological production (photosynthesis and respiration) that were derived from the averaging of measurements conducted at a small selection of sites in Svalbard, Greenland and the European Alps. A constant ablation area was also assumed and so a conservative estimate of the duration of biological production during summer (the ’growing season’) was then employed to avoid overestimating total carbon fixation by inclusion of periods of snow cover (when photosynthesis is prevented by shading). Hodson and others (2010) then estimated the carbon flux from the Greenland ice sheet (GrIS) only, but allowed spatial and temporal evolution of the ablation zone using an adapted version of the runoff retention model of Janssen and Huybrechts (2000). This was used to simulate the extents of slush and bare ice at a monthly time-step. Photosynthesis was allowed to occur in the bare-ice zone, but not in the slush zone (where respiration was assumed to be the only process occurring). Within the bare-ice zone the surface coverage by cryoconite, mass per unit area of cryoconite, and rates of primary production and community respiration were all assumed to be spatially uniform, using field measurements from a short transect (0–5 km from the ice margin, near Kangerlussuaq). For the entire ice sheet, carbon fluxes of the order of ~10¹-10² GgCa^-1 due to photosynthesis and respiration were extrapolated. However, photosynthesis measurements were obtained late in the ablation season, hence rates of photosynthesis were low (~20 GgCa^-1) and assumed to represent a low end-member scenario. Global average rates of photosynthesis taken from previous work (Anesio and others, 2009) were therefore incorporated into the model to provide a high end-member scenario (~180 GgCa^-1), resulting in a carbon balance that was dominated by net carbon fixation at the surface because in each case respiration only amounted to 25 GgCa^-1.

In the present study, a field campaign was undertaken in August 2010 to collect field data that better describe cryoconite coverage, mass of cryoconite per unit area and productivity along a linear transect crossing the entire ablation zone. During this field campaign, a near-ubiquitous surface algal community was also observed between cryoconite holes (where present) and was therefore included in the sampling regime. The aims of this paper are to improve our understanding of the carbon fluxes that are associated with supraglacial microbial production upon the GrIS. We achieve this by modelling the duration and magnitude of biological activity as carbon fluxes associated with: (1) spatially uniform cryoconite distribution (as employed in previous research);(2) spatially variable cryoconite coverage and biological production, as sampled along our transect;and (3) spatially variable cryoconite and surface algal coverage and biological production, defined using our field observations.

Methodology

Field Methods

Fieldwork was undertaken between 22 July and 10 August 2010 along a 79 km transect of the GrIS beginning at the terminus of Leverett Glacier (67°04′17.1" N, 50°08′45.2″ W to 67°09′10.8″ N, 48°22′14.6″ W) (Fig. 1). The change in elevation along the transect was 1047 m (399 m a.s.l. at 0 km from the ice margin to 1446 m a.s.l. at 79 km from the ice margin). On 2 and 3 August 2010, a helicopter was used to access field sites at 7.5, 17, 34, 51 and 79 km along the transect. Field sites at 0, 2, 4, 5.7 and 7.5 km were accessed on foot on 5 and 6 August 2010. The 79 km site was covered by fresh snow, so the measurements taken there were not included in this study. Carrying out this fieldwork in August provided a bare-ice zone which approximated the year’s maximum for that area and diurnal light cycles at surface irradiance levels of ~1800mm photonsm^-1 s^-1, which were considered sufficient to maintain high rates of primary production. Until more measurements are done we cannot know if the rates we measured were representative of those normally attained. At each site a series of five 1m × 1m quadrats were randomly sampled, within which the following data were collected.

Fig. 1. Map showing the locations of the sampling points along the 79 km transect. The inset depicts the location of the transect within Greenland. Adapted from Telling and others (2011).

Percentage coverage of the ice surface by cryoconite

Nine 10cm × 10cm ‘calibration quadrats’ were randomly emplaced at each site on the transect. The entire mass of wet cryoconite contained within each calibration quadrat was removed using a 100mL syringe and injected into a measuring cylinder. The wet volume of cryoconite was recorded. The cryoconite was sealed in 50 mL sample containers and returned to the laboratory at Kangerlussuaq International Science Station (KISS) within 8 hours of collection and dried in an oven for 24 hours at 40°C (after Hodson and others, 2010) to record dry mass. Photographs of each calibration quadrat were analysed in Adobe Photoshop. The images were cropped and resized to standardize pixel numbers in the 10 cm × 10cm area of interest. The ‘select colour’ tool was used to isolate pixels that represented cryoconite. This was based upon the assumption that dark pixels represented cryoconite, whereas light pixels represented bare ice;however, some chromatic adjustment was necessary to prevent shadows and reflections from being selected. Typically a high red balance aided automatic selection of pixels representing cryoconite. Additionally, some images required manual removal of reflections of cryoconite sediment, which was achieved using the ‘clone stamp’ tool. The total number of selected pixels was displayed in Photoshop’s ‘Info Palette’ and provided a measure of percentage areal coverage. The relationship between areal coverage and dry mass for the 10 cm × 10 cm calibration quadrats is shown in Figure 2. This relationship was upscaled to provide a measure of areal coverage and dry mass from photographs of the 1 m × 1 m quadrats, which were emplaced at each site along the transect.

Fig. 2. Relationship between areal coverage and dry mass for calibration quadrats. Note that calibration quadrats refer to 10cm x 10 cm quadrats emplaced at the 2 km site for which the mass of cryoconite present was measured, not the 81 test quadrats emplaced along the transect.

Surface algae were too small and their distribution too extensive for the percentage coverage and mass per unit area to be estimated using the image analysis method. Instead, and due to significant time restraints with the helicopter, each quadrat was qualitatively assigned to one of three algal coverage categories: light, medium or heavy. Quadrats which had ‘light’ coverage were assumed to have 5% of their surface covered by algae, ‘medium’ quadrats were associated with 40% surface cover and ‘heavy’ quadrats were associated with 90% coverage, as determined from manual pixel counts from quadrat images. In order to model the percentage surface cover, the number of quadrats in each category at each site was counted and a weighted average of the coverage estimates was calculated for each site along the transect.

For quadrats in each coverage category, a 10 cm layer of ice was removed by carefully chipping away the ice surface using an ice axe, placed in a 1 L Whirl-Pak® bag and immediately returned to the laboratory at KISS. Here chlorophyll-a assays were undertaken for each sample within 6 hours of collection, from which the mass per unit area of algae was determined. Chlorophyll-a counts were assumed to represent 5% of the biovolume of autotrophic microorganisms after Nicholls and Dillon (1978) who showed that chlorophyll-a measurements represented 0.19.7% of autotrophic biovolume in laboratory experiments. Biovolume was then converted to algal biomass by assuming a specific gravity of 1 gcm^-3 after Nauwerck (1963). Average mass per unit area was determined for each coverage category, and then weighted means calculated for each site.

Primary production and respiration measurements

Net ecosystem production (NEP) and respiration (R) rates for samples of cryoconite sediment were measured using the ΔTDIC (change in total dissolved inorganic carbon) method developed by Hodson and others (2010) and described further in the context of the present study by Telling and others (2011). Since we were motivated by an interest in the biologically mediated carbon fluxes across the air-ice interface, we define NEP as R –PP (PP represents photosynthesis), so that positive values indicate potential CO₂ transfer to the atmosphere (net respiration) and negative values denote net carbon fixation by the ice surface (net autotrophy). From two large cryoconite holes at the 2 km site, cryoconite and cryoconite hole water were added to 60 mL glass BOD (Biological Oxygen Demand) bottles at a ratio of ~1:60. The BOD bottles were naturally illuminated apart from those used for dark incubations, which were completely covered in aluminium foil. All BOD bottles were submerged in cryoconite holes. Incubations were then carried out in situ for 24 ± 1 hours. Gross photosynthesis (PP) rates were calculated as NEP-R. The rates were normalized for the dry weights of sediment in the bottles, determined by drying and weighing the sediment. The detection limits were 1.2 mgCg^-1 d^-1 for NEP and R, and 1.7 mgCg^-1d^-1 for PP. Respiration rates were assumed to be equal for bare ice and snow-covered ice, since respiration is not light-limited, and were measured using the dark incubations. Surface algal communities were found to exist in a thin layer on the ice surface rather than occupying an aqueous habitat like a cryoconite hole. Further, time constraints meant that algae could not be incubated in suspension over the course of 24 hours using the cryoconite ΔTDIC method. Rapid (<1 hour) incubations were therefore employed using¹⁴C label incorporation after melting ice and algae samples (Telling and others, 2011). Owing to the very short incubation time, productivity measurements derived from label incorporation represent in situ rates of gross photosynthesis and we assume that the effects of incubating a suspension (rather than a ‘dry’ algal layer) are minimal. A shortcoming of this approach is that no algal respiration data could be obtained. A dark chamber could have been used to measure algal respiration; however, time restraints were prohibitive. Ideally, cryoconite would also have been incubated using label incorporation, but only a limited volume of isotope was available. Estimates of photosynthesis derived using label incorporation and total carbon were recently compared by Telling and others (2011). However, label incorporation measurements were shown to underestimate rates of photosynthesis by more than an order of magnitude compared with the total carbon method and were close (although not identical) to rates of NEP measured using total carbon. This is because the duration of the incubations was 1 day, allowing time for autotrophic respiration of the label. Therefore the 1 hour incubations were used in the present study. At each site on the transect, a minimum of four replicates for the light incubations and two replicates for the dark incubations were used. At each location, the BOD bottles were all emplaced in a large cryoconite hole following sediment addition to maintain natural illumination and temperature conditions. These incubations were performed between 2 and 3 August 2010.

Cryoconite hole frequency, grain size and sediment arrangement

The average grain diameter for cryoconite granules in cryoconite holes was estimated for each quadrat at each site by removing samples with a 100 mL syringe and spreading the granules across laminated graph paper. The graph paper was marked with 1 mm squares, allowing rapid estimates of grain size to within 0.5 mm for a large number of grains. This method is subject to the least error furthest away from the ice margin, where grain sizes were larger. Close to the ice margin <1 mm grains were observed, for which the error associated with the graph paper method was most likely in excess of ~50%. The number of individual cryoconite holes occurring in each quadrat was also counted and the presence or absence of multiple grain layers noted.

Modelling methods

Hodson and others’ (2010) model of carbon fluxes on the GrIS was based upon three equations:

(1)

(2)

(3)

CR^snow represents the carbon flux from respiration beneath snow or slush; CR^ice represents the carbon flux from respiration on bare ice; CP^ice represents the carbon flux from photosynthesis on bare ice; A bio represents the area of the ice sheet covered in biologically active cryoconite; /is the proportion of A_bio that is not snow-covered;M _d is the mass per unit area of the cryoconite (gcm^-2);and c is the percentage cover of cryoconite within A_bio. R^snow and R^ice are respiration rates before and after snow removal, and PP^ice is the rate of photosynthesis after snow-cover removal (all in mMCg^-1 d^-1). The term d_m represents the time-step and k is a constant used for unit conversion. There is no CP^snow variable because it is assumed that snow cover inhibits photosynthesis due to shielding of cryoconite sediment from solar irradiance.

Hodson and others (2010) used an adapted monthly version of Janssens and Huybrechts’ (2000) runoff retention model, which was run at a resolution of 5 km × 5 km using downscaled European Centre for Medium-Range Weather Forecasts (ECMWF) operational analysis meteorological data (Hanna and others, 2011) in order to grow A_bio and i for the period 2000–08. The runoff model is described in Hanna and others (2005) and has been used for several ice-sheet studies (e.g. Hanna and others, 2008;Sundal and others, 2009, 2011;Murray and others, 2010). In the carbon flux model of Hodson and others (2010), fixed values were used for M_d (0.2 g cm^-2), c (3%), R^snow, R^ice (both 1.74 μmCg^-1 d^-1) and PP^ice (1.56 μmCg^-1d^-1), as derived from average values measured in Kronprins Christian Land, Kangerlussuaq and Thule ramp (all in Greenland). Owing to fieldwork being undertaken late in summer, values for PP^ice were likely to have been anomalously low due to seasonally declining solar intensity. To account for this, Hodson and others (2010) used a second scenario in which a PP^ice value of 14.7 μm Cg^-1 d^-1 was used: the average rate of photosynthesis from non-lidded cryoconite holes for a selection of Svalbard, Greenland and Austrian glaciers from Anesio and others (2009). The field measurements of Hodson and others (2010) yielded a net autotrophic carbon flux of 5.3GgCa^-1, and the second scenario yielded 152 Gg C a^-1 of net carbon fixation. There is no evidence to suggest the rates of autotrophic production suggested by Anesio and others (2009) are representative, and rates of PP used by Hodson and others (2010) have been identified as anomalously low. High uncertainty therefore permeates these carbon flux estimates.

Model 1

Model 1 represents a repeat of the model of Hodson and others (2010) with two additional years of meteorological data (2009 and 2010) driving the adapted runoff retention model of Janssen and Huybrechts (2000), and updated parameter values. Fixed values for all parameters were used; however, these values were derived from a 79 km transect covering the entire ablation zone, instead of the original 5 km transect used by Hodson and others (2010). Field measurements made at 0 and 79 km from the ice margin were disregarded in this model run, since at the margin the surface coverage was 100% and composed largely of sediment of subglacial origin (low photosynthetic activity), while at 79 km fresh snowfall obscured any surface impurities. Model 1 represents a single simulation, in contrast to the maximum and minimum scenarios modelled by Hodson and others (2010), since the new field values are more representative of true averages than the late-season measurements of Hodson and others (2010). The new M_d value was 1.02 g cm^-2. This value is higher than other Kangerlussuaq measurements (e.g. Cook and others (2010) found 0.21 gcm^-2), likely due to previous measurements being limited to locations <2 km from the ice margin, where coverage by cryoconite is demonstrably lower and grain sizes smaller than elsewhere on the ice sheet. The percentage surface cover by cryoconite (c) was 2.61%, slightly lower than Hodson and others’ (2010) estimate of 3%. This estimate excludes the marginal zone, which would have significantly raised the average percentage surface cover. Rates of respiration under snow (R^snow) and on bare ice (R^ice) were assumed to be equal by Hodson and others (2010) because respiration is not light- limited; however, respiration by photoautotrophs, changes to habitat hydrochemistry and carbohydrate availability in the presence of active autotrophs may complicate the situation. Unfortunately, data are still not available to elucidate this, and the present paper also assumes R^snow=R^ice.

Model 2

Model 2 comprised a repeat of model 1 applied to a 1600km² section of the GrIS: a small area for which the assumption of spatial uniformity is better justified. The coordinates of the field transect were used to delineate a rectangular section of ice sheet 20 km wide and 80 km long, which the transect line bisects. The transect was almost parallel to lines of latitude and to the direction of snow retreat. This allowed a rectangular section of ice sheet to be selected in which uniform snow retreat could be assumed, and the calculation of bare-ice area resulting from a given distance of retreat was simple. In order to force the transect to be precisely horizontal the latitudes were adjusted by ~1°, resulting in a major increase in computational efficiency.

Model 3

In model 3, the fixed values for R^snow, R^ice, pp^ice, M_d and c used in models 1 and 2 were replaced with variable estimates collected along the transect, which were applied to the 1600km² section of ice sheet employed by model 2. Only carbon fluxes associated with cryoconite were included. For each parameter, field-measured values were plotted against distance from the ice margin, and linear interpolation (point-to-point) was used to estimate parameter values between each pair of points (Fig. 3). The interpolated relationships were subsequently used to update a cumulative average of the parameter values for the entire bare-ice zone as it expanded and contracted seasonally. An approximation of the spatial variability of the model parameters is therefore included in model 3. The position of the snowline was determined using melt area output from the modified runoff retention model of Janssen and Huybrechts (2000), adapted to provide the total bare-ice area for the 1600 km² section of the GrIS at Leverett Glacier. A marginal zone was also included, which was assumed to be 100 m wide and with constant R^snow, R^ice, pp^ice, M_d and c, as measured on the Leverett Glacier transect. The resultant model equations are

(4)

(5)

(6)

Fig. 3. (a) Percentage of ice surface covered by cryoconite with distance from ice margin (linear fit shown); (b) mass of cryoconite per square centimetre with distance from ice margin (linear fit shown); (c) respiration rates for cryoconite with distance from ice margin (third-order polynomial fit shown); (d) rate of primary production for cryoconite with distance from ice margin (no error bars because PP was calculated directly from NEP and R; third-order polynomial fit shown); (e) percentage cover of ice surface by surficial algae with distance from ice margin (logarithmic fit shown); (f) mass of algae per square metre with distance from ice margin (logarithmic fit shown); (g) primary production of surface algae against distance from ice margin (third-order polynomial fit shown); (h) respiration rate for surface algae against distance from ice margin (no error bars as this was assumed to be a constant proportion of algal PP; third-order polynomial fit shown)

Model 4

Model 4 included carbon flux contributions from surficial algal communities in addition to cryoconite. Gross photosynthesis measurements deduced from the short-term radiolabel incorporation incubations were used and it was assumed that respiration and bacterial growth in the community could be accurately estimated as a fixed proportion of gross photosynthesis, according to Eqns (7) and (8). R^alg refers to the respiration rate in the algal community, which includes algal and bacterial respiration. Since data were available for bacterial production but not bacterial respiration, bacterial respiration was calculated using an assumed bacterial growth efficiency of 10% (after Anesio and others, 2010) and was therefore equal to nine times the bacterial growth rate. It is also assumed that 50% of the carbon fixed by algal photosynthesis is used in algal respiration (after Williams and Lefevre, 2008) which is therefore defined as 0.5 times the rate of gross primary production. Owing to these uncertainties, the model calculations provide only a first-order estimate of the potential carbon fluxes of ice algae:

(7)

(8)

where BR = 9 × bacterial growth rate and AR = 0.5 (gross primary production rate).

The existence of the algae in a thin film influenced the model formulae, since mass per unit area was negligible. The productivity of the algae was therefore directly measured per square metre of algal coverage, and the term M_d was omitted from the model equations. However, a measure of percentage cover (c) of the ice surface by algae was included in the model calculations, which accounts for the different spatial density of algal colonization observed in different regions on the GrIS. In addition, rates of photosynthesis and respiration in cryoconite were normalized to surface area instead of mass per unit area in model 4 in order to facilitate direct comparison between cryoconite and algal carbon fluxes. Further, the surface area of a given mass of sediment can vary and is closely coupled with rates of photosynthesis (Cook and others, 2010), suggesting normalization of rates of primary production for sediment mass may be inappropriate. The model equations for algal carbon flux were therefore

(9)

(10)

(11)

In model 4, the carbon flux for both cryoconite and surficial algal ecosystems was calculated for the section of ice sheet described in models 2 and 3. Production in a 100m marginal zone was also included. The resulting model equations were

(12)

(13)

(14)

Model summary

In order to elucidate for the reader the differences between each of the models we employed, Table 1 provides a summary of the properties of each model.

Table 1. Summary of the properties of models 1–4

Climate scenarios

The carbon flux from the GrIS was simulated for present-day climatic conditions using model 4. Two additional hypothetical climate scenarios were also modelled in order to provide the first approximation of the response of supra- glacial carbon fluxes to predicted future climate change for the period 2100–10. The current literature includes some predictions of climate change in coastal regions of the GrIS. For example, Huybrechts and others (2004) predicted a warming over the GrIS of +4.8°C and +4.1 °C from two separate models, with associated precipitation increases of +42% and +21%. The present study therefore used the average temperature and precipitation changes from the models of Huybrecht and others (2004) to predict a climate warming of 4.5°C and a precipitation increase of 32% by the year 2100. This was interpolated for the years 2100–10 to provide input data for model 4, in a simulation hereby referred to as the ‘warmer GrIS scenario’.

Hanna and Cappelen (2003) showed that coastal areas of the GrIS can cool concurrently with general climate warming, suggesting a –1.3°C temperature change between 1958 and 2001 relating to North Atlantic Oscillation changes, although this trend was subsequently reversed (Hanna and others, 2008). Bromwich and others (1993) also predicted a ~15% decrease in precipitation between 1963 and 1988, and cooling between 1940 and 1991 was simulated from reanalysis data by Smith (1999). Extrapolating the data from Hanna and Cappelen (2003) predicted a temperature change of a further –3°C and 15% decrease in precipitation in 2100, which was used to drive a ‘cooler GrIS scenario’. This represents a sensitivity experiment for model 4 rather than a realistic simulation of future carbon fluxes, since the bulk of the literature suggests that the GrIS is unlikely to cool during the 21st century. The carbon flux in both 2100–10 predictions assumes that temperature and runoff changes as a result of climate change do not influence the coverage or productivity of the organisms inhabiting the supraglacial environment, a limitation which should be considered in future studies.

Results

Field results

Figure 3 shows strong spatial variability in all measured parameters along the study transect. Maximum coverage, mass per unit area, respiration rate and rate of primary production for both cryoconite and surface algae occurred beyond 5 km and all parameters in Figure 3 exhibited a positive correlation with distance (note that primary production was measured as a reduction in carbon, therefore more negative values represent higher rates of photosynthesis). Cryoconite hole frequency was found to peak at 6 km, although an overall negative correlation (p <0.05) was evident over the entire transect (data not shown). A single grain layer arrangement was observed in every cryoconite hole, while the grain size increased with distance from the ice margin (data not shown). This is indicative of an increase in mass per unit area with distance from the ice margin. There is a positive linear relationship (r ² = 0.89;p <0.05) between the percentage of the surface covered by cryo- conite and distance from the ice margin (Fig. 3), despite the cryoconite hole frequency decreasing along the transect. This suggests that the surface area of each cryoconite hole increased with distance from the ice margin, which was confirmed by measuring cryoconite hole diameters. The mass per unit area was much higher for cryoconite than surface algae at all sites, by ~3–4 orders of magnitude. However, the total coverage and productivity rates were much higher for surface algae.

Model results

Table 2 shows that model 1 predicted a net carbon flux approximately twice that predicted by Hodson and others (2010). Figure 4a provides a comparison between the rates of primary production predicted by model 1 and Hodson and others (2010). The overall positive trend evident in model 1’s output is not apparent in the data from the model of Hodson and others (2010), which also does not appear to resolve the interannual variability predicted by model 1. The respiration datasets (Fig. 4b and c) also indicate a failure to resolve interannual variability by the model of Hodson and others (2010).

Fig. 4. (a) CP^ice, (b) CR^snow and (c) CR^ice as predicted by model 1 and Hodson and others (2010). Negative fluxes indicate carbon fixation by autotrophs.

Table 2. NEP as predicted by each of the carbon flux models

The outputs from models 2 and 3 constrain carbon flux estimates to the 1600 km² section of ice sheet. Table 2 shows that although model 2 includes field measurements from a greater range of locations than previous models, a significantly more negative carbon flux (representing higher rates of photosynthesis) is predicted by model 3. The mean annual carbon flux estimate was ~30 times higher from model 3 than from model 2 (-0.000069 GgC a^-1 for model 2 and –0.0022 GgC a^-1 for model 3).

Table 2 and Figure 5 show significant interannual variability in the NEP predicted by model 3, which is not resolved by the model of Hodson and others (2010). In 2001, model 3 predicted a particularly strongly net autotrophic NEP;however, in the same year the model of Hodson and others (2010) predicted net heterotrophy. An unexpected feature of the model was that between 2000 and 2006 a shift in NEP from net autotrophy towards balance was observed. Interannual variability in NEP was comparable between the two models except between 2001 and 2003, during which period a maximum discrepancy of >4 kg C km^-2 is observed.

Fig. 5. NEPas predicted by model 3 and by Hodson and others (2010).

Figure 6a and b depict the carbon flux resulting from autotrophy and community respiration on bare ice as predicted by model 3 and Hodson and others (2010) in kgCkm^-2. In both, the model employed by Hodson and others (2010) estimates greater carbon fluxes per square kilometre. For example the carbon flux resulting from autotrophy in 2003 predicted by Hodson and others (2010) is about three times greater than predicted by model 3. The only significant convergence of the two models occurs in 2002, when the CP^ice values are approximately the same. Similarly, the carbon flux resulting from cryoconite community respiration in 2003 is about seven times greater when the model of Hodson and others (2010) is used compared with model 3. The greatest discrepancy between NEP values predicted by the two models occurred in 2001, which coincides with the greatest bare-ice extent. There is a general trend (r ² = 0.65) of convergence between the two model predictions as the bare-ice extent decreases.

Fig. 6. (a) Carbon flux resulting from cryoconite community autotrophy as predicted by model 3 and Hodson and others (2010); (b) carbon flux resulting from cryoconite community respiration as predicted by model 3 and Hodson and others (2010).

Model 4 estimates that the 1600 km² section of ice sheet (representing about 0.007% of the total ice-sheet area) actually fixes ~16.4 GgCa^-1, whereas the entire ice sheet was estimated to fix ~5GgCa^-1 by Hodson and others (2010) (Table 2). This was the result of ~19.5GgCa^-1 respired by heterotrophs and ~35.9GgCa^-1 fixed by autotrophs. The majority (~92%) of this carbon fixation was concentrated in algal communities which photosynthe- size ~33GgCa^-1, outperforming cryoconite by ~11 times. Similarly, the heterotrophic consumption by cryoconite was 2.1GgCa^-1, whereas the estimated algal community respired 17.4GgCa^-1. The maximum rates were observed in 2002 when the NEP reached 37.3 GgC of carbon fixation, comprising 68.7GgC autotrophic carbon fixation and 31.4GgC heterotrophic carbon liberation. This coincided with the maximum modelled melt rates and bare-ice zone extent. Marginal sediment contributed only 0.001% (0.00017 GgC a^-1) of the total NEP.

Yearly fluctuations in NEP for cryoconite and algae in Figure 7 correlate strongly (correlation coefficient 0.97). There is a negative correlation (-0.88) between the total NEP and the bare-ice area. Correlations between NEP and algal autotrophic carbon flux (0.95), algal heterotrophic carbon flux (0.98), cryoconite autotrophic carbon flux (0.77) and cryoconite heterotrophic carbon flux (0.82) were found, all significant at p <0.05.

Fig. 7. Carbon flux contributions from algae and cryoconite between 2000 and 2010 predicted using model 4.

Climate scenarios

The bare-ice areas and carbon fluxes resulting from the warm and cool GrIS scenarios are presented in Figure 8, along with the 2000–10 areas for comparison. The warmer climate scenario predicted bare-ice areas comparable with those simulated for 2000–10 in the warmer years of the decade: a maximum bare-ice extent of 875 km² was predicted for the first 3 years of both simulations. However, in the cooler fifth and sixth years of the simulations, smaller maximum bare-ice areas were predicted by the warm climate scenario than by the 2000–10 simulation;the warm GrIS scenario predicted maximum bare-ice areas of 725 km² for both 2104 and 2105 compared with 750 km² simulated for 2004 and 775 km² in 2005. However, the predicted maximum bare-ice extent for 2106 exceeded that simulated for 2006 by 250km². In all subsequent years a greater bare-ice extent was modelled for the warm climate scenario than for 2000–10: 3.125% more of the total ice-sheet sector was exposed in 2107 than 2007, 18.75% more in 2108 than 2008, 21.88% more in 2109 than 2009 and 9.38% more in 2110 than 2010. Bare ice persisted for longer each year under warm GrIS conditions, and in 2112–13 a bare-ice zone exceeding 400 km² even persisted throughout winter.

Fig. 8. Bare-ice areas as predicted by the warm GrIS, cool GrIS and present-day simulations.

Figure 9 shows that overall the supraglacial carbon flux was more net autotrophic in warmer climates (mean NEP of –21 GgC a^-1 for the warm GrIS simulation compared with –16GgCa^-1 for the 2000–10 simulation) and a much greater volume of carbon was fixed by primary producers during the warm GrIS simulation (41.3 GgC a^-1) than for the 2000–10 period (35.94GgCa^-1). In the cool GrIS simulations, the net ecosystem production was consistently closer to balance than the other simulations (mean NEP –4.8Gg Ca^-1) and even shifted into a net heterotrophic state during the particularly cold period 2114–17. In the same years the contribution to the total carbon flux by cryoconite was very high, peaking at 59.91% in 2115.

Fig. 9. Net carbon flux as simulated by the present-day, warm GrIS and cool GrIS scenarios.

Discussion

Spatial heterogeneity in ecosystem characteristics

The diameter of cryoconite granules was observed to increase with distance from the ice margin (data not shown), indicative of longer residence times and less grain disintegration due to melt runoff. Large grains have small surface-area/volume ratios, so a much greater surface area will be produced by small grains in a single grain layer (Cook and others, 2010). Since the surface area of a given mass of sediment can vary, normalizing primary production rates for mass may not be appropriate. Additionally, the mass of surface algae in a particular area is usually negligible. Sediment surface area is related to primary production in a more predictable manner and was therefore used in preference to mass for normalizing primary production values for both cryoconite and surface algae in model 4. Coverage, mass per unit area and rates of biological activity for cryoconite sediment all peaked beyond 5 km (the extent of the field transect of Hodson and others (2010)), with most parameters exhibiting a positive correlation with distance from the ice margin. Similarly for surficial algal communities, coverage, mass and productivity were observed to vary with distance from the ice margin. This indicates that for realistic quantification of these variables, field data must be obtained from transects that cross the entire ablation zone. This is supported further by Stibal and others (2012) who suggest that organic carbon accumulation increases with distance from the ice margin. High spatial variability in all measured parameters suggests that fixed values are unlikely to yield realistic estimates of supraglacial carbon fluxes.

Incorporating new field data into carbon flux models

The carbon flux predicted by model 1 was approximately twice that predicted by Hodson and others (2010), probably because Hodson and others (2010) did not include the maximum cryoconite coverage and highest rates of respiration and photosynthesis found beyond the ice-sheet margin in the present study. Model 1 differs from the model of Hodson and others (2010) only in the spatial extent of the field data, yet interannual variability in carbon flux is apparent in the output of model 1, which was not resolved by the model of Hodson and others (2010). This, combined with the large discrepancy in total carbon flux, highlights a necessity for field data derived from transects crossing the entire ablation zone for supraglacial carbon flux modelling.

Model 2 estimates carbon fluxes using the same formulae as model 1;however, only a 1600 km² section of ice sheet was considered. Although model 2 includes field measurements from our 79 km transect, a greater carbon flux is predicted by the spatially variable model 3, which was also constrained to the 1600 km² ice-sheet sector (model 3 predicted a carbon flux ~30 times greater than model 2 for the same area). This indicates that using spatially variable parameter values instead of fixed values has a notable influence upon carbon flux predictions. This could be exacerbated if the model is upscaled to predict carbon fluxes for the entire ice sheet, highlighting the requirement for field measurements from a wider range and greater number of locations to inform future models. NEP predicted by model 3 shows significant interannual variability which is not resolved by the model of Hodson and others (2010) (Fig. 5). In 2001, model 3 predicted a particularly strongly net autotrophic NEP;however, in the same year the model of Hodson and others (2010) predicted net heterotrophy. An unexpected feature of the model was that between 2000 and 2006 a shift in NEP from net autotrophy towards balance was observed, which broadly coincided with decreasing bare-ice extent, suggesting that bare-ice area exerts a control over NEP on the GrIS, which is obscured in spatially constant models.

For both photosynthesis and respiration, the model of Hodson and others (2010) estimates greater carbon fluxes per square kilometre than model 3 (Fig. 6). For example, Hodson and others (2010) predicted the photosynthetic carbon flux in 2003 to be about three times greater than did model 3, while the carbon flux resulting from community respiration was seven times greater. This illustrates the dependence of carbon flux estimates upon spatially variable model parameters. The greatest discrepancy between the spatially uniform model of Hodson and others (2010) and model 3 coincided with the greatest bare-ice area, and the models tended to converge as bare-ice area decreased. This suggests that in the event of future climate warming, spatially constant models are likely to predict the carbon flux from the GrIS ever more inadequately.

The spatially variable model 3 represents the most sophisticated model of cryoconite-only carbon fluxes from the GrIS to date. The model output suggests significant interannual variability in carbon fluxes linked to bare-ice extent and therefore ultimately global temperatures. The introduction of spatially variable model parameters and field data spanning the entire ablation zone has been shown to increase the interannual variability and resulted in carbon flux estimates of much greater magnitude.

A further major outcome of this work is that model 4 predicts strong net autotrophy when surface algae are included, such that the magnitude of total carbon fluxes far exceeds previous estimates based upon cryoconite only. Our field observations showed that rates of photosynthesis and respiration were actually lower in the algal samples than in cryoconite samples (mean primary production rates were 360 mgCm^-2d^-1 for algae and 8.7mgCm^-2d^-1 for cryoco- nite);however, algae are able to colonize large areas of ice, in contrast to cryoconite which is concentrated in cryoconite holes. The surficial algal ecosystem has been reported in previous studies (e.g. Yoshimura and others, 1997;Uetake and others, 2010);however, the contribution to global carbon fluxes has hitherto been overlooked. Modelling results presented here confirm the importance of surficial algae for global carbon fluxes, suggesting for the first time that surficial algae contribute a far greater carbon flux than cryoconite. Additionally, the inclusion of carbon fluxes associated with surface algae makes the section of the GrIS strongly and robustly net autotrophic. The increased estimates of total carbon fluxes when surface algae are included indicate that supraglacial carbon fluxes should become significant components of global climate models.

Marginal sediment contributes only 0.001% (0.00017GgCa^-1) of the total NEP, suggesting that it could justifiably be omitted from carbon flux models in future. This is due to the narrow (~100m) marginal zone and low productivity. Even a hypothetical 1 km marginal zone contributed only 0.01% (0.0017 GgCa^-1) of the total carbon flux from the ice-sheet section.

The correlations between NEP for cryoconite and algae in Figure 7 suggest external as opposed to community-specific controls upon the productivity of the entire supraglacial zone. This is emphasized by the distinct difference in habitat of the two ecosystems, with cryoconite being sheltered from meteorology and high light intensities by melting into the ice, whereas the surficial algal ecosystem is not. A close coupling between climate and carbon flux is therefore evident. There is a negative correlation (R = –0.88; p < 0.05) between the total NEP and bare-ice area, indicating stronger net autotrophy with greater bare-ice areas. Strong correlations between bare-ice area and algal autotrophic carbon flux (R = 0.95;p<0.05), algal heterotrophic carbon flux (R = 0.98;p<0.05), cryoconite autotrophic carbon flux (R = 0.77;p <0.05) and cryoconite heterotrophic carbon flux (R = 0.82;p<0.05) suggest that the extent of the bare- ice area exerts a first-order control upon the carbon flux from the GrIS, also supported by the results of multivariate analysis (Stibal and others, 2012). This correlation is most likely stronger for algal communities because cryoconite holes are sheltered and adjust to increasing solar irradiance by melting into the ice (e.g. McIntyre, 1984). Hodson and others (2010) showed how this decouples the rate of photosynthesis from variations in the amount of light available at the ice surface at timescales longer than a few hours.

The correlation between NEP and bare-ice area suggests that climate warming may promote carbon fixation on the GrIS. It is therefore feasible that this may initiate a feedback whereby the fixation of greater volumes of atmospheric CO₂ may encourage biomass increase, lowering ice-sheet albedo and promoting enhanced melting upon the ice sheet. Additionally, there is likely to be a threshold melt rate at which biota is washed from the supraglacial zone, ‘cleaning’ the ice-sheet surface (Stibal and others, 2010) and altering nutrient concentrations (Uetake and others, 2010). This suggests that the response of ice surfaces to changes in NEP is complex; however, following further research it might be predictable. In order to further explore the link between climate and supraglacial carbon fluxes, climate predictions for the years 2100–10 were developed.

Climate scenarios

The modelled scenarios presented here represent a first attempt to model carbon fluxes from ice surfaces in future climates. Although simplistic, the climate scenarios do provide interesting outputs which support assertions made elsewhere in the literature that future climate warming could significantly influence carbon fluxes associated with supraglacial ecosystems. Our modelling suggests that a warmer climate generally results in both greater carbon fixation by primary producers and greater carbon metabolism by heterotrophs; however, the system tends to become more net autotrophic. This is largely due to the greater biomass of surface algae, which photosynthesize at a high rate but have relatively low heterotrophic productivity. This is supported by the greater contribution from surface algae in the warmer climate simulations compared with the cooler or present-day simulations. NEP, which is more strongly autotrophic, may therefore be associated with biomass increase, so a decrease in albedo upon the surface of the bare-ice zone seems likely

Conclusions

We have presented new carbon flux estimates for biological production upon the GrIS using more sophisticated models than have been published previously. New field measurements from a 79 km transect along Leverett Glacier (near Kangerlussuaq, Greenland) allowed minimization of upscaling errors and modelling of the carbon flux from a small 1600 km² section of the entire ablation zone that surrounded the transect. Spatial variability was included in input data for the first time, with the effect of increasing the predicted annual carbon flux by about four times that predicted using fixed values for the same area. Surficial algal communities (which have hitherto been ignored) have also been found to be the primary contributors to supraglacial carbon fixation. In our model, the supraglacial algal ecosystem fixes about 11 times more carbon per year than cryoconite under contemporary climate conditions. The GrIS is therefore thought to be in a relatively stable state of net autotrophy and only likely to reach balance or net heterotrophy in the event of severe global cooling. A strong link between bare- ice area and carbon fixation is evident, so climate predictions for the year 2100 show that should the GrIS warm over the next 100years, greater carbon fixation could be expected. This means that carbon fluxes mediated by supraglacial microorganisms should be included in regional models of the carbon cycle in polar regions and that potential feedbacks upon the albedo of ice-sheet surfaces resulting from increased biomass should be explored.