Arctic mid-winter phytoplankton growth revealed by autonomous profilers

Sensor payload. ProIce floats are equipped with a CTD (conductivity-temperature-depth) Seabird41 unit, an Aanderaa 4330 oxygen optode, a WETlabs remA [combining an ocean color radiometer (OCR504) (380, 412, and 490 nm) and PAR] and an Environmental Characterization Optics (ECO) Triplet for the observation of chlorophyll a fluorescence, colored dissolved organic matter, and particle backscattering at 700 nm (bbp), as well as a Satlantic SUNA V2 sensor for nitrate measurements.

Design and mission management. The Argo float model used in this study is the so-called ProIce float, manufactured by the French company NKE Instrumentation. It is based on the Provor CTS5 developed under a collaboration between NKE and LOV (34). This float has the same mechanical characteristics as the previous version (CTS4) and can carry the entire suite of biogeochemical sensors. It does not require preballasting as a function of mean seawater density, a substantial advantage in Polar Regions where low-density surface layers are present (35). Mission parameters are highly flexible and can be changed both in real time by Iridium communication and in advance by means of a script file. The latter has been used in this study to modify the float sampling schedule under ice without any satellite communication.

Adaptation of the CTS5 float to Arctic conditions was carried out in a close collaboration between LOV and Takuvik Laboratories in the framework of the Novel Argo Ocean Observing System (NAOS) project. Takuvik lead field trials in a frozen lake close to Québec city and at the Green Edge ice camp in Qikiqtarjuaq (Nunavut) to validate the reliability of the float and sensors in polar conditions. Lagrangian trajectories were simulated to choose the best dropping zones, combined with observations from historical climatology and ice charts.

Ice avoidance. Argo floats must surface for geolocalization, data transmission, and command reception using satellite networks (35). Sea ice, which can damage floats, must be avoided to make the floats operational in the Arctic Ocean (36). If sea ice is present at the surface, Argo floats need to postpone surfacing. It then has to perform several consecutive profiles storing data during wintertime.

ProIce floats use a combination of three technologies to prevent surfacing in the presence of sea ice: an upward looking Teledyne PSA-196 altimeter, an ice sensing algorithm (ISA) (37) based on local temperature, and a simple date criterion. The altimeter was used to compare the depth to the distance to the surface to detect ice cover. Because of the uncertainties on both the depth and the distance, a threshold of 5 m was used, which corresponds mainly to the presence of an iceberg. This triggering never occurred for floats deployed in the present study. We used a substantial database of temperature and salinity profiles and ice cover observations in Baffin Bay to locally adapt ISA’s parameters. The floats were hence programmed to abort surfacing whenever the median temperature between 30 and 10 dbar was below −1.3°C. As Arctic surface waters exhibit high variation of salinity, we assume that these thresholds are highly local and should be adapted before being applied to other Arctic seas.

Last, we prevented surfacing between November and July when the probability of sea ice was very high. This was assessed on the basis of a survey of climatological sea ice cover data. An optical device, specifically designed by Takuvik to detect sea ice based on laser polarimetry (38), was tested with promising results but not used in operation.

Since floats did not sample shallower than 10 m when ice was present, some phytoplankton biomass was likely missed in July, which explains the fact that the spring/ice-edge bloom is not more pronounced in Fig. 2E.

Sampling schedule. Briefly, the floats acquired daily profiles in the ice-free summer period, every fourth week during winter and every 10th day in spring (table S2). Profiles were terminated at 15:00 UTC each, i.e., approximately local noon.

In mid-November, the floats entered a winter cycle of one profile every 14 (float 020b) or 28 days (other floats) ending in late June the following year. During this period, the up-casts were limited to 10 dbar to prevent any collisions with sea ice on the surface. Starting in July, profiles were collected every 10th day; the floats surfaced after each profile if not prevented by the ice detection system. Daily profiling was programmed from 1 August until 15 November.

Float sensor validation. RemA sensors (OCR504 and ECO Triplet) data were collected under dark conditions before deployment to evaluate their offset drift since factory calibration. A night profile is programmed as an in situ check of this offset down to 1000 m to meet relatively large temperature dynamics, helping to see the response of OCR to temperature.

