Storage and export of microbial biomass across the western Greenland Ice Sheet

The Greenland Ice Sheet harbours a wealth of microbial life, yet the total biomass stored or exported from its surface to downstream environments is unconstrained. Here, we quantify microbial abundance and cellular biomass flux within the near-surface weathering crust photic zone of the western sector of the ice sheet. Using groundwater techniques, we demonstrate that interstitial water flow is slow (~10−2 m d−1), while flow cytometry enumeration reveals this pathway delivers 5 × 108 cells m−2 d−1 to supraglacial streams, equivalent to a carbon flux up to 250 g km−2 d−1. We infer that cellular carbon accumulation in the weathering crust exceeds fluvial export, promoting biomass sequestration, enhanced carbon cycling, and biological albedo reduction. We estimate that up to 37 kg km−2 of cellular carbon is flushed from the weathering crust environment of the western Greenland Ice Sheet each summer, providing an appreciable flux to support heterotrophs and methanogenesis at the bed.


Introduction
This Supplementary Information documentation includes additional descriptions of the ice surface hydro-meteorological assessments employed in our study, the flow cytometric enumeration protocols and gating procedures employed in our flow cytometric analyses, the platform employed to acquire digital imagery of the ice surface in western Greenland, and a detailed explanation of our upscaling approach to estimate microbial cell and carbon fluxes over a wider region of the Greenland Ice Sheet.

Supplementary Method S1: Derivation of ice surface aerodynamic roughness lengths (z0) and hydraulic conductivity (K)
The aerodynamic roughness length (z0) of a glacier's surface imparts a control on the turbulent heat transfer from the atmosphere to the ice. Therefore, for estimates of glacier surface melt rates, many numerical models require z0 as a parameter. Physically plausible parameterizations of z0, based on the geometry of glacier surface topography (or roughness) and associated effective drag for turbulent flow, have been effective S1-3 . These traditional microtopographic methods typically rely on measurements of surface elevation along a transect of length L, perpendicular to dominant wind direction. The linear set of surface elevation measurements are detrended and broken into sections protruding above the mean transect elevation value. These sections represent obstacles to air flow, with a frequency (f) defined as the number of continuous groups of positive elevation deviations above the mean elevation over length L. The method S1-2 estimates the resistance to air flow based on the average obstacle height (h*) and the vertical silhouette area (s) and specific area measured over the horizontal plane (S) of the average obstacle over the transect. Because s and S are both defined by the ratio of L : f ref.S1 , the parameter z0 can be derived as: where, as an approximation, the mean obstacle height (h*) is estimated as twice the standard deviation of elevation (z) for the detrended transect. For transect lengths between 3 and 15 m, the estimate of z0 is independent of length S3 .

Supplementary Method S2: Derivation of near-surface weathering crust hydraulic conductivity (K)
Near-surface hydraulic conductivity is commonly assessed using piezometer tests that, following a disturbance to the water level, quantify the hydrological recovery of a shallow auger hole S4 (Amoozegar and Warwick, 1986). In the ice surface environment, positive disturbance through introduction of water to an auger hole would artificially raise the local water table. Recognising the bulk density gradient in the weathering crust, such a disturbance would likely induce water flow through higher porosity ice and an over-estimate of hydraulic conductivity. Therefore, a negative disturbance in which water is evacuated and hydraulic recovery monitored (the bail-recharge method) was considered more appropriate for this study.
We utilized a Kovacs ice auger to drill shallow holes of up to 36 cm depth in the weathering crust, and a biOrb TM syphon was employed to evacuate (bail) water from the hole. The bailed auger holes were immediately instrumented with bespoke electronic capacitance piezometers, recording the rising water level (recharge) at 2 s intervals with a Track-It USB data-logger (Monarch Instruments, USA). Our piezometers function by recording the voltage derived from the frequency of an electrical current applied across the auger hole, within the piezometer housing, which scales in proportion to capacitance: output voltage increases as the water level rises and capacitance is reduced.
Water recharge curves, plotted as water volume against time, were characterised by three stages.
The first, rapid quasi-linear recharge stage, relates to pressure driven water flow into the auger hole (Stage 1). The second, non-linear stage (Stage 2), reflects the flow of water through an undisturbed weathering crust. The final, linear stage with a gradient of zero is the point at which the auger hole water level equilibrates with the weathering crust's water table (Stage 3). From these curves, hydraulic conductivity (K) can be approximated empirically following the solution proposed by Bouwer and Rice S5 where length scales are given in centimetres: where Qr represents the recharge water flow rate (cm 3 s -1 ), defined by the auger hole geometry and rate of water column recovery measured by the piezometer and derived from recharge curve Stage 2; dr is the depth of the auger hole through which water enters, as defined by Stage 3, and rh is auger hole radius (2.5 cm). Dwc is the depth between the water table and the impermeable base of the weathering crust: for this study, we set Dwc = 40 cm, and such that [(Dwc -dr) / rh ] ≤ 6 ref. S5 . A and B defined here as 2 and 0.25, respectively, are dimensionless constants related to the ratio of hr:rh. The final variable, Dwt, is the depth between the observed water surface and equilibrium water table, approximated to be immediately below the static water table such that Dwt > 0. Numerical simulations accounting for error arising from measurement precision and user determination of recharge stages revealed a mean uncertainty in K of 4.8% S6 . Recharge curves in instances where either water level failed to recover following disturbance, or Stages 2 and/or 3 were indistinguishable were defined as 'incomplete' and reliable derivation of K was impracticable.

