in

Sea-ice derived meltwater stratification slows the biological carbon pump: results from continuous observations

Table S1 lists the data used in this paper, the instruments that it is based on, the data repositories, and in which figures the data are used.

Global data sets

Bathymetry

Bathymetric data was taken from the International Bathymetric Chart of the Arctic Ocean (IBCAO 30 sec V3)64 available at https://www.ngdc.noaa.gov/mgg/bathymetry/arctic/grids/version3_0/.

Sea ice concentration

We use data derived from the Advanced Microwave Scanning Radiometer sensor AMSR-2 for the years 2013–18 processed in accordance with65 and downloaded from https://seaice.uni-bremen.de/sea-ice-concentration-amsr-eamsr2/66. At each grid point the sum of days during all April/May/June of 2013–2018 when the sea ice concentration at the grid point was >20% was divided by the total number of days with data in those months to obtain the percentage of days with ice concentration >20% (Fig. 1). For separate 7-day periods in April/May/June 2017 and 2018 the mean ice concentration over those 7 days was calculated and the 20% contour of this mean was plotted separately for each of those 7-day periods. For each mooring and each day, the ice concentration at the grid cell closest to the mooring was calculated (Fig. 4a and S1a), and if the ice concentration at the mooring was below 20%, the shortest distance to grid cells where the ice concentration exceeded 20% was calculated (Fig. 4a and S1a). If the ice concentration at the mooring exceeded 20%, the shortest distance to grid cells where the ice concentration was below 20% was calculated and the distance was defined as negative.

Sea ice velocity and sea ice area export

Ice area flux estimates in Fig. 2a are calculated using CERSAT (Center for Satellite Exploitation and Research, France) motion estimates together with CERSAT ice concentration information67. Fluxes are estimated along a zonal gate positioned at 82°N between 12°W and 20°E and a meridional gate at 20°E between 80.5°N and 82°N (Fig. 1) for the period 1994–2020 (January–May). The ice area flux at the gate is the integral of the product between the meridional and zonal ice drift and ice concentration. For a more detailed description we refer to ref. 68. Arctic-wide sea ice velocity anomalies (Fig. 2b, c) were computed from the OSI-405-c motion product provided by the Ocean and Sea Ice Satellite Application 635 Facility (OSISAF)69.

Satellite chlorophyll

Surface chlorophyll concentrations measured with the Sentinel 3 A OLCI (Ocean and Land Colour Instrument) were downloaded from https://earth.esa.int/web/sentinel/sentinel-data-access. The 8-day satellite data were averaged for the time series over grid points within boxes of 60 km by 60 km around the moorings.

Atmospheric reanalysis

ERA-Interim reanalysis70 data at the surface on a 0.25° latitude by 0.25° longitude grid at 12 hourly resolution was downloaded from https://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/. Incoming shortwave radiation (ssr) and outgoing longwave radiation (str), sensible heat flux (sshf), and latent heat flux (slhf) were extracted and averaged to daily values.

Physical numerical models

FESOM

In this study, we used model data from the Finite-Element Sea ice-Ocean Model (FESOM) version 1.471. FESOM is a sea ice-ocean model that solves the hydrostatic primitive equations for the ocean and comprises a finite element sea ice component. It uses triangular surface meshes for spatial discretization, allowing for a refined mesh in regions of interest, while keeping a coarser mesh elsewhere. In the model configuration used here, a mesh resolution of nominally 1° was applied in the global oceans. The mesh was refined to 25 km north of 40°N, and to 4.5 km in the Nordic Seas and Arctic Ocean. In the wider Fram Strait (20°W-20°E/76°N-82°30′N), the mesh was further refined to 1 km. In this region, the simulation can be considered as eddy-resolving, as the local internal Rossby radius of deformation is about 2–6 km72,73. In the vertical, the model used 47 z-levels with a resolution of 10 m in the upper 100 m, and coarser resolution with depth (with a resolution of ~100 m at 800 m depth). For bottom topography, the RTopo-2 data set was used74. The model simulation covers the period 2010–2018 and has daily model output. It was forced with atmospheric reanalysis data from Era-Interim70, and was initialized with model fields from the simulation described in ref. 75. River runoff (except for Greenland) was taken from the JRA-55 data set76, and Greenland ice-sheet runoff was taken from ref. 77. Tides were not taken into account in this simulation. Here we studied the model data of 2016 to 2018 in Fram Strait for comparison with our observations.