Analysis onboard for cross-validation: A 0- to 1000-m CTD rosette cast was systematically performed after the deployments, and water was sampled to provide an evaluation and characterization of the calibration errors of the biogeochemical sensors (39): 10 levels in the 100-m upper layer for shipboard chl-a measurements, as well as for high-performance liquid chromatography (HPLC) for chlorophyll concentration performed by the SAPIGH (Service d’Analyses de PIGments par HPLC) unit in France for the validation of the chl-a fluorometer, and 10 levels in the water column for FDOM analysis onboard by means of an Ultrapath unit for validation of the CDOM fluorometer.

Chlorophyll a concentration was determined with a WETLabs ECO chlorophyll a fluorometer. Raw chl-a fluorescence was smoothed with a five-point median filter to remove noise and spikes and then dark-corrected based on deep values (40). In the third step, nonphotochemical quenching was corrected (40), and last, the slope (F-factor) was determined on the basis of the comparison between the first float profile and on-board fluorometry. Corrected chl-a fluorescence was calculated as F · (FChlaraw − FChlaDark). Each profile of bbp data was smoothed using a five-point median filter.

Mixed layer depth. Mixed layer depth was calculated as the shallowest depth where density (σ_θ) exceeded the surface density by at least 0.1 kg m⁻³ (41). Surface density is defined as the 0 to 15 m average value to be consistent between summer and winter, when the shallowest sampling depth was 12 m. The density criterion of 0.1 kg m⁻³ is higher than those commonly used at lower latitudes (42) due to strong Arctic stratification.

Light model. A radiative transfer model (43) was used to complement float measurements to compute diurnally integrated PAR profiles along the float trajectory following (44). Briefly, the model provided theoretical clear sky PAR values just below the sea surface with a time increment of 1/60 of the day length. From these values, instantaneous PAR (iPAR) profiles were created using a diffuse attenuation coefficient derived from the chl-a profile measured by the float (45). By comparing the float iPAR profile with the one modeled for the same time and location, the light attenuation due to the cloud cover and the sea ice cover was estimated. The mean, over the vertical profile, of the ratio of float iPAR to modeled iPAR, weighted by absolute values of float iPAR, gave a correction factor (fig. S1) that was subsequently applied to all the modeled profiles of the day. The underlying assumption was that both the average cloud and sea ice cover estimated during the float ascent were representative of the day. The vertical profile of daily integrated PAR was then computed by time integrating all these corrected iPAR profiles over the day length.

Specific phytoplankton growth rates. Net specific phytoplankton growth rates (d⁻¹) were calculated on the basis of surface layer average and integrated values as follows, adapting the methodology developed in (13). For each profile, the surface layer sfc was defined, for each profile, as the larger of z0.4 [the 0.4 mol m⁻² d⁻¹ isolume (12)] and the mixed layer depth, which, in practice, meant the mixed layer depth throughout winter and the isolume depth through summer (fig. S2). chl-a and bbp were then both integrated (∫_sfcdz⋯) and averaged (⟨⋯⟩_sfc.) over the surface layer.

An isolume depth, i.e., the depth of a specific daily irradiance level, is commonly used (12, 44) to bound the vertical extent of net phytoplankton growth. We chose z_0.4 as the value previously found both to vertically constrain the extent of phytoplankton net growth in subtropical oceans (12) and to terminate the Arctic fall bloom in Baffin Bay. We hence assume that this value takes into account grazing rates typical of summer. Because during winter, the mixed layer depth was always deeper than any reasonable choice of an isolume depth, we do not need to determine a corresponding isolume criterion for winter.

Biomass changes were primarily reflected by changes in mean concentration, with the exception of possible convective dilution during winter, when the mixed layer deepened. To calculate the “net specific growth rate,” we stipulate, still following (13), that during mixed layer deepening in winter, the total phytoplankton biomass is represented by the integrated chl-a or bbp to account for dilution by entrainment of low biomass water from below the mixed layer. As chl-a and bbp observed in Arctic winter are extremely low (figs. S3 and S4), biomass below the mixed layer is frequently not significantly less than in the mixed layer, and hence, increases in the integrated biomass only reflect the deeper integration limit not a real signal. Accordingly, changes in the mean biomass should be used to infer growth rates then as well. When the mixed layer is shoaling or the zone of net phytoplankton growth (here stipulated as being bounded by the 0.4 mol m⁻² d⁻¹ isolume) extends deeper than the mixed layer, detrainment (loss) of higher biomass water from the mixed layer is happening, and so, biomass changes are reflected by changes in mean concentration too.