Supplementary Method S3: Flow cytometry protocol for the enumeration of microbial cells in supraglacial meltwaters
Flow cytometry (typically abbreviated to FCM, as here) has become a well-established and valued tool for the study of microbial communities in aquatic ecosystems S7 and has been employed successfully in glaciological studies S8-S10 . Such analyses can enumerate total particle loads, provide estimates on microbial abundance and size distributions, and combined with specific molecular or biological probes and/or autofluorescence is ideal for the detection of individual cell types against the backdrop of organic and inorganic particulates. Here, each sample was divided into two 1 mL aliquots, one of which was stained using SYBR Gold (1 × final concentration) whilst the other remained unstained, acting as a negative control. Both aliquots were stored in the dark at 20 °C for 20 to 240 minutes prior to enumeration. The FCM enumerations were undertaken using a Sony SH-800EC cell sorter, configured with 488 nm forward and back scatter detectors (FSC and BSC, respectively) and six fluorescence channels, including fluorescein isothiocyanate (FITC: 525 ± 30 nm).
The instrument was aligned and calibrated daily using 8-peak reference beads (Sony Biotechnology, Japan). An instrumental trigger threshold was established in the FSC-Area (FSC-A) channel using 1 μm size-calibration beads (Molecular Probes, USA) suspended in Ultrapure Water (Elga), ensuring that micron-scale particles were detected. The threshold and gain of the trigger channel was set at the lowest possible point at which no noise was observed whilst ensuring the detection of the 1 μm beads. A log-log BSC-A against fluorescence channel 2 area (FL2: FITC-A) plot was used to define the bacterial gate, via the comparison of unstained (Supplementary Fig. 1a) and stained ( Supplementary   Fig. 1b) sample pairs. Size gates were defined from a non-fluorescent size calibration kit (Molecular Probes, USA) enabling discrimination of microbial cells into seven size categories from the FSC channel, and bacterial size distribution was measured using these gates, as per the manufacturer's instructions ( Supplementary Fig. 1c). Gates were kept consistent throughout FCM sample runs. Fig. 1: Illustration of the gating protocol used to enumerate microbe abundance and size fractions with example unstained sample (a), the same sample with SYBR Gold stain (b) and size fractions defined by the FSC-A (c).

