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 Reference TellingTelling 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 Reference HodsonHodson 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 Reference Nicholls and DillonNicholls 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 Reference NauwerckNauwerck (1963). Average mass per unit area was determined for each coverage category, and then weighted means calculated for each site.

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.

Reference HodsonHodson and others (2010) used an adapted monthly version of Reference Janssens and HuybrechtsJanssens 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 (Reference HannaHanna and others, 2011) in order to grow A_bio and i for the period 2000–08. The runoff model is described in Reference Hanna, Huybrechts, Janssens, Cappelen, Steffen and StephensHanna and others (2005) and has been used for several ice-sheet studies (e.g. Reference HannaHanna and others, 2008;Reference Sundal, Shepherd, Nienow, Hanna, Palmer and HuybrechtsSundal and others, 2009, 2011;Reference MurrayMurray and others, 2010). In the carbon flux model of Reference HodsonHodson 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, Reference HodsonHodson 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 Reference Anesio, Hodson, Fritz and Psenner Rand SattlerAnesio and others (2009). The field measurements of Reference HodsonHodson 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 Reference Anesio, Hodson, Fritz and Psenner Rand SattlerAnesio and others (2009) are representative, and rates of PP used by Reference HodsonHodson and others (2010) have been identified as anomalously low. High uncertainty therefore permeates these carbon flux estimates.

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 Reference AnesioAnesio 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 Reference Williams and LefevreWilliams 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 (Reference Cook, Hodson, Telling, Anesio, Irvine-Fynn and BellasCook 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, Reference Huybrechts, Gregory, Janssens and WildHuybrechts 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’.

Reference Hanna and CappelenHanna 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 (Reference HannaHanna and others, 2008). Reference Bromwich, Robasky, Keen and BolzanBromwich 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 Reference SmithSmith (1999). Extrapolating the data from Reference Hanna and CappelenHanna 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.

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 Reference HodsonHodson and others (2010). Figure 4a provides a comparison between the rates of primary production predicted by model 1 and Reference HodsonHodson and others (2010). The overall positive trend evident in model 1’s output is not apparent in the data from the model of Reference HodsonHodson 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 Reference HodsonHodson and others (2010).

Fig. 4. (a) CP^ice, (b) CR^snow and (c) CR^ice as predicted by model 1 and Reference HodsonHodson 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 Reference HodsonHodson and others (2010). In 2001, model 3 predicted a particularly strongly net autotrophic NEP;however, in the same year the model of Reference HodsonHodson 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 Reference HodsonHodson 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 Reference HodsonHodson and others (2010) in kgCkm^-2. In both, the model employed by Reference HodsonHodson and others (2010) estimates greater carbon fluxes per square kilometre. For example the carbon flux resulting from autotrophy in 2003 predicted by Reference HodsonHodson 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 Reference HodsonHodson 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 Reference HodsonHodson and others (2010); (b) carbon flux resulting from cryoconite community respiration as predicted by model 3 and Reference HodsonHodson 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 Reference HodsonHodson 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.

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 (Reference Cook, Hodson, Telling, Anesio, Irvine-Fynn and BellasCook 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 Reference HodsonHodson 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 Reference Stibal, Telling, Cook, Mak, Hodson and AnesioStibal 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 Reference HodsonHodson and others (2010), probably because Reference HodsonHodson 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 Reference HodsonHodson 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 Reference HodsonHodson 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 Reference HodsonHodson and others (2010) (Fig. 5). In 2001, model 3 predicted a particularly strongly net autotrophic NEP;however, in the same year the model of Reference HodsonHodson 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 Reference HodsonHodson and others (2010) estimates greater carbon fluxes per square kilometre than model 3 (Fig. 6). For example, Reference HodsonHodson 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 Reference HodsonHodson 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. Reference Yoshimura, Kohshima and OhtaniYoshimura and others, 1997;Reference Uetake, Naganuma, Hebsgaard, Kanda and KohshimaUetake 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 (Reference Stibal, Telling, Cook, Mak, Hodson and AnesioStibal 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. Reference McIntyreMcIntyre, 1984). Reference HodsonHodson 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 (Reference Stibal, Lis, Lawson, Mak, Wadham and AnesioStibal and others, 2010) and altering nutrient concentrations (Reference Uetake, Naganuma, Hebsgaard, Kanda and KohshimaUetake 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