In summary, the net specific growth rate r was calculated from the integrals

if and only if all three of the following conditions are fulfilled: (i) if the mixed layer was deepening (ΔMLD > 0, where MLD stands for mixed layer depth), hence if entrainment can have taken place at all; (ii) the mixed layer was deeper than the 0.4 mmol photons m⁻² d⁻¹ isolume (MLD > z0.4); and (iii) there was a statistically significant (P < 0.05) difference between mixed layer average biomass (bbp or chl-a) and the biomass averaged from the mixed layer base to 15 m below the mixed layer (see computer code).

In all other cases, r was calculated using the difference of the means

To recapitulate, criteria (i) and (ii) are following a published methodology (13), while we have added the restrictive criterion (iii) to account for weak winter biomasses that may render invalid the original reasoning. To account for values that may be at or below the noise levels, we further discarded any growth rate where the surface layer mean biomasses in both involved profiles was below the environmental background (see below for the calculation of environmental background variability level).

The entire growth rate calculation procedure is illustrated in figs. S5 to S12. Biomasses had to be averaged in bins before calculating the growth rates (see below for a more detailed explanation). Taking as an example a bin width of 14 days, starting on 20 July 2017, we have a total of 180 estimates of the specific change of phytoplankton biomass (14 time periods × 4 floats × 2 biomass proxies, i.e., chl-a fluorescence and bbp, but not every float was reporting during every bin). Of these, 12 were rejected because the associated biomasses were not significantly different from the environmental background variability. For the remaining 168 values, the net growth rate r was calculated on the basis of the surface layer integrals 33 times or 19% of the time, while the rest (81%) were calculated on the basis of the differences in mean values.

This methodology does not take into account that active mixing may not extend to the base of the mixed layer at all times (49), but most temperature and salinity profiles during winter exhibited steep density gradients at the mixed layer base (figs. S13 to S15), indicating strong mixing throughout winter.

As common parameterizations for biomass are linear in bbp and chl-a, respectively (50), the specific biomass growth rates do not depend on the numeric values of the conversion factors between chl-a, bbp, and biomass (in g C m⁻³). Further exploring the robustness of growth and decay patterns, we also wrapped the entire time series into a single annual cycle, smoothing it using a generalized additive model with 12 splines and calculating the growth rates as described above. The results can be seen in fig. S16 and Fig. 2H.

Sensitivity of specific net growth rates to biomass binning. The above algorithm calculates growth rates from successive differences of either surface layer integrated or surface layer averaged chl-a or bbp. Especially in summer, time series are very noisy, necessitating aggregating those surface layer mean or integrated biomasses further in, e.g., 7-daily or 28-daily bins. To address the sensitivity in the growth rate calculation to the exact choice of how we binned each time series, we varied the bin edges to conduct a parameter sweep. Specifically, bin width was, in turn, set to 7, 14, 21, or 28 days, and for each bin width, the bin edges were successively positioned in 1-day intervals. [That means, for example, for a bin width of 7 days, bin edges were placed on (1 July 2017, 8 July 2017, 15 July 2017, …) for one iteration, (2 July 2017, 9 July 2017, 16 July 2017, …) for another, and so on.] See fig. S17 for an illustration.

As the focus is on the winter period, when each float sampled at 14- or 28-day intervals, we also investigated the difference in growth rates based on biomasses calculated from biomass averaged in 14- and 28-day bins. (In practice, this means calculating each winter growth estimate from either successive profiles or from successive means of each two profiles.) For these a Kruskal-Wallis test confirmed that calculated growth rates did not significantly vary from one method to another (fig. S18).

Sensor noise levels, significance of low chlorophyll a and backscattering measurements, and other noise sources potentially contaminating winter biomass growth rates. Given the low biomasses measured in winter, one has to worry about the statistical significance of any results derived from them. In this section, we show that most noise sources were not strong enough to potentially contaminate the signal and shadow the real growth rate. Three noise sources have to be considered: noise stemming from the sensor itself, noise stemming from (random or small-scale) environmental variability, and environmental spatial trends (such as a potential steady drift into higher chlorophyll water masses during winter). We address each of these in turn.