Supplementary Method S4: Image acquisition using an Unmanned Aerial Vehicle (UAV)
The fixed-wing UAV used in this study was identical to those described in detailed by  .
The UAV platform comprised a 2.1 m wingspan airframe, a 10 Ah 16.8 V LiPo battery, an Arduino autopilot module, and included a downward-facing Sony NEX-5N digital camera, totalling a platform weight of ~ 4 kg. The autonomous control system of the Arduino navigation and flight computer is updated in real-time. A horizontal displacement threshold was used to automatically trigger the camera shutter. The camera has a 16 mm fixed-focus lens, which yields a field-of-view of 53.1° by 73.7°, and was pre-set for a fixed exposure of 1/1000 s, ISO-100 and F-stop of 8. At an elevation of 350 m above the surface, the imaging footprint was approximately 525 x 350 m, with a pixel equating to a horizontal ground sampling resolution of 0.11 m. The digital RGB images, recorded in RAW format, were georeferenced using the latitude, longitude, altitude and attitude data recorded by the onboard flight controller. Pre-planned flight lines over the ice sheet, were uploaded to the UAV using the ArduPilot Mission Planner ground station software (http://ardupilot.org/planner/) and the image acquisition sortie was conducted on 8 August, 2014.

Supplementary Method S5: Derivation of a seasonal export of cellular carbon from the weathering crust interfluve
The porous weathering crust develops through subsurface melt arising from shortwave radiation penetration and heat transfer by percolating meltwater. However, to date, this process remains poorly constrained outside Antarctica S14 . Therefore, as a first order estimate following Muller and Keeler S15 , we simply assumed that the weathering crust in western Greenland deepens from the initial exposure of bare-ice in early summer when short-wave radiation is elevated; in the later summer, as the shortwave radiation is reduced relative to sensible heat fluxes, the weathering crust thins and declines.
To determine the duration bare ice exposure, we used the data presented by Ryan et al. S16 that used the MODIS MOD09GA and MOD10A1 products to define a bare ice exposure frequency over the 90day summer season (JJA) centred on 15 July (DOY196). From these data, we extracted the mean bare ice duration (bi) for 795 moulin terminating catchments in the low gradient, moderately crevassefree ablation area region ranging from 400 m to 2000 m a.s.l. as delineated by Yang and Smith S17 .
The western Greenland region reported here covers a total area of 13847 km 2 , with a mean internally drained catchment area of 17.4 km 2 . Our ~0.5 km 2 S6 study catchment lay within a larger 43.8 km 2 moulin-terminating catchment that, in 2014, exhibited an areal mean bare ice duration of 68 days.
In the absence of any data pertaining to weathering crust development, we restrict its thickness to 0.4 m with a maximum hydraulically active depth (Wmax) of 0.29 m, as defined by our recharge experiments and observations of the near-surface water table. To model the seasonal evolution of the hydraulically active layer we used a cosine function. To define this function, we employed the 68-day bare ice duration for the wider S6 study catchment, and our observation period from DOY 204-211 throughout which Wmax was assumed constant (Supplementary Fig. 2). This approach yielded a 30-day period with the hydraulically active layer assumed to be Wmax, and 19-day phases of weathering crust growth and decay. For bare ice durations less than 38 days, we assumed a linear increase in the depth of the weathering crust's hydraulically active layer based on bi.
To define the water discharged from the supraglacial interfluve area over the bare-ice period (bi), we integrate the depth function over t days, and multiply by the weathering crust effective porosity ( = 0.22), our derivation of throughflow velocity (vt) and the surface fraction (Sf) we expect to represent ice adjacent to a stream: Based on our study catchment, the interfluve stream bank area across the region was kept constant at Sf = 14%, and the bare ice duration that sees the hydraulically active weathering crust layer peak but not remain at Wmax (bp) was 38 days.
Supplementary Fig. 2: Illustration of the evolving hydraulically active weathering crust depth for varied bare ice durations at 15-day intervals used to simulate the seasonal growth and decay of the porous surface ice layer. The 90-day summer (JJA) season is centred on 15 July (DOY 196). The 2014 field campaign, from DOY204-211, is shown in grey with the hydraulically active weathering crust depth as derived for the 68-day bare-ice duration observed in the region of our study catchment. The maximum mean bare-ice duration observed in the catchments in the western sector of the ice sheet was 80 days.  Supplementary Table 1). Our approach here provides recognition of the latitudinal and elevational variation in weathering crust development, and associated drainage evolution. The derivation of an active stream network from remote sensing data will be conditioned by the resolution of the imagery employed, as well as the seasonal evolution of the surface drainage density; these variables are beyond the scope of this study. However, it is clear more detailed modelling of both weathering crust development S14 in the context of Greenland, and determination of seasonal evolution in drainage density as well as the temporal and spatial variation in both the intrinsic and extrinsic controlling factors influencing the hydraulic conductivity of the weathering crust S6 are required to further refine our estimates of carbon fluxes at wider regional scales.

Supplementary Method S6: Numerical analyses and figure design
Commercially available software utilised in the numerical analyses, and development of illustrations presented in the manuscript and Supplementary Information Table 1: Seasonal estimates of microbe abundance and associated mass of cellular carbon (in kg C) exported from the weathering crust to the supraglacial stream at the S6 study site in 2014, and to the subglacial environment in the western sector of the ice sheet based on moulin terminating catchments extending over an area of 1.38 ×10 -4 km 2 . Results are also shown for low (2006) and high (2012) melt years. Biomass estimates are made for the ≤ 15 µm size classes and for all classes including large (>15 µm) algae, using low S18,S19 and high S20,S21 cell constants, allometric S22 and constant S23,S24 ratios.