1-dimensional mixed layer depth model

The PWP78 1-dimensional mixed layer model simulates the response of the ocean to surface fluxes. It ignores horizontal gradients and horizontal advection. This allows to judge whether certain surface flux conditions can on their own explain observed conditions. We ran the PWP model (as implemented for Matlab by http://www.po.gso.uri.edu/rafos/research/pwp/) with four different scenarios (Fig. S6: P17-M17, P17-M18, P18-M17, P18-M18) where: P17: An idealized initial profile based on the observed profiles (Fig. 3) representing the conditions in 2017: constant temperature of 2 °C in the vertical, linear salinity gradient from 30.5 at the surface to 35 at 50 m and another linear salinity gradient from 35 at 50 m to 35.1 at 200 m. P18: An idealized initial profile based on the observed profiles (Fig. 3) representing the conditions in 2018: Same as P17 except that the surface to 50 m salinity gradient is from 34.8 to 35. M17: A time series of the the meteorological forcing (10 m wind velocity, heat fluxes, and evaporation minus precipitation) from the ERA-Interim reanalysis (Fig. 4b) at the grid point closest to mooring HG-IV for the period 15-May-2017 to 01-Aug-2017. M18: Same as M17 but for the period 15-May-2018 to 01-Aug-2018. M17 and M18 are provided in Supplementary Data 1.

Shipboard CTD data

Shipboard CTD casts of a standard dual sensor Seabird 911+ CTD-rosette were occupied in spatial and temporal vicinity to the moored observations (Tab. S2) on three cruises: PS107 in 2017 (https://doi.org/10.1594/PANGAEA.894189), PS114 in 2018 (https://doi.org/10.1594/PANGAEA.898694) of RV Polarstern, and JR17005 in 2018 (https://doi.org/10.5285/84988765-5fc2-5bba-e053-6c86abc05d53) of RRS James Clark Ross. The data were processed according to standard routine79. Additionally, we use underway CTD data from an OceanScience underway CTD collected during PS107 in 2017 (https://doi.org/10.1594/PANGAEA.886146) and processed according to ref. 21.

Mooring data

The mooring data discussed in this paper is from two mooring clusters in the central and eastern Fram Strait (named “HG-IV” at ~79°N 4°20’E and “F4” at ~79°N 7°E) where moorings were located as close to each other as possible (the horizontal separation was equal to the water depth) in order to enable more measurements than could be fit physically onto a single mooring. Tab. S2/S3 list the deployment and recovery details of the moorings including the exact latitudes/longitudes as well as the individual instruments on the moorings. Note that all data shown in this paper from ~30 m depth and the temperature/salinity/oxygen data from ~55 m is from the HG-IV-S-* and F4-S-* moorings, while all other data is from the HG-IV-FEVI-* and F4-* moorings. The AZFP data is from F5-17 located roughly half way between the two clusters. All sensor based mooring raw data (except for the ASL AZFP data) is available at ref. 80.

It is known that conversion factors for biogeochemical sensors (e.g., chlorophyll fluorescence) change over the seasons, depths, and regions81,82. In order to make as few assumptions as possible, we used the following approach: we could have determined the conversion factors from the instance when the ship was there with the CTD-rosette, but these conversion factors might not be appropriate for the majority of the time series. Hence, simply using the manufacturers’ calibrations, as we do here, introduces fewer uncertainties. Where we have different estimates of the same parameter, we present them together and demonstrate that they agree qualitatively and also mostly quantitatively (e.g., Fig. 5b). In particular the timing of events is robust.

At some locations, the target variables were not measured the whole time or the measurements failed, hence we present what is available. The vertical location of the instruments (Fig. 4c and S1c) varied substantially (intermittently up to 200 m) as a result of mooring blow downs caused by strong intermittent ocean currents. Time series have not been corrected for this vertical motion, but data are not used during blow downs in order not to bias the time series interpretation by temporal changes introduced by instruments traversing through vertical property gradients.

Physical sensor measurements

The physical sensors (for pressure, temperature, conductivity, and oxygen) were pre-cruise manufacturer calibrated and processed similar to ref. 83; the processed data is also available at ref. 80.

Mixed layer depth (MLD)

Since there are no autonomous vertically profiling measurements available, we can only determine the minimum value of the mixed layer depth. At each hourly time step, the potential density difference (Δσ) between the uppermost (~30 m) temperature/salinity recorder and the underlying temperature/salinity recorders is calculated. The 0.5th percentile of each Δσ time series is added to the Δσ time series for the different deployments. This fixes slight offsets in the temperature and/or conductivity calibrations which result in too negative or too positive density differences. The minimum estimate of the mixed layer depth at hourly resolution is then determined as the depth of the deepest instrument where Δσ < 0.05 kg m−3. If Δσ > 0.05 kg m−3 for all depths at a time step, then the minimum mixed layer depth can only be determined as 0 for that time step. Daily values of the MLD were defined as the depth at which three hourly realizations of MLD were shallower within a 24 h time span and at which the remaining 21 MLD realizations were deeper. This biases the daily MLD estimate towards situations where phytoplankton is kept in the surface ocean rather than also being mixed down for some amount of time.

Stratification estimated between 30 m and 55 m

Based on the temperature and salinity time series observed at ~30 m and ~55 m, we estimate the buoyancy frequency as ({N}^{2}=frac{-g}{{rho }_{0}}frac{Delta rho }{Delta z}) where g is the acceleration due to gravity, Δσ is the potential density difference over the vertical distance of Δz = 25 m, and ρ0 is the average density. The contributions to stratification due to temperature (N2T) and salinity (N2S) are estimated as ({N}_{T}^{2}=g*alpha frac{Delta T}{Delta z}) and ({N}_{S}^{2}=-g*beta frac{Delta S}{Delta z}), respectively, where ΔTS are the temperature/salinity differences and α/β are the thermal expansion/haline contraction coefficients estimated from the average temperature/salinity at the two measurement depths.

Apparent oxygen utilization (AOU)

Oxygen concentration from the microcats was calculated using the pre-cruise manufacturer calibrations. AOU was calculated as the atmospherically equilibrated oxygen concentration (calculated from measured pressure, temperature, and salinity with sw_satO2 from the Seawater toolbox available at http://www.cmar.csiro.au/datacentre/ext_docs/seawater.htm) minus the measured oxygen concentration.

Light

Polar night/polar day

The length of day (hours per 24 h that the sun is above the horizon) was calculated from the sunrise equation as implemented for Matlab by https://de.mathworks.com/matlabcentral/fileexchange/55509-sunrise-sunset.

Photosynthetically available radiation (PAR)

The WetLabs Eco PAR measured PAR for 5 (in 2016–2017) or 10 (in 2017–2018) individual measurements 1 s apart from each other before it slept for 1 h before repeating the measurement cycle. These 5 or 10 individual measurements are averaged linearly to obtain hourly values at ~30 m depth (Fig. 5a blue). Values below the detection limit are set to a constant of 10−1.32 μmol m−2 s−1. Hourly values are linearly averaged to daily values (Fig. 5a black). The incoming solar shortwave radiation varies as a function of season and latitude as well as cloud cover as represented in the ERA-Interim reanalysis (parameters ssr). Its unit of W m−2 is converted to PAR assuming a constant spectral distribution as 1 W m−2 = 2.1 μmol m−2 s−184. In order to compare the PAR measured at a depth of approximately 30 m to the surface values, we approximate a spectrally averaged diffuse attenuation coefficient for PAR in clear water using the values of85 as kd = 0.02 m−1 and apply it to calculate a constant exponential extinction applied to the reanalysis surface values (Fig. 5a yellow). The average PAR available (PARavailable) to phytoplankton being moved around in the clear water mixed layer of depth MLD was calculated as the depth averaged vertical integral of the clear water extinguished PAR at the surface (PARsurf from the shortwave radiation of ERA-Interim): ({{PAR}}_{{available}}=frac{1}{{MLD}}*{int }_{z=0}^{z={MLD}}{{PAR}}_{{surf}}*{e}^{-{k}_{d}z}{dz}) (Fig. 5a red).

Chlorophyll concentration and optical backscattering

Chlorophyll fluorescence

The WetLabs ECO Triplet measures fluorescence at a “chlorophyll wavelength” and at a “CDOM wavelength” as well as optical scattering at 700 nm. The conversion from fluorescence to chlorophyll a concentration (in μg l−1) follows a manufacturer determined conversion determined for a mono-culture of phytoplankton (Thalassiosira weissflogii), which typically overestimates the chlorophyll concentration. Hence, we applied the community-established calibration bias of 2 for the WetLabs ECO-series fluorometer to these in situ fluorometric chlorophyll values81. This conversion factor may be different in ocean waters of Fram Strait, but it still gives reasonable agreement with independent estimates.

Optical backscattering

The EcoTriplet measured 8 individual measurements 1 s apart from each other before it slept for 1 h before repeating the measurement cycle. For the chlorophyll fluorescence, the individual measurements are averaged to hourly values. For the scattering, times when individual 1-second measurements exceed 0.002 m−1 sr−1 are indicative of strong optical backscattering not due to small particles in the water column, but rather to larger potentially aggregated particles. The times of strong backscattering are marked individually (Fig. 5b red).

Nutrients

Nitrate (SUNA sensor)

Prior to deployment (11 and 15 days for sensors deployed at HG-IV and F4, respectively), the reference spectrum of the sensors were updated as per manufacturer specifications. We first let the sensors cool down for 24 h at 0 °C in a temperature controlled laboratory. Next, the reference spectrum update was achieved by measuring Milli-Q water (i.e., no nitrate present). To verify if this update was successful, solutions with three different nitrate concentrations (3, 7, and 14 μmol l−1) were then measured, with the output being monitored live (expected to be within ±2 μmol l−1 of each concentration). A measuring time of 20 s yields stable results and was thus applied during the deployments with an interval of 6 h. Upon recovery, SUNA data were processed using the SeaBird UCI software package version 1.2.1. Here, temperature and salinity data were used to remove the spectrum of bromide and compensate for temperature dependent absorption using an algorithm developed by ref. 86. This step yields the spectrum of nitrate only, at a precision of ±0.3 μmol l−1. The sensor is characterized by a drift of 0.3 μmol l−1 per hour lamp time. Given the deployment settings, a total operational time of about 8 h was accumulated. Therefore, a linear drift correction of 2.4 μmol l−1 (365 days)−1 was applied. Up to this point, however, accuracy remains at 2 μmol l−1 as per manufacturer specifications. Therefore, an offset correction is then applied based on the in situ concentrations observed at the beginning of the deployment as well as with the RAS (see below) where available, with outliers excluded.

Inorganic nutrients from Remote Access Samplers (RAS)

McLane RAS were programmed to draw two 500 ml samples (1 h apart, starting at noon) approximately every other week. Samples within the RAS were collected in sterile plastic bags and fixed with 700 μl of 50% mercuric chloride solution. Upon recovery, two samples from a given sampling date were combined to yield a volume of 1 l, required for bacterial and phytoplankton genetic analyses (see below), and a 50-ml aliquot destined for the measurement of dissolved inorganic nutrients. Aliquots for nutrient analysis were collected in PE bottles, which were then stored frozen (−20 °C) until analysis on land. Analyses for inorganic nutrients were carried out using a QuAAtro Seal Analytical segmented continuous flow autoanalyser following standard colorimetric techniques. The accuracy of the analysis was evaluated through the measurement of KANSO LTD Japan Certified Reference Materials and corrections were applied accordingly. Finally, we evaluated pressure, temperature, and salinity data from the CTD (SBE37-SMP-ODO) attached to the RAS to determine whether the two samples taken one hour apart on a given date drew water from the same depth and with consistent properties.

Carbonate system

pCO
2 and pH

The calibration of SAMI pH and SAMI CO2 sensors was carried out by the manufacturer, approximately 2 months prior to deployment. The calibration certificates specify accuracy and precision of ±0.003/±0.001 pH units and ±3/<1 μatm, respectively. pCO2 sensors measure and yield pCO2 in μatm. For the pH sensor, raw absorption data were converted to pH (total hydrogen ion scale) in combination with temperature and salinity (SBE37-SMP-ODO) using the quality control tool (QC_pH) supplied by the manufacturer.

Upon assessment of pCO2 data, values from the HG-IV mooring in 2017–2018 were deemed to be biased high by approximately 130 μatm (with a step jump at the turn-around of the moorings), therefore a constant value of this magnitude was subtracted from that record. The pH sensors ran out of battery towards the end of the deployments, resulting in interrupted records. At HG-IV, some erratic data before the sensors stopped recording were excluded for the first deployment after 16-Feb-2017 and for the second deployment after 02-Aug-2018. At F4 in 2017–2018 pH values below 8 were excluded.

Carbon takeup is estimated from the change in dissolved inorganic carbon (DIC) between the beginning of the bloom and the time when the minimum pCO2 is reached. In turn, DIC is calculated from pCO2 and alkalinity87 as well as measured temperature, measured salinity, phosphate concentration (0.5 μmol l−1), and silicate concentration (5 μmol l−1). Alkalinity is taken from the relationship Alk = 736 + 45.2 * S in ref. 88.

Microbial communities

Relative abundances of bacteria and microbial eukaryotes based on 16 S and 18 S rRNA gene sequences

The methodology followed89 which we briefly summarize here: the same water samples as for the inorganic nutrients (see above) were used, in which mercuric chloride resulted in fixation of microbes for long-term preservation90. After recovery, the ~1 l per sampling event was immediately filtered through 0.22 µm Sterivex filter cartridges (Millipore, Burlington, MA). Filters were immediately frozen at −20 °C until DNA extraction in the home laboratory.

DNA was extracted using the DNeasy PowerWater kit (Qiagen, Germany) according to the manufacturer’s instructions and quantified using a Quantus fluorometer (Promega, Madison, WI). Obtained DNA quantities ranged between 0.01 and 11 ng (µl)−1. Bacterial 16 S and eukaryotic 18 S rRNA gene fragments were amplified using primers 515F–926 R91 and 528iF–964iR92, respectively, according to the 16 S Metagenomic Sequencing Library Preparation protocol (Illumina, San Diego, CA). Amplicon gene libraries were sequenced using Illumina MiSeq instruments at CeBiTec (Bielefeld, Germany; bacteria) or AWI (eukaryotes). After primer removal using cutadapt93 reads were classified into amplicon sequence variants (ASVs) using DADA294.

After singleton removal, we obtained a mean of 60,000 bacterial and 119,000 eukaryotic reads per sample that sufficiently covered community composition89. Bacterial and eukaryotic reads were taxonomically classified using the Silva v138 and PR2 v4.12 databases, respectively.

Normalized mean volume backscattering (MVBS) as proxy for zooplankton biomass

MVBS

RDI Workhorse Longranger ADCPs were deployed, using a four-beam, convex configuration with a beam angle of 20° and frequencies of 76.8 kHz. The number of bins was set to 70 with a bin length of 8 m. The sampling interval was set to 20 pings per ensemble with a ping rate of about 20 pings every 60 min. The instruments were moored at a nominal depth of ~400 m in upward-looking mode and measured horizontal and vertical currents and acoustic backscatter intensity. Instrument heading, pitch and roll and temperature data were also collected. The echo intensities were given in an automatic gain control count scale of 0 to 255. Following95, they were converted to MVBS. We used beam-averaged data because the four beams together gave a better signal-to-noise ratio than individual beams. First, the noise level of all four ADCP beams was determined from the minimum values of RSSI (Received Signal Strength Indicator) counts obtained in the remotest depth cell, when the sea surface was outside the ADCP range. Sound velocity c and sound absorption coefficient a were considered variable with depth and time and calculated according to the UNESCO formulas from an interpolation in time of 2 CTD profiles collected at the beginning and the end of the deployment. The mean at each depth for the 20 single pings comprising an hourly burst was calculated.

An ASL Acoustic Zooplankton Fish Profiler (AZFP) was moored at F5 (79°N, 5°40′E) with the transducer faces pointing upward in ~150 m depth and operated at four different frequencies (38, 125, 200 and 455 kHz). The sampling interval was set to 30 s with a pulse duration of 0.5 ms (38 kHz) and 0.17 ms (125 and 200 kHz), respectively. Pitch and roll were measured with each ping. Data from the 455 kHz transducer were omitted from this analysis. AZFP data were post-processed and integrated with Echoview software version 11.0.239 (Echoview Software Pty Ltd, Hobart, Australia). Background noise was removed through time-varied thresholds for each transducer, and after removing noise and unwanted signals originating from e.g., other mooring devices that temporally became visible as backscatter in the echograms, the MVBS was integrated and exported for 24-hr bins.

For both ADCP and AZFP data the 50–100 m vertical average was calculated and for the AZFP also the 15–100 m and 30–100 m vertical averages. The median for each mooring deployment of the daily mean values of the vertical MVBS means was calculated. This median was subtracted from the MVBS to obtain the normalized MVBS, which corrects for possible hardware differences between the different deployments96.

Note that no AZFPs existed at moorings HG-IV and F4 and the ADCPs did not return data shallower than 50 m. Therefore, we use the comparison in Fig. S5 to show that the ADCPs deliver similar qualitative statements compared to the shallower reaching multi-frequency AZFPs. Apart from one exception, there are no large qualitative differences between the averaging in the different depth layers. That means that, most of the time, observations between 50–100 m do not miss a large part of the biomass. The exception occurred in August 2017 when the biomass above 30 m appeared to be much stronger. Likewise, there were no qualitative differences in the progression of the curves between the different frequencies. Therefore, the 50–100 m ADCP appears to be a reasonable proxy of zooplankton and fish biomass and we use it at HG-IV and F4 (Fig. 6c and S3c).

Particle and POC flux

Sediment volume flux in water column

Sediment traps were located 200 m and 1200 m below the sea surface at mooring HG-IV and at 200 m at mooring F4. The collector cup opening times ranged from 9 to 59 days depending on the season (lower resolution in winter). The cups had an interior diameter of 4 cm and the height of the sedimented layer on the bottom was measured from photos of the cups. From these the sediment volume was calculated and it was normalized by the 0.5 m2 collection area and the opening duration. We assume that this sediment volume flux is approximately proportional to total particulate matter (i.e., POC and lithogenic matter) flux and use it only to infer qualitative differences and the timing of events.

POC flux from sediment trap on bottom lander

Aliquots of the sedimented material collected at 2.5 m above the seafloor in sediment traps on the bottom lander at HG-IV were sieved through 500-µm mesh size to remove larger zooplankton swimmers or benthic organisms, then filtered on pre-combusted Whatman GF/F glass fiber filters, acidified with 0.1 N HCl, and dried at 60 °C. POC concentrations were determined with a CaloErba CN-analyzer. These were then normalized by the split factor as well as to the 0.25 m2 collection area and the cup opening duration. As the sediment trap was deployed 2.5 m above the bottom, it also collected resuspended material in addition to settling material.

POC flux inferred from benthic oxygen consumption

Benthic carbon mineralization was estimated from sediment oxygen consumption rates. The benthic crawler called TRAMPER97 was deployed at HG-IV for one year and moved 15 m along the seafloor every 7 days. Upon arriving at the new location, it measured oxygen concentration profiles through the top 15 cm of the sediment. From the shape of the oxygen profile, the benthic oxygen consumption, i.e., the oxygen flux from the water column into the sediment was calculated. Assuming a Respiratory Quotient of 1.0 (i.e., that O2 consumed via benthic diagenesis is balanced by a corresponding molar production of CO2) oxygen consumption rates were converted to POC fluxes reaching the benthos.

Seafloor imagery

Detritus seafloor areal coverage

An underwater camera (VTLC: Video Time Lapse Camera; AquaPix, USA) fitted to a benthic lander at HG-IV took 5-second video sequences twice a day for the 2017–2018 deployment year. Frame grabs were extracted and divided into 100 boxes referring to equal sized areas of the seafloor. At each time, the number of boxes covered by bare seafloor, by accumulating small white material, and green material was calculated. The white material, though very fine in size, accumulated over a larger area of seafloor than the much more massive green material identified as algal detritus. This is presumably related to hydrodynamics associated with the microtopography of the seafloor.

Number of megafauna present

All eelpout fish (Lycodes frigidus) and visible epibenthic megafauna including shrimps (Bythocaris spp), isopods (Saduria megalura), holothurians (Kolga hyalina), and gastropods (Mohnia spp.) observed in each frame grab were counted. Additionally, at the very bottom edge of the screen, a purple sea anemone and branches of the sponge Cladorhiza gelida were present, but not included in the data.


Source: Ecology - nature.com

Strategic Forest Reserves can protect biodiversity in the western United States and mitigate climate change

New visions for better transportation