First, let us consider sensor noise. WETLabs/Seabird ECO sensors (that measure bbp and chl-a) state a resolution of one digital count for both types of sensors of the ones mounted on our floats. This is defined as the SD of a sensor sampling a static target at 60 Hz for 60 s (i.e., a sample size of n = 60). For chl-a, the factory calibration is that 1 count = 0.0073 mg chl-a m⁻³, which after quality control using on-board HPLC drops by a factor of approximately 2 down to 0.0037 mg chl-a m⁻³. [This is the slope factor mentioned in table S3. This factor is also in line with a global validation of the ECO sensor (51) that recommended applying a factor of two to reduce all measurements using this sensor.] A float samples the ECO sensors at every 10 cm of depth and therefore the typical “surface layer averaged chl-a” value is composed of n = 300 samples for a surface layer of, e.g., a depth of 30 m. The standard error

, that is, the expected deviation of the sample mean from the “true” population mean, is therefore much lower than one count. The lowest mixed layer chl-a values we measured, on the other hand, were an order of magnitude larger and therefore very much significant. For bbp, 1 count corresponds to around 1.23 · 10⁻⁵ m⁻¹. The smoothed time series of averaged bbp displayed in new Fig. 2H climbs from 1.9 · 10⁻⁴ to around 2.5 · 10⁻⁴ from early March to mid-April, again a much larger increase than sensor noise.

To estimate the chlorophyll equivalent of such an increase in 700-nm backscattering, one can use literature values of the chlorophyll-to-backscattering ratio. Concretely, at 700 nm, backscattering can be decomposed into

, where

is the backscattering ratio at 700 nm and

is the chl-a specific scattering coefficient at 700 nm. The backscattering ratio at 700 nm is in the range of 10⁻⁴ to 3 · 10⁻³ (52), and the chl-a specific scattering coefficient at 700 nm is 0.1 to 0.3 m² (mg chl-a)⁻¹ (53). This means that one chl-a count (0.0037 mg chl-a m⁻³) corresponds to a backscattering increase as little as 3 · 10⁻⁸ to 3 · 10⁻⁶ m⁻¹, demonstrating that the chl-a signal is much more sensitive to the phytoplankton concentration than the backscattering sensor is. The environmental (nonphytoplankton) background of the backscattering signal can hence be expected to play a more important role when biomasses are low, as the next paragraph details.

To estimate the environmental background levels and their variability, we calculated the mean chl-a and bbp over a depth of 750 to 800 m of every profile (800 m is the maximum depth of Davis Strait and therefore the deepest depth that we can expect to be reasonably well ventilated and not potentially accumulate an excess of detritus), we find for chl-a an SD of 0.0025 to 0.0035 mg m⁻³ depending on the float, just below one digital count, again much smaller than the lowest average mixed layer values. For bbp, the SD is 1.5 to 8.7 · 10⁻⁵ m⁻¹, which is a more appreciable fraction of the winter increase in bbp but still relatively small. It also reflects the fact that bbp has a higher nonbiological background value than chl-a, as is also evidenced by the fact that mixed layer bbp values essentially floor at the bbp baseline in the new Fig. 2E (which is not the case for chl-a, which always stays much above its baseline). These 750 to 800 m averages are 2 · 10⁻⁴ ± 5 · 10⁻⁵ m⁻¹ for bbp and 0.0059 ± 0.0029 mg chl-a m⁻³ for chl-a (Fig. 2E). Overall drift of these background values (as determined by an ordinary least squares regression) were mostly statistically insignificant, and the associated drift in values over a year was always less than the SD.

Third, larger-scale changes in the environment have to be considered, such as the floats drifting, during the course of winter, into higher biomass waters. As Argo floats are never 100% Lagrangian, this may be an explanation for some short-term changes, besides spatial patchiness, but extremely unlikely to affect growth rates over several months: (i) Inferred growth rates are consistent between two consecutive years, reducing the chances of a statistical fluke; (ii) because water mass properties as sampled by the entirety of our Argo fleet did not systematically change over the course of the year (Fig. 1D); and (iii) because all floats exhibited increasing biomasses during most of winter even individually.

In summary, it is true that during the dead of winter, environmental variability in the backscattering data renders many of the backscattering-based growth rates unreliable, and these have been excluded from the analysis as described above. However, some remain above this noise level and give positive growth levels, albeit later than chl-a, which is always above noise levels and demonstrates growth as early as February (seen over both years). In other words, for biomasses as low as we have observed, bbp is not a sensitive enough proxy, and the zero growth rates are likely merely an artifact rather than part of the phenology.

Photosynthetic parameters. In addition to using literature values for summer Baffin Bay photosynthetic parameters (15), photophysiological experiments were conducted off the fast ice in Baffin Bay in 2015 and 2016 during the Green Edge ice camps (54). Methods are described in (55), but briefly, we determined E_k, the light saturation irradiance for photosynthesis, and P_{b, max}, the chl-a specific light saturated growth rate, using the models P = P_s · (1 − e^−αE) · e^−βE + P₀ (when photoinhibition was apparent) and P = P_s · (1 − e^−αE) + P₀ (when no photoinhibition was apparent). Here, P is the rate of photosynthesis and E is the light intensity; the rest are free parameters. P_{b, max} was then calculated by normalizing

to chl-a, and

. Values were averaged for winter (November through May), spring (June), and summer (July through October), see fig. S19.

Cell division rates. To investigate the relationship between phytoplankton net growth r and cell division rates μ, we modeled light-field modulated variations in phytoplankton cell division rates based on the previously calculated underwater irradiances and photosynthetic parameters. We calculate photosynthesis [in g C (g chl-a)⁻¹ d⁻¹] as (56)

and then μ as μ_b/θ, where θ is the carbon-to-chlorophyll ratio, here assumed to be 30 (g C) (g chl-a)⁻¹ (1, 57).

The original formula (56) gives an integral from sunrise to sunset, assuming that light is not sufficient for photosynthesis when the Sun is below the horizon. We use an integral over an entire cycle of 24 hours, weighting daily PAR by the sine of the solar elevation angle but an exponential decrease during twilight and night as described above. During winter, when PAR did not exceed noise values, we used the noise PAR profiles (derived above) to determine the upper limit of cell division rates that could have been possible given the combinations of iPAR noise floor at noon, shallowest measurement depth, and Sun angle.

As a reality check, we compared modeled cell division rates μ with the specific net phytoplankton growth rates r, showing good agreement for March through May when averaged over all floats (fig. S20), but considerable scatter when considering each float on its own (fig. S21).

Mesozooplankton. Mesozooplankton stratified net sampling in the Amundsen Gulf during the International Polar Year/Circumpolar Flaw Lead (IPY-CFL) 2007–2008 overwintering expedition (fig. S22) was described in (28). Biomass concentration was calculated from mesozooplankton counts and species-specific length-carbon content relationships. From three stratified nets (0 to 20, 20 to 40, and 40 to 60 m), biomass was integrated (g C m⁻²) and then divided by the 60-m thickness of the layer, expressing biomass as g C m⁻³). Figure S23 shows that Baffin Bay and Amundsen Gulf are similar in both abundance and mesozooplankton community structure.

Grazing experiments. Zooplankton samples for estimation of fecal pellet production were collected from the fast ice once every 2 to 7 days between 24 May to 7 July 2016 as a part of the Green Edge 2016 field campaign in the Baffin Bay using a 200-μm mesh net hauled vertically from 100 m (before 25 June) and from 30 m (after 25 June) depth to the surface and incubated following (58). Fecal pellets samples (100%) were counted, and its shape was measured under a stereomicroscope (×20 to ×40 magnification). Fecal pellet production rate (mg C m⁻² d⁻¹) was estimated on the basis of fecal pellet volume, gut clearance time of 33 min for Calanus glacialis stage CV and CIV (adult female) in Northern Baffin Bay (59), and a conversion factor from fecal pellet volume to particulate organic carbon contents as 0.048 mg C mm⁻³ (60). To remove the effects of variation in the volume filtered due to the difference of net towing depth (30 or 100 m), fecal pellet production rate was standardized for those data as towing distance of 100 m (i.e., the production rate per square meter = raw data × 100 per towing depth).