Changes in gross oxygen production, net oxygen production, and air-water gas exchange during seasonal ice melt in Whycocomagh Bay, a Canadian estuary in the Bras d'Or Lake system

Sea ice is an important control on gas exchange and primary production in polar regions. We measured net oxygen production (NOP) and gross oxygen production (GOP) using near-continuous measurements of the O2/Ar gas ratio and discrete measurements of the triple isotopic composition of O2, during the transition from ice-covered to ice-free conditions, in Whycocomagh Bay, an estuary in the Bras d’Or Lake system in Nova Scotia, Canada. The volumetric gross oxygen production was 5.4+2.8 −1.6 mmol O2 m−3 d−1, similar at the beginning and end of the time series, and likely peaked at the end of the ice melt period. Net oxygen production displayed more temporal variability and the system was on average net autotrophic during ice melt and net heterotrophic following the ice melt. We performed the first field-based dual tracer release experiment in icecovered water to quantify air–water gas exchange. The gas transfer velocity at > 90 % ice cover was 6 % of the rate for nearly ice-free conditions. Published studies have shown a wide range of results for gas transfer velocity in the presence of ice, and this study indicates that gas transfer through ice is much slower than the rate of gas transfer through open water. The results also indicate that both primary producers and heterotrophs are active in Whycocomagh Bay during spring while it is covered in ice.


Introduction
The annual cycle of sea ice formation and melt regulates primary production and CO 2 uptake and ventilation in polar regions. Ice alters the rate of air-water gas exchange, reduces the penetration of light into surface water, changes stratification and mixing processes, and harbors microbes and biogenic gases including CO 2 (Cota, 1985;Loose et al., 2011a;Loose and Schlosser, 2011).
The question of whether climate change will increase or decrease Arctic Ocean carbon uptake is a topic of considerable debate (Bates et al., 2006;Bates and Mathis, 2009;Cai, 2011). Global warming is causing dramatic reductions in sea ice cover and increases in freshwater inflow and organic carbon supply to the Arctic Ocean, which impacts ecosystems (ACIA, 2004;Vaughan et al., 2013;Macdonald et al., 2015). Because conducting field work in the Arctic is chal-lenging, measurements of productivity and gas exchange are limited. Biogeochemical time-series observations resolving seasonal changes in productivity are particularly scarce in the Arctic (MacGilchrist et al., 2014;Stanley et al., 2015). Measurements at Palmer Station in Antarctica show a strong seasonality in biological productivity and carbon uptake associated with changes in light, physical mixing, and grazing, and demonstrate the benefits of high-frequency sampling for quantifying CO 2 uptake in seasonally ice-covered waters (Ducklow et al., 2013;Tortell et al., 2014;Goldman et al., 2015).
The parameterization of gas exchange in the presence of ice also remains highly uncertain. Many investigators have assumed that there is negligible gas transfer through ice, and therefore the gas transfer velocity can be linearly scaled as a function of the fraction of open water, multiplied by the open water gas transfer velocity Legge et al., 2015;Evans et al., 2015;Stanley et al., 2015). A recent field study by Butterworth and Miller (2016) in the Southern Ocean at 0 %-100 % ice cover verified this approach. However, other studies report that gas exchange is reduced or enhanced in the presence of sea ice relative to a linear scaling based on the fraction of open water, including some studies measuring higher transfer velocities in ice-covered waters than in open water (Fanning and Torres, 1991;Else et al., 2011;Rutgers van der Loeff et al., 2014). Additional studies show that gas exchange in ice-covered waters cannot be predicted from wind speed alone, which may be a cause of the wide range of results (Loose et al., 2009(Loose et al., , 2016Lovely et al., 2015). An additional factor reducing air-water exchange in ice-covered waters is that ice significantly reduces fetch for wave generation, and therefore wind-driven near-surface turbulence (Squire et al., 1995;Overeem et al., 2011). However, other processes may enhance near-surface turbulence in the presence of sea ice including convection associated with heat loss and brine rejection (Morison et al., 1992;Smith and Morison, 1993), boundary layer shear between ice and water (McPhee, 1992;Saucier et al., 2004), and wave interactions with drifting ice (Kohout and Meylan, 2008).
In this study, we measured productivity and gas exchange over a 1-month period during and following ice melt in Whycocomagh Bay, an estuary in the Bras d'Or Lake system on Cape Breton Island in Nova Scotia, Canada (see Sect. 1.1 for a full site description and map). We performed two dual tracer release experiments at different fractions of ice cover to quantify air-water gas exchange by adding 3 He and SF 6 to the mixed layer. We measured net oxygen production (NOP) and gross oxygen production (GOP) at Little Narrows, while Whycocomagh Bay transitioned from completely ice-covered to ice-free conditions. GOP is the total amount of O 2 generated by autotrophic microbes as a result of photosynthesis. NOP is GOP minus respiratory consumption of O 2 by autotrophs and heterotrophs. The NOP/GOP ratio provides an estimate of the export efficiency, i.e., the fraction of autotrophic production available for export to waters below the mixed layer. The time-series approach allowed us to quantify non-steady-state O 2 fluxes, which can be a significant fraction of the total O 2 flux in many regions Tortell et al., 2014;Manning et al., 2017b). We quantify GOP with discrete measurements of the triple oxygen isotopic composition of O 2 (Luz et al., 1999;Luz and Barkan, 2000;Juranek and Quay, 2013) and NOP with nearcontinuous measurements of the O 2 /Ar saturation anomaly (Cassar et al., 2009).
To our knowledge, there are no other published field-based experiments where the dual tracer technique was used in the presence of ice, and this study adds to a limited number of in situ measurements of NOP and GOP during ice melt (Goldman et al., 2015;Stanley et al., 2015;Eveleth et al., 2017).

Site description
The Bras d'Or Lake system consists of several interconnected channels and basins and has a total surface area of 1070 km 2 and an average depth of ∼ 30 m (maximum 280 m) (Petrie and Raymond, 2002;Petrie and Bugden, 2002). The Bras d'Or Lake system exchanges water with the Atlantic Ocean primarily through the Great Bras d'Or Channel at the northeastern region of the estuary; this channel has a shallow (16 m deep) and narrow (0.3 km wide) sill at the mouth (Petrie and Raymond, 2002). We conducted the work for this study in Whycocomagh Bay, an enclosed embayment approximately 13 km long and 3 km wide, at the northwestern end of the estuary, approximately 60 km from the open ocean ( Fig. 1). Whycocomagh Bay is separated from the rest of the Bras d'Or Lake system by Little Narrows, a channel which is approximately 0.2 km wide and 0.5 km long. In Whycocomagh Bay, there are two basins (western and eastern basin) that are up to 40 m deep. Little Narrows channel is ∼ 15-20 m deep . The annual maximum ice cover is typically reached in early March. Ice disappears rapidly during April until its total melt by the first week of May (Petrie and Bugden, 2002).
The Bras d'Or Lakes system exhibits standard estuarine circulation with freshwater flowing outward at the surface and denser saltwater flowing inward in the deeper layers . The deep water in western Whycocomagh Bay is strongly isolated due to the presence of a sill at the bottom of Little Narrows channel, and very little mixing occurs between the surface and sub-surface waters (Petrie and Bugden, 2002). Whycocomagh Bay receives a relatively large amount of freshwater compared to other regions of the bay. This freshwater forms ice in winter and the stratification remains very stable, preventing vertical mixing (Krauel, 1975). In general, excluding the narrow channels in direct contact with the ocean, the Bras d'Or Lakes system does not exhibit significant tidally driven variability in temperature, salinity, and sea surface level (Krauel, 1975;Petrie, 1999;Petrie and Bugden, 2002;  Dupont et al., 2003). A 21 d time series in western Whycocomagh Bay showed no detectable diurnal or semidiurnal tides (Dupont et al., 2003).
To our knowledge, there are no previously published measurements of water-column chemistry in Whycocomagh Bay in ice-covered conditions. Measurements from July 1974 showed that the western portion of Whycocomagh Bay became anoxic in the isolated deep waters below 25 m depth, whereas measurements in the eastern portion of the basin (closer to Little Narrows) showed the water column had an O 2 saturation of 61 % at the bottom depth of 30 m (Krauel, 1975;. Measurements collected from 1995 to 1997 (from late April to late September) showed O 2 concentrations in Whycocomagh Bay from 1 to 5 m depth were near equilibrium (94 % to 112 % saturation) throughout the bay. In the deeper waters, O 2 concentrations in eastern Whycocomagh Bay ranged from 69 % at 25 m depth on 28 April 1996 to 54 % at 13 m depth on 26 Septem-ber 1995 and 30 % at 30 m depth on 23 September 1996, and the western basin was persistently anoxic (Strain et al., 2001).
2 Field and analytical methods 2.1 Sampling setup at Little Narrows Table 1 presents key events that occurred during the experiment. We continuously sampled water at Little Narrows channel (Fig. 1) during a 33 d time series (25 March-28 April 2013). We moored a Goulds SB Bruiser 5-18 GPM submersible pump with intake at ∼ 0.5 m depth (within the mixed layer), placed inside a mesh filter bag to prevent large particles from clogging the pump, and a conductivity and temperature (CT) sensor (Sea-Bird Electronics SBE37) at ∼ 0.5 m depth. From the submersible pump, water flowed through flexible high-pressure PVC tubing submerged un-derwater to a 3-port pressure-relief valve (on shore) and was then pumped along shore (∼ 50 m) to our sampling apparatus. The water passed through three 10 canister filters (100, 20, and 5 µm nominal pore size) and then into a sampling manifold containing valves for distributing water for measurement of O 2 /Ar (continuously, by mass spectrometry) and for discrete sampling. As discussed below, we sampled discretely for SF 6 , 3 He, and O 2 /Ar and the triple oxygen isotopic composition of O 2 , and near-continuously for O 2 /Ar. Excess water flowed through tubing back into the bay. We covered the tubing on shore and the filter canisters in foam insulation to minimize temperature changes in the water.
We deployed a Nortek acoustic Doppler current profiler (ADCP) at ∼ 4 m depth in the middle of the channel from 7 to 28 April. The mean current speed at 0.5 m depth was 3.4 km d −1 (3.9 cm s −1 ) with an orientation toward 31 • (approximately along the channel axis), indicating the transit time through Little Narrows is relatively short. The ADCP data are shown in Supplement Fig. S1. Our measurements (salinity, O 2 /Ar, etc.) did not display any correlation with tidal cycles, consistent with previous studies of the Bras d'Or Lakes system indicating tides have a negligible influence on water properties within Whycocomagh Bay (Krauel, 1975;Petrie, 1999;Petrie and Bugden, 2002;Dupont et al., 2003).
The CT sensor was initially placed closer to shore than the water pump because the cable was not long enough to reach the pump, but after obtaining a longer cable, we were able to co-locate the CT sensor with the water pump (beginning 12 April). For the discrete samples, we used the CT sensor temperature and salinity measurements beginning on 12 April (when we moved the CT to the same location as the pump) to calculate the equilibrium concentration of each gas. Prior to 12 April, we collected measurements with a YSI temperature and salinity probe from the water on shore and used these measurements as the temperature and salinity for the discrete samples. We determined the average warming through the underway line to be 0.37 ± 0.22 • C based on comparisons between the temperatures from the CT sensor (in situ) and the YSI probe (on shore) after 12 April and applied this offset to all YSI temperature measurements. For the continuous O 2 /Ar data we used a temperature record from a thermocouple located in the sampling bucket because it had fewer gaps in time. We calibrated the thermocouple using an Aanderaa 4330 optode sensor which contains a temperature sensor (accuracy ±0.03 • C) and then decreased the temperature by 0.37 • C to correct for warming effects.
We installed the gas chromatograph (for measurement of SF 6 ) and mass spectrometer (for measurement of O 2 /Ar) inside a garage connected to the Little Narrows Ferry building. The majority of the wet equipment was set up outside the garage adjacent to a window on the garage, which was used for connecting equipment and power cables between the outdoors and indoors.
We deployed the water pump and in situ CT sensor to the east (oceanward) of the Little Narrows cable ferry which periodically crosses the channel and operates 24 h per day. We found no correlations between our temperature, salinity, and geochemical measurements and the position of the ferry within the channel. We collected conductivity, temperature, and depth (CTD) profiles with a SBE 19Plus sensor at Little Narrows, usually by boat using a winch, but occasionally by lowering the CTD by hand on a rope from the Little Narrows cable ferry. This CTD package was also equipped with an O 2 sensor, but unfortunately the O 2 sensor malfunctioned throughout the experiment. Therefore, we can characterize vertical structure of salinity and temperature but not O 2 .
The GPS-equipped boat enabled us to map out the location of the ice edge nearest to Little Narrows, to perform the second tracer injection, and to sample after the tracer injection.

Tracer injections
The approach in this study was to dissolve the tracer mixture ( 3 He/SF 6 ) in Whycocomagh Bay, and sample continuously at Little Narrows, a constriction at the mouth of the bay. The net surface flow within Whycocomagh Bay, Little Narrows, and St. Patrick's Channel is toward the ocean due to the substantial freshwater inputs to the bay from surface runoff and precipitation Yang et al., 2007), and therefore tracer dissolved within the bay at the surface will eventually pass through Little Narrows, or be ventilated to the atmosphere. Two tracer injections occurred during the time series, resulting in estimates of the gas transfer velocity for two extremes: near-complete ice cover and essentially ice-free conditions. Injection 1 occurred through a hole in the ice from 30 to 31 March, near MacInnis Island (Fig. 1a). Approximately 0.11 mol SF 6 and 4.0 × 10 −4 mol 3 He were diluted by a factor of 50 with N 2 and then bubbled using a manifold within the mixed layer, over a 21 h period. The ice thickness was ∼ 0.3 m and the injection occurred in the upper 0.5 m of the water column. Because the tracer was added at a fixed location, it was added very slowly to increase the fraction of the gas that dissolved. We sampled for initial 3 He and SF 6 concentrations at the injection site, immediately before starting the injection, and after terminating the tracer addition, by drilling several holes > 10 m away from the injection site (because the bubbling action would have perturbed gas concentrations at the injection site). Subsequently, we sampled the tracer as it flowed through Little Narrows from 7 to 11 April. From 31 March to 11 April, the bay was nearly completely full of ice, with a small opening near Little Narrows (Fig. 1a). We also collected under-ice samples for O 2 /Ar and O 2 triple oxygen isotope (TOI) composition immediately before and after the injection but they displayed a wide range in values for O 2 /Ar, from −14 % to 2 %. Therefore, we were unable to define any initial under-ice values for these parameters.
Injection 2 occurred by boat on the morning of 19 April. By this time, the bay was nearly ice-free and all tracer from End of continuous water sampling and ADCP measurements at Little Narrows, and last CTD profile at Little Narrows the previous experiment had passed through and/or been ventilated to the atmosphere such that the tracer concentrations at Little Narrows were below detection. While the boat was moving, we used the same manifold as for Injection 1 to bubble approximately 4.4 mol SF 6 and 0.021 mol 3 He, diluted by a factor of four with N 2 into the mixed layer (Fig. 1a). The injection lasted 40 min. We detected the tracers at Little Narrows beginning on 20 April 23:30 UTC and measured the change in the ratio between 20 and 23 April as the tracer patch flowed through Little Narrows. We used a higher quantity and concentration of SF 6 and 3 He during this injection because we expected the tracer would be ventilated more rapidly due to higher gas exchange rates given the open water conditions. We were able to inject the tracer more rapidly because it was distributed over a large area instead of through a small hole in the ice.

Measurement of O 2 /Ar and the triple oxygen isotopic composition of O 2
We set up an equilibrator inlet mass spectrometer (EIMS) for measurement of O 2 /Ar similarly to the system described in Cassar et al. (2009). However, we used a larger membrane contactor cartridge (Liqui-Cel MiniModule 1.7 × 5.5) because the design is more robust than the Liqui-Cel Micro-Module 0.75×1 used by Cassar et al. (2009). The water flow rate through the cartridge was ∼ 1.5 L min −1 . We attached a custom female Luer-Lok fitting paired to a capillary adapter to the upper headspace sampling port and the lower sampling port was left closed.
For O 2 /Ar and the triple oxygen isotopic composition of O 2 , we collected samples in pre-evacuated, pre-poisoned glass flasks from a spigot in the water pumped to shore, or for shipboard sampling, using a small submersible water pump. Analysis occurred within ∼ 6 months of flask evacuation and 4 months of sample collection at the Woods Hole Oceanographic Institution with a Thermo Fisher Scientific MAT 253 isotope ratio mass spectrometer, following the method of Barkan and Luz (2003) with modifications as described in Stanley et al. (2010Stanley et al. ( , 2015.
The precision of the discrete samples, calculated based on the standard deviation of equilibrated water samples run daily along with the environmental samples was 0.011 and 0.020 ‰ for δ 17 O and δ 18 O, respectively, 5.6 per meg for 17 , and 0.07 % for O 2 /Ar. The mean difference between the EIMS and discrete samples was 0.05 %, and the mean magnitude of the difference was 0.35 %; given the small mean offset, the EIMS data were not calibrated using discrete samples.

Measurement of SF 6
For SF 6 , we collected ∼ 20 mL of water in 50 mL glass gastight syringes, then added ∼20 mL of nitrogen and allowed the samples to be shaken for 10 min to achieve equilibration between the headspace and water (Wanninkhof et al., 1987). After the water equilibrated to room temperature, we injected 1 mL of the headspace into an SRI-8610C gas chromatograph with an electron capture detector, heated to 300 • C (Lovely et al., 2015). We calibrated the detector response using a 150 parts per trillion SF 6 standard (balance N 2 ).
We also operated an automated gas extraction system at Little Narrows which sampled nearly every hour. The system recirculated a 118 mL loop of water through a membrane contactor, and the headspace of the membrane contactor was under vacuum, causing the dissolved gas to be extracted from the water into the headspace. When the automated system showed measurable SF 6 in the water, we began collecting discrete samples for SF 6 and 3 He. SF 6 equilibrium solubility concentrations are calculated following Bullister et al. (2002) and diffusivity is from King and Saltzman (1995). Precision of the system, based on the standard deviation of measurements of the 150 ppt standard, was 7 %. We assume a dry atmospheric mole fraction of 8 ppt for SF 6 (Bullister, 2015). For 3 He analysis, we collected samples in copper tubes, mounted in aluminum channels and sealed the samples at each end using clamps. Sample analysis occurred at the Lamont Doherty Earth Observatory, using a VG5400 mass spectrometer for 3 He and 4 He concentration (Ludin et al., 1998). Error for 3 He sample analysis (combined precision and accuracy) was ≤ 2 % of the measured 3 He-excess concentration above equilibrium. We used He solubility from Dempsey E. Lott III and William J. Jenkins (personal communication, 2015) and diffusivity from Jähne et al. (1987a). The Lott and Jenkins solubility values are ∼ 2 % higher than published data from Weiss (1971). The He solubility is for bulk He and we calculate the 3 He solubility using an atmospheric mole ratio M( 3 He)/M( 4 He) = 1.399×10 −6 (Mamyrin et al., 1970;Porcelli et al., 2002), although some more recent results indicate the current ratio may be slightly lower, 1.390 × 10 −6 (Brennwald et al., 2013). We use the He equilibrium isotopic fractionation as described in Benson and Krause (1980b).

Calculations, results, and discussion
The three goals of the experiment were to (1) quantify gas transfer velocity by dual tracer release, (2) quantify gross oxygen production from the triple isotopic composition of O 2 , and (3) quantify net oxygen production from O 2 /Ar. We begin by discussing the hydrographic characteristics of the study area and then describe the calculations, results, and interpretation for each of the three goals, in sequence.

Hydrography
We began sampling O 2 at Little Narrows on 25 March, when Whycocomagh Bay was nearly (> 95 %) fully covered by ice, and completed sampling on 28 April, when the bay was completely ice-free (Fig. 2). The surface ice cover retreated most rapidly between 16 and 18 April and was completely gone by 22 April or perhaps even earlier. Figure 2 shows 18 April and 23 April; MODIS satellite images on 22 April were also ice-free (but more blurry, and so are not shown in the figure) and we estimated ice cover to be < 10 % in the bay during surveys by boat during daytime on 20 April. The ice was likely melting even at the beginning of the time series since the surface water temperature was always above the freezing temperature of water (Fig. 3). Changes in surface ice cover and total ice volume are both important factors during the study; changes in ice volume/thickness will affect stratification and convection in the mixed layer as well as light penetration through the ice, and the surface ice cover affects the rate of gas exchange (Smith and Morison, 1993;Butterworth and Miller, 2016;Loose et al., 2016).
CTD profiles at Little Narrows channel near the water pump intake showed substantial changes in stratification dur-ing the time series (Fig. 3). From 25 March through 8 April the water column was strongly stratified and we estimated the mixed layer depth to average 0.8(0.3) m. During this period, it was often difficult to determine the exact mixed layer depth because the mixed layer depth was similar to the length of the CTD instrument and obtaining a stable CTD response so close to the surface was challenging. Following 8 April, the mixed layer deepened and warmed and its salinity increased, likely due to convection and heating following sea ice melt. For this period, we defined the mixed layer depth as the first depth where the density is 1 kg m −3 greater than the value at 1 m. The mixed layer reached a maximum of 3.0 m on 23 April and then shoaled by the end of the time series on 28 April. On most days, the density profile and mixed layer depth were driven by stratification in salinity, but for the final profile on 28 April the mixed layer depth was determined by a combination of temperature and salinity stratification due to heat uptake by the surface water following the ice melt. These changes in mixed layer depth must be considered in order to interpret the O 2 measurements and to quantify the gas transfer velocity.

Calculation
We calculate the gas transfer velocity using the dual tracer approach (Watson et al., 1991;Wanninkhof et al., 1993). For each experiment, we dissolved a mixture of 3 He and SF 6 in the water, both of which are normally present at very low ambient concentrations, and then measured the change in the ratio of the two gases as a function of time. Measuring two tracers enables correction for dilution and mixing, which reduces the excess concentrations of both gases (relative to air-water equilibrium) but does not change their ratio. Over time, the concentrations of both gases decay toward air-water equilibrium as gas is ventilated to the atmosphere through gas exchange. Because the molecular diffusivity of 3 He is 8-9 times higher than SF 6 , 3 He is lost to the atmosphere more rapidly than SF 6 , and therefore the 3 He:SF 6 ratio decreases with time. The ratio of the gas transfer velocity for the two gases is expressed as where k is the gas transfer velocity (m d −1 ) and Sc is the Schmidt number (unitless), defined as the kinematic viscosity of water divided by the molecular diffusivity of the gas in water, and n is an empirical exponent, typically between 0.5 and 0.67 (Jähne et al., 1984;Liss and Merlivat, 1986). Using a time series of measurements of the two gases, the gas transfer velocity for 3 He is calculated as (2) Using this equation and a cost function, we can find the value of k3 He that minimizes the error between the measured and modeled [ 3 He] exc / [SF 6 ] exc . Once k3 He is known, we can calculate k for any other gas using Eq. (1) by substituting Sc SF 6 for the Schmidt number of the gas of interest. For example, for Sc = 600 (the Schmidt number for CO 2 at 20 • C in freshwater), For this study, we use a Schmidt number dependence of n = 0.5 which is appropriate for wavy, unbroken water surfaces (no bubble entrainment) (Jähne et al., 1984(Jähne et al., , 1987aLiss and Merlivat, 1986;Ho et al., 2011b). At Little Narrows, we observed that tidal currents generated surface waves, even at low wind speeds. These waves produce near-surface water turbulence which is the ultimate driver of air-water gas exchange (Jähne et al., 1987b;Wanninkhof et al., 2009). The Schmidt number at a salinity of 4 PSS and temperature of 2 • C is 305 for 3 He and 2684 for SF 6 . The initial Sc SF 6 /Sc3 He ratio was 8.7 for Injection 1 and 8.2 for Injection 2.

Results
The gas transfer velocity was much lower for Injection 1, which was sampled while the basin was essentially full of ice (31 March-10 April), compared to Injection 2, which was sampled when the basin was nearly ice-free (20-23 April). We used Eq.
(3) to model the measurements. We started the model at the time of the first measurement, initialized it with an initial excess concentration ratio and ran it through time for the duration of the injection. We selected the value of k3 He yielding the smallest root mean square deviation (RMSD) be-   (22,26,29,and 31 March,and 4,6,8,12,16,20,23,and 28 April). The black dot-dashed line shows the estimated mixed layer depth. The data in this plot were binned into 0.2 m depth bins from 0.4 to 3.0 and 0.3 m depth bins from 3.3 to 9.9 m. The shallowest depth varies between casts due to challenges in getting stable CTD data in the upper 1 m. tween the measured ratio and modeled ratio for each injection. The model ran 1000 times using a Monte Carlo simulation where the measured excess concentration ratios, including the initial condition, are varied with a Gaussian distribution, with the standard deviation being the estimated measurement error in the ratio (7.3 %).
We assume a constant Sc3 He /Sc SF 6 and mixed layer depth (h) for each injection. In actuality, during each injection Sc3 He varies by ∼1 % and the ratio of the Schmidt numbers varies by less than 0.2 %, so this is a small source of error. For Injection 1, we assume a mixed layer depth (h) of 0.8(0.3) m. This depth is consistent with the salinity profiles at Little Narrows ( Fig. 3a) between 31 March and 8 April (between 0.6 and 1.0 m depth), as well as measurements with a handheld temperature probe at the site of Injection 1 which indicated that the mixed layer depth was between 0.75 and 1 m. The mixed layer depth may have increased slightly between 8 and 10 April, but was likely 1.3 m or less (Fig. 3).
For Injection 1, the excess SF 6 and 3 He concentrations were reduced by 2 orders of magnitude by the time the tracer reached Little Narrows (7-11 d after injection). The tracer ratio did not display a consistent decrease over the 3 days we sampled it at Little Narrows, which may be due to the substantially lower gas transfer velocity, as well as the very low tracer concentrations potentially increasing mea- . Measured and modeled excess tracer ratios. Measured and modeled ratio of excess 3 He/SF 6 , normalized to the initial measured ratio for each injection. The modeled excess ratio is calculated using the k3 He that minimizes error between the model and measurements. Model results are shown for the model initialized with the initial measured concentration (solid lines), and starting 1 standard deviation above or below the measured initial concentration (based on an error of 7.3 % in the tracer ratio, dashed lines). surement uncertainty. The best fit to all the measurements was k 600 = 0.0457(0.0051) m d −1 , in the presence of ice and a shallow mixed layer, with the uncertainty the standard deviation of the distribution of k 600 from the Monte Carlo simulation (Fig. 4). The mixed layer appeared to deepen between the CTD profiles on 8 April 16:08 UTC and 12 April 19:13 UTC (Fig. 3), and it is possible that the mixed layer depth on 9-10 April may have been slightly deeper than the estimate of 0.8(0.3) m. However, if mixed layer deepening had a significant influence on our gas transfer velocity estimate, we would expect k 600 (calculated assuming a constant mixed layer depth) to be lowest when calculated over the longest time period, using the sample collected on 10 April. Instead, the gas transfer velocity was actually the lowest when integrated to 9 April (the excess tracer ratio appears above the best-fit line) and second lowest on 8 April. Since the gas transfer velocities for Injection 1 integrate over 7-10 d, any change in mixed layer depth during the last 1-2 d will have a small effect on the calculated k.
For Injection 2, we use a mixed layer depth of 2.7(0.3) m based on CTD profiles at Little Narrows on 20 and 23 April, which had mixed layer depths of 2.4 and 3.0 m, respectively (Fig. 3a). Calculation of the gas transfer velocity for Injection 2 was relatively straightforward as the ratio of excess 3 He/SF 6 steadily decreased over the five measurements (Table 2). The best fit to all four measurements was k 600 = 0.71(0.13) m d −1 for open water (Fig. 4).
Using the published He solubility from Weiss (1971) instead of the unpublished data of Dempsey E. Lott III and William J. Jenkins (personal communication, 2015) results in a gas transfer velocity that is 8 % lower in the presence of ice (for Injection 1) and 0.4 % lower in open water (for Injection 2).

Discussion
The gas transfer velocity calculated for Injection 1 is the effective gas transfer velocity (k eff ); it averages the gas transfer velocity through ice (k ice ), weighted by the time the tracer spent under ice, and the gas transfer velocity for open water (k), weighted by the time the tracer spent in open water . In partially ice-covered waters, the effective gas transfer velocity is sometimes calculated as where f is the fraction of open water Lovely et al., 2015). If k ice is negligible, then k eff = (f )k (Loose and Jenkins, 2014;Crabeck et al., 2014;Butterworth and Miller, 2016). For Injection 2, we determined k, the value for open water. We expect k ice to be lower than k because the ice acts as a physical barrier to gas exchange. The rate of gas molecular diffusion in water (Jähne et al., 1987a;King and Saltzman, 1995) is higher than gas diffusion through ice (Gosink et al., 1976;Ahn et al., 2008;Loose et al., 2011b;Lovely et al., 2015). However, the exact rate of gas diffusion through saltwater ice (and by extension the value of k ice ) is not well constrained and likely varies based on the physical properties of the ice such as brine volume and temperature (Golden et al., 2007;Loose et al., 2011b;Zhou et al., 2013;Moreau et al., 2014;Lovely et al., 2015).
To evaluate these results within the framework of Eq. (5), we must estimate the fractional ice cover during Injection 1.
Visual surveys along the shoreline of Whycocomagh Bay and satellite data indicated that the bay was nearly fully covered with a continuous sheet of ice from 31 March to 10 April, except for an opening close to Little Narrows ( Fig. 2a-b). Beginning between 7 and 12 April, a small region of water appeared to open up along the shoreline northwest of the site of Injection 1; however, by this time the tracer patch had moved eastward close to Little Narrows and likely was not significantly affected by this open water (Fig. 2c). We mapped out the location of the ice edge closest to Little Narrows by boat on 26 March and 7 and 12 April (Fig. 1a). Using these surveys and shoreline data, we calculate that for the surface area of the bay between the injection site and Little Narrows, f = 0.01 on 26 March, 0.05 on 7 April, and 0.08 on 12 April. The f experienced by the tracer patch during Injection 1 is likely between 0.02 and 0.08 because the tracer was injected on 30-31 March and flowed through the open water near Little Narrows between 6 and 11 April. If k for open water is the same for both injections (k = 0.71 m d −1 ), then the results yield k eff = (f )k with f = 0.064, which is consistent with the fraction of open water we estimate for Injection 1. Thus we conclude that k ice was negligible, compared to k. For example, if k ice was even 10 % of k for open water, then k eff for Injection 1 would have been ∼ 0.11 m d −1 , more than double the observed value of 0.0457(0.0051) m d −1 . One source of uncertainty in estimating the correct value of f is that the transit speed between the injection site and Little Narrows was non-constant. The mean current velocity at Little Narrows channel was 3.4 km d −1 , but the tracer took approximately 8 d to flow 7 km from the injection site to Little Narrows.
In calculating GOP and NOP by oxygen mass balance, we apply the tracer-based gas transfer velocities estimated by dual tracer release throughout the time series, since there is no consensus on the best treatment of gas transfer in lakes and estuaries (Clark et al., 1995;Cole and Caraco, 1998;Crusius and Wanninkhof, 2003;Ho et al., 2011a), nor on the parameterization of gas transfer in the presence of ice (Else et al., 2011;Lovely et al., 2015;Butterworth and Miller, 2016). Additionally, if bottom-derived turbulence (e.g., from tidal flow) is a significant contributor to air-water gas exchange, a parameterization based on wind speed alone may not be appropriate. This method of calculating one average k 600 for each injection does not enable the development of a wind speed-dependent parameterization for the gas transfer velocity.
Because the k 600 values for Injection 1 and Injection 2 are very different, the treatment of the gas transfer velocity in between the two injections strongly affects the productivity estimates for this period. We use k 600 = 0.0457(0.0051) m d −1 from the beginning of the time series until the end of 15 April. Figure 2d, collected on 16 April, is the first satellite image showing substantial open water within Whycocomagh Bay, but the open water is primarily in the western half of the bay, far from Little Narrows. We use k 600 = 0.71(0.13) m d −1 from 21 April until the end of the time series on 28 April. Surveys by boat on 19 and 20 April indicated < 10 % ice cover on these days, and we collected the first tracer measurements following Injection 2 on 20 April 23:30 UTC. Between 16 and 20 April, we apply a linear interpolation of the k 600 for Injection 1 and Injection 2 as a function of time. The gas transfer velocity is most uncertain during the period when the ice cover rapidly decreased because we do not have any measurements of gas transfer at intermediate ice cover. However, because the ice cover retreated rapidly, only 4 days of the productivity estimates (out of a 33 d time series) are affected by the uncertainties in gas transfer at intermediate ice cover.

Comparison with other estimates
To compare the gas exchange estimates with other published studies, we use wind speed data measured at 10 m height (u 10 ) at Eskasoni Reserve, 27 km northeast of Little Narrows (Fig. 1b) and archived by the Government of Canada (http://climate.weather.gc.ca, last access: 27 August 2019). The archived data are 2 min averages measured once per hour. During the period that samples with tracer from Injection 2 were collected (between 20 April 23:00 and 23 April 11:00), the average wind speed was 2.6(1.4) m s −1 and the median was 2.2 m s −1 . The calculated k 600 over this time period from dual tracer data is 0.71(0.13) m d −1 . Cole and Caraco (1998) found k 600 = 0.636(0.029) m d −1 (95 % confidence interval) and their estimate is independent of wind speed in a lake with daily wind speeds of 1.39(0.06) m s −1 (95 % confidence interval); this k 600 is consistent within uncertainty with the Injection 2 result. Standard open ocean parameterizations that use a quadratic dependence on wind speed predict k 600 0.5-0.6 m d −1 with uncertainties of ∼ 20 % or ∼ 0.10 m d −1 (Ho et al., 2006;Sweeney et al., 2007;Wanninkhof, 2014). Crusius and Wanninkhof (2003) found that in a lake gas exchange can be estimated nearly equally well using three different parameterizations. Using their parameterizations with the wind speed record during our time series, we calculate transfer velocities of 0.32-0.66 m d −1 , and the velocity is most similar to our result when we use a constant gas transfer velocity k 600 = 0.24 m d −1 for u 10 < 3.7 m s −1 and k 600 = 1.23u 10 − 4.30 for u 10 ≥ 3.7 m s −1 . However, Crusius and Wanninkhof (2003) parameterized the gas transfer using instantaneous (e.g., 1 min averaged) winds measured throughout the time series, not once per hour, and emphasize the importance of including the variability in short-term winds when quantifying gas exchange at low wind speeds. If gas transfer velocity has a nonlinear dependence on wind speed, then short-term wind speed measurements will more accurately represent the gas transfer than wind speed values averaged over longer periods (Livingstone and Imboden, 1993;Crusius and Wanninkhof, 2003). Since we only have 2 min averages measured once per hour (for a total of 60 measurements during Injection 2), the wind record we use may not fully represent the variability in winds during the Injection 2 measurement period.
A source of error in comparisons with published results is that the wind speed data come from a different location than the study area. Although Eskasoni Reserve is adjacent to the Bras d'Or Lake, the local topography and bathymetry are different near the reserve and in Whycocomagh Bay. Thus, it is likely that the wind speed and momentum stress at the air-water interface differ at Whycocomagh Bay compared to Eskasoni Reserve (Ortiz-Suslow et al., 2015).
The measurements are in agreement with other studies showing that gas transfer velocity is significantly reduced under near-complete ice cover (Lovely et al., 2015;Butterworth and Miller, 2016) and contrast with studies showing enhanced gas transfer under > 85 % ice cover (Fanning and Torres, 1991;Else et al., 2011). We find that k eff = (f )k for > 90 % ice cover, but we cannot evaluate whether the same equation holds at intermediate ice cover because there was no injection at a lower fractional ice cover. In this study, the ice cover was near-continuous across the entire bay during Injection 1 and likely did not contain the ice leads that are common in Arctic and Antarctic sea ice; differences in gas transfer behaviour are expected based on the nature of the ice pack.

Calculation
The triple oxygen isotopic composition of O 2 is an effective tracer of gross photosynthetic O 2 production . Due to reactions in the upper atmosphere that impart a small mass-independent isotopic signature on atmospheric oxygen, dissolved O 2 derived from air-water exchange has a unique triple isotopic signature compared to O 2 generated by photosynthesis and O 2 consumed by respiration. We characterize the oxygen isotopic composition using We report 17 with λ = 0.5179, the ratio of the fractionation factors for respiratory O 2 consumption in 17 O relative to 18 O (i.e., λ= 17 / 18 , where is the isotopic fractionation of O 2 due to respiratory consumption) (Luz and Barkan, 2005). This value for λ is selected so that 17 is nearly unaffected by respiratory O 2 consumption and reflects the proportion of O 2 that is derived from photosynthesis relative to air-water gas exchange (Hendricks et al., 2005;Juranek and Quay, 2013;Nicholson et al., 2014). We report 17 in per meg (1 per meg = 0.001 ‰) due to the small range of values in natural waters, typically 8-242 per meg Manning et al., 2017a).
Two key constraints in the calculation of GOP from measurements of the triple isotopic composition of O 2 are the isotopic composition of O 2 derived from air-water exchange and the isotopic composition of photosynthetic O 2 . The composition of photosynthetic O 2 is dependent on the triple oxygen isotopic composition of H 2 O, the substrate for photosynthetic O 2 , and the isotopic fractionation associated with photosynthetic O 2 production. In oceanic studies, a common assumption is that the H 2 O isotopic composition is equivalent to VSMOW (standard mean ocean water), but in brackish systems it is necessary to estimate the isotopic composition of water and incorporate this into the GOP calculation (Manning et al., 2017a Manning et al. (2017a). We assume that the waters in the estuary represent a mixture of two endmembers: seawater and meteoric (precipitation-derived). We define the salinity and δ 18 O-H 2 O for the two endmembers and then calculate δ 18 O-H 2 O for each water sample collected during the time series as a linear function of salinity. A similar approach is applied for δ 17 O-H 2 O.
For the seawater endmember, we use compilations of δ 18 O-H 2 O and salinity Bigg and Rohling, 2000) available from an online database . We included all 19 near-surface samples (< 5 m depth) between 44-48 • N and 58-64 • W in the database.
with λ w = 0.528 and all isotopic compositions referenced to VSMOW-SLAP. Equations (7) and (8) (Luz and Barkan, 2010). The value of λ w = 0.528 is well established for meteoric waters and seawater (Meijer and Li, 1998;Landais et al., 2008;Luz and Barkan, 2010). Spatial variability in the 17 O-excess of natural waters is less well understood due to the currently limited observations at sufficient accuracy to resolve the small excess (Luz and Barkan, 2010;Li et al., 2015). To calculate the freshwater and seawater endmembers for δ 17 O-H 2 O, we use the average values of 17 O-excess of 33 per meg for meteoric water and −5 per meg for seawater (Luz and Barkan, 2010 In this equation, k O 2 is the gas transfer velocity for O 2 (m d −1 ), [O 2 ] is the O 2 concentration (mol m −3 ), h is the mixed layer depth (m), X 17 = r( 17 O/ 16 O), and λ = 0.5179 (Eq. 7). The subscripts eq and p refer to O 2 at air-water equilibrium and produced by photosynthesis and terms without a subscript ([O 2 ], X 17 , and 17 ) are the measured mixed layer values. The first term on the right-hand side of Eq. (9) is the steady-state GOP term, and the second term is the nonsteady-state GOP term. If the system is at steady state with respect to 17 (i.e., there is no change in 17 with time), then the second term on the right-hand side of Eq. (9) equals zero and can be eliminated.
This non-steady-state equation for GOP incorporates temporal variability in both O 2 concentration and isotopic composition (see Eq. 5 in Prokopenko et al., 2011). It does not account for time variability in mixed layer depth (e.g., entrainment of deeper waters into the mixed layer), nor does it account for vertical diffusion of O 2 across the mixed layer. In order to account for vertical processes in the O 2 mass balance, we would need to have triple oxygen isotope measurements below the mixed layer (Castro Morales et al., 2013;Munro et al., 2013;Wurgaft et al., 2013;Nicholson et al., 2014). Since we do not have these data, we cannot correct for vertical processes affecting the O 2 budget.
We calculate X 18 eq based on Krause (1980a, 1984) and X 17 eq using 17 eq = 8 per meg (Reuer et al., 2007;Stanley et al., 2010), which is consistent with the daily measurements of distilled water equilibrated at room temperature that were analyzed along with the environmental samples (8.1 per meg with a standard error of 1.6 per meg, n = 12), as well as prior measurements of distilled water equilibrated at < 5 • C (Rachel H. R. Stanley, unpublished data). We calculate X 18 p and X 17 p using the salinity-dependent isotopic composition of H 2 O defined above and isotopic fractionation factors for photosynthetic O 2 with respect to H 2 O based on data in Luz and Barkan (2011) for average phytoplankton (α 18 p = 1.0033890 and α 17 p = 1.0017781). The Matlab code used to calculate GOP and the triple oxygen isotopic composition of water (from two-endmember mixing of δ 18 O-H 2 O and salinity) is available online (Manning and Howard, 2017).
We calculate gross oxygen production using samples collected at Little Narrows from 25 March to 27 April (Fig. 6). Visual inspection of the 17 data indicated that 17 changed during the time series, and therefore the calculation includes a non-steady-state GOP term. The non-steady-state term in Eq. 9 is h[O 2 ]∂ 17 /∂t. To calculate the rate of change in 17 with time, we first averaged the data into 24 h bins (beginning and ending at 19:30, local sunset) to avoid overweighting times when samples were collected at higher frequency. We calculated the average 17 and sampling time for all samples collected each day. Next, we separated the data into two periods: one period began on 25 March and ended 19 April 07:30, and the second period covered the remainder of the time series (ending 27 April). A linear regression of 17 versus time for the two time periods yielded a slope of 0.67 per meg d −1 (r 2 = 0.47) for the first period and −2.99 per meg d −1 (r 2 = 0.94) for the second period. The approximate timing for the change between periods was determined by visual inspection and then adjusted to maximize the r 2 and so that the equations of the two lines gave very similar 17 values at 19 April 00:00 (within 1 per meg). Splitting the period from 25 March to 19 April into two separate regressions (or one period where 17 increased and one period where it was constant) yielded much lower r 2 values and a discontinuous 17 record (different 17 values at the end of one period and the start of another), so a single regression was used for this period.
The other two variables in the non-steady-state GOP term are the mixed layer depth ( (Cassar et al., 2011). Using Eq. (9), GOP is estimated for each sample, using an isotopic composition for H 2 O based on the salinity of the sample, and a mixed layer depth, rate of change in 17 with time, and gas transfer velocity based on the sampling time (values shown in Fig. 6), and then calculate the daily average GOP from all samples on a given day (beginning and ending at 19:30, local sunset). The mean number of samples per day was two, the maximum was four, and a few days had no measurements. The uncertainty in GOP on each day is calculated by propagating uncertainty in k O 2 (11 % for Injection 1, 18 % for Injection 2), uncertainty in the mixed layer depth (from 10 % to 38 %, 0.3 m), uncertainty in the rate of change in 17 with time (22 % and 9 % where 17 is increasing and decreasing, respectively), and uncertainty in the photosynthetic endmember (discussed below). Measurement uncertainty in the isotopic composition of O 2 is excluded from the error calculation because it is a random error, rather than a systematic error (meaning that by taking many measurements of 17 over several days, the measurements with high and low 17 will average out) and because the measurement error is smaller than most other sources of error. All uncertainties are expressed as the standard deviation.
The GOP calculated with the local isotopic composition of H 2 O is 46 %-97 % higher (mean 74 % higher) than GOP calculated assuming the water's isotopic composition is equivalent to VSMOW. Using the local isotopic composition of water instead of VSMOW is particularly important in this study because the system is not pure seawater. However, even in some oceanic regions such as the Arctic, the isotopic composition of H 2 O can be substantially different from VSMOW (LeGrande and Schmidt, 2006). The definition and importance of the photosynthetic endmember for GOP calculations in different environments warrant further review (Manning et al., 2017a).
We calculate the mixed layer-integrated GOP (mmol O 2 m −2 d −1 ) and the volumetric GOP (mmol O 2 m −3 d −1 ), which is the mixed layer-integrated GOP divided by the mixed layer depth. For this time series GOP is only calculated for the mixed layer because there are no O 2 measurements below the mixed layer. The average errors in the daily GOP are +77 −49 % and +52 −31 % for the volumetric and mixed layer-integrated GOP, respectively. which has a variable Schmidt number based on temperature and salinity. The shaded grey area is the period where we do not calculate GOP due to uncertainties related to the rapidly melting ice cover and non-steady-state term.

Results and discussion
We do not present GOP estimates during the period where ice cover decreased most rapidly (16-19 April, grey shaded area in Fig. 6) due to the large uncertainties in GOP at this time. We do not have an estimate of k O 2 at intermediate or varying ice cover, nor do we have estimates of ice cover on 17 April or 19 April, because clouds obscured the satellite images on those days. Additionally, the exact timing of the change in the sign of the non-steady-state term is unclear. For example, on 18 April the calculated GOP assuming 17 is increasing at 0.67 per meg d −1 (the slope of the first best-fit line in Fig. 6b) is 58 mmol O 2 m −2 d −1 and the calculated GOP assuming 17 is decreasing at 2.99 per meg d −1 (the slope of the second best-fit line in Fig. 6b)) is 27 mmol O 2 m −2 d −1 on 18 April. The concentration of photosynthetic O 2 is highest during this time ( 17 = 48-54 per meg on 16, 18, and 19 April), despite the increasing gas transfer velocity. It is likely that GOP peaked during or soon after the end of the ice melt period and then declined to values similar to the be-ginning of the time series, but the uncertainties in GOP from 16 to 19 April are large.
Before interpreting the GOP results, we emphasize that the impact of changing mixed layer depth on the O 2 mass balance is not accounted for in our calculations. The mixed layer deepened from 8 to 23 April and then appeared to shoal from 23 to 28 April. If the 17 value below the mixed layer were higher (lower) than 17 in the mixed layer, this would cause us to overestimate (underestimate) GOP as the mixed layer deepened.
Overall, the rate of volumetric mixed layer GOP was relatively constant throughout the time series at 5.7 +3.1 −1.9 mmol O 2 m −3 d −1 . The mixed layer-integrated GOP showed larger changes with time that are related to the influence of changes in the mixed layer depth on the non-steady-state term. From 25 March through 8 April, mixed layer-integrated GOP was 4.6 +3.7 −2.6 mmol O 2 m −2 d −1 . Beginning after 8 April, the mixed layer depth began to increase and the non-steady-state calculation showed a substantial increase in mixed layerintegrated GOP to 7.2 +5.4 −3.3 mmol O 2 m −3 d −1 on 14 April, a 56 % increase. The non-steady-state GOP term is multiplied by the mixed layer depth, and therefore it increases linearly as the mixed layer deepens, causing the total mixed layerintegrated GOP to decrease. However, the non-steady-state GOP term is constant in volumetric units because the mixed layer depth cancels out of the equation. After the ice was gone, 17 began to decrease, and so did the mixed layerintegrated GOP. On the last 4 days of the time series, the mixed layer-integrated GOP was 8.7 +3.6 −6.3 mmol O 2 m −2 d −1 , 89 % higher than the average value prior to 9 April, but the mixed layer was also more than twice as deep at the end of the time series compared to the beginning.
The result that volumetric GOP was similar at the beginning and end of the time series indicates that ice-free conditions are not a pre-requisite for phytoplankton growth in this system. Currently, ecosystem dynamics within and below ice formed from fresh and brackish waters are not well understood (Salonen et al., 2009;Bertilsson et al., 2013;Hampton et al., 2015). Other investigators have shown that photosynthetic microbes can inhabit the interior, upper surface, and lower surface of ice, and tend to be most concentrated on the bottom surface (Welch et al., 1988;Cota et al., 1991;Frenette et al., 2008;Boetius et al., 2013). Traditionally, investigators have argued that ice-associated communities are most prevalent in ice formed from seawater; as salinity increases, the volume of unfrozen brines within the ice that the microbes can inhabit increases, and the bottom surface of the ice becomes more uneven, increasing bottom algal settlement efficiency (Legendre et al., 1981;Gosselin et al., 1986). However, more recently, investigators have also found algae growing within and on the bottom of freshwater ice in lakes and rivers, including locations in Canada such as the Great Lakes and the St. Lawrence River (Bondarenko et al., 2006;Frenette et al., 2008;Twiss et al., 2012;D'souza et al., 2013).
Phytoplankton can also grow in the water column beneath ice, especially under thin first-year ice (Legendre et al., 1981;Mundy et al., 2009;Arrigo et al., 2012). Bare ice transmits more light to surface waters than snow-covered ice, and melt pond-covered ice transmits 4 times as much light as bare ice (Light et al., 2008;Arrigo et al., 2012;Light et al., 2015). First-year ice in the Arctic (0.5-1.5 m thick) transmits ∼ 47 %-75 % of incident light through melt pond-covered ice and ∼ 13 %-25 % of incident light through snow-free ice (Arrigo et al., 2012;Light et al., 2015). The ice in the Bras d'Or Lake near the site of Injection 1 was ∼ 0.3 m thick on 29 March, and therefore likely similar or greater fractions of light were transmitted through the ice. Ice transmitting just 2 % of surface irradiance may support high rates of photosynthetic activity if the phytoplankton are acclimated to lower light levels (Cota, 1985). We observed melt ponds in Whycocomagh Bay during tracer injections on 31 March and frequently during visual surveys in April. The shallow mixed layer prior to ice melt (∼ 0.8 m from the beginning of the time series until 8 April) would have kept phytoplankton in the mixed layer close to the surface and therefore receiving light that penetrated through the ice.
Our O 2 mass balance techniques will record GOP by freefloating phytoplankton in the water column below the ice, as well as GOP by ice-associated phytoplankton if the O 2 they produce diffuses into the water rather than into the ice surface or atmosphere. Bottom-associated algae likely release much of their O 2 into the water column, especially for filamentous forms such as the diatoms frequently observed in Lake Erie and the Arctic (D'souza et al., 2013;Boetius et al., 2013).
We believe that the O 2 measurements at Little Narrows from the beginning of the time series through 16 April (the date when the fraction of open water in Whycocomagh Bay began to increase rapidly) are generally representative of the under-ice rates due to the rapid transport of water through Little Narrows channel. The mean current speed at Little Narrows channel at 0.5 m depth was 3.4 km d −1 . The distance between the ice edge and the sampling intake ranged from ∼ 0.4 to 1 km between the beginning of the time series and 12 April. Assuming a gas transfer velocity for O 2 of 0.5 m d −1 once the water mass encountered open water and a mixed layer depth between 0.8 and 2 m, the residence time of O 2 in the mixed layer was 1.6-4 d. Assuming a current speed of 3.4 km in open water, the water mass would have only been exposed to open water for between 0.1 and 0.3 d (20 % or less of the residence time of O 2 ). However, we recognize that the current speed near the ice edge may have been somewhat slower than the speed at the constriction at Little Narrows; therefore, the transit time in open water may have been somewhat longer than 0.1-0.3 d.

Calculation
We quantify non-steady-state NOP, incorporating the observed changes in O 2 /Ar during the time series. We estimate NOP as where k O 2 is the real-time gas transfer velocity (m d −1 ), [O 2 ] eq is the equilibrium O 2 concentration (mol m −3 ), and h is the mixed layer depth (m) . The first term on the right-hand side of Eq. (11) is the steady-state NOP term, and the second term is the non-steady-state NOP term, which is dependent on the rate of change in (O 2 /Ar) with time. To calculate the rate of change in (O 2 /Ar) with time, we resampled the data to a fixed 5 s interval (each scan of all masses took 5-6 s) and filled in gaps with a linear interpolation. Then we applied a third-order lowpass Butterworth filter with a cutoff frequency of 0.3 d −1 to generate a smooth O 2 /Ar record (Roberts and Roberts, 1978) (Fig. 7). We selected the cutoff frequency to remove the short-term variability from tides and diel changes in photosynthesis and respiration, and to minimize the number of times the inflection of the curve changed while capturing the overall trends in O 2 /Ar. Below we discuss the sensitivity of the calculated NOP to the choice of cutoff frequency. Finally, we calculated the derivative of (O 2 /Ar) with respect to time using the filtered record. We applied the same filtering method to the in situ salinity and thermocouple temperature data and used the filtered data to calculate the [O 2 ] eq and k O 2 , to prevent shortterm fluctuations in salinity and temperature from producing apparent changes in NOP. We calculate the daily NOP (from 19:30 to 19:29 local time) using the average ∂ (O 2 /Ar)/∂t based on the filtered record and the average (O 2 /Ar) (using the raw, unfiltered data).
As for the GOP calculations, our NOP calculations do not account for the impact of changing mixed layer depth, nor the impact of vertical diffusion, on the O 2 mass balance because we do not have any measurements of O 2 below the mixed layer (Munro et al., 2013;Castro Morales et al., 2013). We expect that the O 2 concentration below the mixed layer was lower than the concentration within the mixed layer, based on past measurements showing undersaturated O 2 in the subsurface waters of Whycocomagh Bay (Krauel, 1975;Strain et al., 2001). Therefore, vertical processes would likely cause us to underestimate NOP. For example, a deepening of the mixed layer between 8 and 23 April may have entrained waters with a lower O 2 concentration into the mixed layer, causing a decrease in O 2 (and O 2 /Ar). Using Eq. (11), we would interpret this decrease in O 2 as a decrease in NOP.
To calculate the uncertainty in NOP, uncertainty in k O 2 (11 % for Injection 1, 18 % for Injection 2), mixed layer depth (0.3 m, 10 %-38 %), and the non-steady-state term are propagated. Uncertainty in the non-steady-state term is based on the results using different filtering methods. A conservative 13 % error associated with the cutoff frequency choice is included in the estimates of daily NOP. Uncertainty in (O 2 /Ar) (< 0.1 %, based on the mean offset between the EIMS and the discrete samples) has a negligible impact on NOP, relative to the other sources of error. The average errors in the daily NOP are ±0.8 mmol O 2 m −2 d −1 (34 %) and ±0.3 mmol O 2 m −3 d −1 (23 %) for the mixed layerintegrated and volumetric NOP, respectively. All uncertainties are the standard deviation.
Finally, we calculate the ratio of NOP to GOP for each daily estimate (the export efficiency). This ratio is similar to an f-ratio or an e-ratio (Dugdale and Goering, 1967;Laws et al., 2000) and provides information on the fraction of GOP that is available for export out of the mixed layer (Fig. 7f). The uncertainties in the NOP/GOP ratio on each day are quite large. In a steady-state NOP and GOP calculation, the gas transfer velocity k O 2 cancels out of the equation for the NOP/GOP ratio, and therefore it is not a source of uncertainty; however, in the non-steady-state term the k O 2 does not cancel out.

NOP results and comparison of NOP and GOP
As for GOP, we do not present NOP estimates during 16-19 April due to the uncertainties associated with the rapid decrease in ice cover. There are additional uncertainties in NOP from 12 to 23 April due to the increasing mixed layer depth, which may have entrained lower-O 2 waters into the mixed layer. Historical O 2 data from Whycocomagh Bay showed that there is a strong decrease in O 2 with depth below the mixed layer; however, to our knowledge, no prior measurements during ice melt have been published (Krauel, 1975;Strain et al., 2001).
Based on the non-steady-state NOP estimates, the ecosystem was on average net autotrophic as the ice was melting, from the beginning of the time series through 15 April (mean volumetric NOP of 1.9(2.1) mmol O 2 m −3 d −1 , median 2.5 mmol O 2 m −3 d −1 ). During the (nearly) ice-free period from 20 April through the end of the time series, the community was on average net heterotrophic but with a smaller magnitude than during the start of the time series (mean volumetric NOP of −0.7(0.9) mmol O 2 m −3 d −1 , median −0.7 mmol O 2 m −3 d −1 ). When the bay was nearly full of ice cover (from the beginning of the time series until ∼ 16 April), NOP was dominated by the non-steady-state term, and this term was positive except for between 31 March and 3 April when it was negative but small in magnitude. As the ice cover retreated, the non-steady-state term decreased and became negative on 18 April. From 18 to 22 April the steady-state term was roughly equal in magnitude but opposite in sign to the non-steady-state NOP term. The NOP is more strongly negative from 23 to 25 April (volumetric NOP of −1.6(0.5) mmol O 2 m −3 d −1 ) and then on the last 2 days of the time series (O 2 /Ar) was close to 0, and so was the rate of change in (O 2 /Ar) with time (volumetric NOP of −0.7 mmol O 2 m −3 d −1 ). If the time series had continued for longer, it would have been possible to observe whether the NOP value eventually stabilized near 0 following the dynamic ice melt period or whether it continued to oscillate between periods of net autotrophy and net heterotrophy.
The ratio of NOP/GOP has large uncertainties but qualitatively follows the trends of NOP, since the mixed layerintegrated GOP was relatively constant throughout the time series, whereas the mixed layer-integrated NOP was far more variable. At the start of the time series, there are two dates where NOP/GOP > 1, which by definition is not possible. The high NOP/GOP values could be due to uncertainty in the isotopic composition of water, which enters into the GOP calculation but not the NOP calculation, and/or the non-steady-state terms for GOP and NOP. Vertical diffusion across the mixed layer is another possible cause of the estimated NOP/GOP ratios exceeding 1; because we do not have measurements below the mixed layer, we cannot compute a correction for vertical mixing. The gradients in 17 and (O 2 /Ar) are likely different with depth, leading to different magnitudes of the vertical mixing correction for NOP Figure 7. Net oxygen production at Little Narrows. Net oxygen production at Little Narrows and data used in the calculation. (a) (O 2 /Ar) and (b) temperature measurements. The blue lines are the raw data and the black line is the filtered data. NOP in (c) mixed layer-integrated and (d) volumetric units. The green diamonds and blue squares show the values of the two terms in the GOP calculation (steady state, SS, and non-steady state, NSS). (e) Mixed layer depth. (f) Export efficiency ratio (NOP/GOP). The shaded grey area is the period where we do not calculate NOP due to uncertainties related to the rapidly melting ice cover. and GOP, which could potentially lead to errors in the estimated NOP/GOP ratio (Sect. 3.6). We recommend that researchers planning future studies using O 2 -based tracers include depth profiles of O 2 /Ar and 17 in their experimental design (Castro Morales et al., 2013;Munro et al., 2013;Wurgaft et al., 2013).
The GOP results demonstrate that GOP was similar when Whycocomagh Bay was full of ice and when it was ice-free (on a volumetric basis). On the other hand, the calculated volumetric NOP was generally positive when the bay was full of ice, but began to decrease toward negative values around 12 April. The mixed layer depth was ∼ 0.8 m from the beginning of the time series until 12 April, deepened from 12 to 23 April and then increased until the end of the time series. The apparent decrease in volumetric NOP beginning around 12 April could be due to a number of factors, including an increase in respiration and recycling of organic carbon by autotrophs and/or heterotrophs, and/or vertical processes that are not taken into account in our calculations. We note that the mixed layer NOP was actually lowest and O 2 /Ar was decreasing between 23 and 27 April, a period where the mixed layer actually appeared to shoal. Therefore a deepen-ing mixed layer cannot be the only driver of the decrease in NOP following ice melt.

Comparison to other productivity estimates
To our knowledge, our data are the only published estimates of NOP and GOP in the Bras d'Or Lake system. Geen (1965) measured primary production by 14 C uptake in the Bras d'Or Lake during summer 1962-1964 and found an average uptake rate at 5 m depth, the depth of maximum photosynthesis, of ∼ 4 mmol C m −3 d −1 (50 mg C m −3 d −1 ) (Geen and Hargrave, 1966). These rates are based on 6 h daytime incubations and ignore nighttime respiratory loss of 14 C, and therefore approximate net primary production (NPP, gross primary production minus autotrophic respiration) or something between NPP and gross primary production (GPP) Marra, 2002Marra, , 2009Quay et al., 2010). Assuming a C : O 2 ratio of 1.1-1.4 (Laws, 1991), the equivalent O 2 production based on the 14 C incubations is 4.4-5.6 mmol O 2 m −3 d −1 . Thus, the 14 C-PP is between the average NOP and GOP values, as expected.
In the Bras d'Or Lake, Hargrave and Geen (1970) found, based on summertime incubations, that zooplankton grazing was sufficient to consume all of the daily primary production, indicating the estuary ecosystem metabolism may be close to balanced (NOP ∼ 0) on a daily basis. We obtained a more dynamic record of NOP, with an average volumetric rate of 1.1 ± 2.0 mmol O 2 m −3 d −1 over the entire time series.
Comparisons with in situ gas tracer-based productivity estimates in other environments are challenging because mixed layer-integrated rates are most commonly calculated, and the mixed layer in many other systems is much deeper than 0.8-3 m. In the Beaufort Gyre (Arctic Ocean) Stanley et al. (2015) estimate a steady-state GOP of 16(5) and 38(3) mmol O 2 m −2 d −1 in summer 2011 (higher ice cover) and 2012 (lower ice cover), respectively, and NOP of 3 mmol O 2 m −2 d −1 in both summers. Mixed layer depths were ∼ 10 m. Mixed layer-integrated GOP increases as the mixed layer depth increases, but trends in NOP are less clear.

Effect of physical processes on productivity estimates
There are a number of additional environmental processes that may affect the O 2 isotope and O 2 /Ar mass balance in the mixed layer but cannot be directly quantified from the time series. However, in some cases we can determine whether these processes would tend to increase or decrease our NOP and GOP estimates. The isotopic composition of freshwater within the bay may have varied during our time series, which would impact the GOP estimates. Using model results from , we calculate that the residence time of water in Whycocomagh Bay is ∼ 0.7 years for surface waters (0-10 m) with respect to freshwater input and ∼ 2 year for deep waters (10 m to bottom) with respect to exchange with the surface waters. Thus we expect the isotopic composition of the water in the mixed layer of Whycocomagh Bay to reflect some average over multiple months. For example, if a substantial portion of the meltwater entering the estuary is derived from snow rather than from ice that freezes from water within the bay, its isotopic composition will be more reflective of seasonal precipitation. If we calculate an amountweighted δ 18 O-H 2 O for meteoric water at Truro, NS, using only the months when ice was present in Whycocomagh Bay in 2013 (January-April), δ 18 O-H 2 O = −11.0(3.6) ‰ versus VSMOW, which is within the uncertainty of our annually averaged value (−9.3(3.1) ‰). In general, a lower value of δ 18 O-H 2 O will increase GOP estimates.
The freezing and melting of ice in saline waters will generate a nonlinear salinity-δ 18 O-H 2 O relationship because the δ 18 O-H 2 O value of sea ice is similar to the water from which it formed (within ∼ 2-3 ‰) (Tan and Strain, 1980;Macdonald et al., 1995Macdonald et al., , 1999, but the salinity in sea ice is substantially lower due to brine rejection (O'Neil, 1968;Weeks and Ackley, 1986). We are not able to accurately quantify the triple oxygen isotopic composition of H 2 O in the ice and in water, nor can we quantify the timing and volume of ice freezing and melt within Whycocomagh Bay (although we know when the ice cover decreased most rapidly, the ice volume was decreasing throughout our time series). The volume contribution and isotopic composition of other sources of water inputs (e.g., riverine input and melting snow) are another source of uncertainty in the calculations.
In addition to affecting the isotopic composition of H 2 O, ice melt and riverine inflow may affect the NOP and GOP calculations in other ways. If the ice melted at the upper surface (in contact with the atmosphere) and then drained through brine channels in the ice, it likely had an isotopic composition and gas ratio similar to air-equilibrated water ( 17 8 per meg and (O 2 /Ar) 0 %). Thus, water in melt ponds that was added to the water column would tend to decrease GOP, as 17 always exceeded 8 per meg in the mixed layer (Fig. 6a) and either increase or decrease NOP, since (O 2 /Ar) transitioned from negative to positive as the ice was melting. If the ice melted at the bottom (in contact with the water), its effect on NOP and GOP estimates is less clear. During sea ice formation, approximately 40 %-55 % of the O 2 and Ar originally dissolved in the water is retained in the ice matrix (i.e., within the ice itself, in gas bubbles, or in brine pockets) and the remainder is excluded, generating supersaturations of the gases in the water below the ice (Top et al., 1988;Hood et al., 1998;Loose et al., 2009). Photosynthesis and respiration both occur in sea ice (Loose et al., 2011a;Zhou et al., 2014) and will change the O 2 /Ar and 17 signatures within the ice, and it is difficult to predict what proportions of the O 2 within brine pockets in the ice remained within the brines, migrated into the water column, or migrated into the atmosphere prior to the complete melting of the ice. We measured an ice thickness of ∼ 0.3 m near the injection site on 29 March, and thus if the mixed layer depth after ice melt was 2.5 m deep, the ice could contribute ∼ 11 % of the mixed layer volume, or ∼ 5 % of the mixed layer O 2 (assuming [O 2 ] in ice is ∼ 45 % of the equilibrium [O 2 ] in water). Thus bottom-ice melt would likely be a minor influence on the oxygen mass balance. We observed bare ice and melt ponds at the surface of the ice (Figs. 2-4) and the water temperature throughout our time series was above freezing, which would stimulate bottom melt. Therefore, both surface and bottom melt likely occurred during the time series. The volume, O 2 concentration and isotopic composition of runoff and river water during our time series are also poorly constrained, and thus these water sources are another uncertainty in our NOP and GOP calculations.
It is likely that GOP occurred below the mixed layer but was not quantified by our methods because we only had measurements within the mixed layer. In open water, the Secchi depth at the ice edge on 26 March (just west of Little Narrows) was 1.9 m, yielding a euphotic zone depth of ∼ 5 m (defined as the depth where 1 % of surface photosynthetically active radiation penetrates) and the Secchi depth at Little Narrows on 7 April was ∼ 4.5 m, giving an approximate euphotic zone depth of 12 m (Idso and Gilbert, 1974).
Mixed layer depths during our time series ranged from 0.8 to 3.0 m. Even in field studies where O 2 measurements below the euphotic zone are available, it is challenging to quantify this subsurface productivity because the biological O 2 fluxes below the mixed layer are small and the physical fluxes are large and highly uncertain, as they are driven by lateral and vertical mixing rather than air-water gas exchange (Giesbrecht et al., 2012;Munro et al., 2013;Manning et al., 2017b).
We also cannot correct our results for the effect of entrainment of deeper waters into the mixed layer (Hendricks et al., 2004;Munro et al., 2013;Wurgaft et al., 2013). We were not able to sample below the mixed layer, but we hypothesize that 17 and O 2 /Ar both likely decreased below the mixed layer when the ice was melting and the mixed layer was deepening. Therefore, entrainment of these waters into the mixed layer would tend to decrease NOP and GOP estimates. Respiration increases with depth, causing [O 2 ] and (O 2 /Ar) to decrease with depth in oceanic systems (Spitzer and Jenkins, 1989;Emerson et al., 1991). In some oceanic regions, photosynthesis below the mixed layer generates excess 17 (because there is no process decreasing 17 below the mixed layer), which can then be entrained into the euphotic zone, where 17 is lower because some of the photosynthetic O 2 is ventilated to the atmosphere (Hendricks et al., 2004;Sarma et al., 2005;Juranek and Quay, 2013). In this study, the gas transfer velocity out of the mixed layer was extremely low up until 16 April, and therefore photosynthetic O 2 (i.e., 17 ) would also accumulate in the mixed layer, likely at a greater rate than below the mixed layer due to photosynthesis rates being higher closer to the surface. Once the gas transfer velocity increased to the open water value (20 April), it is not clear whether we would expect to observe an excess of 17 below the mixed layer, and therefore whether entrainment would increase or decrease GOP.
Another source of error in our interpretation of the O 2 data is that we must assume that spatial variability in O 2 has a negligible effect on our calculations. We interpret all changes in O 2 assuming that we are measuring the same water mass; this is an oversimplification because the bay is within an estuary that experiences tidal flows. However, spatial surveys of O 2 /Ar during a pilot experiment in 2011 indicated that spatial variability in O 2 within the Bras d'Or Lake is relatively low.

Conclusions
Using the dual tracer ( 3 He/SF 6 ) technique in the Bras d'Or Lake, we found that at > 90 % ice cover, the gas transfer velocity was 6 % of the open water gas transfer velocity. This result indicates that k ice is negligible.
The volumetric GOP was similar at the beginning and end of the time series, when the basin was full of ice and when it was ice-free. Volumetric NOP was more vari-able with time; Whycocomagh Bay was on average net autotrophic (NOP > 0) while the bay was ice-covered, and net heterotrophic (NOP < 0) but with a smaller magnitude after the bay was ice-free. These results indicate that an algal bloom (increasing NOP) can occur in an ice-covered estuary, similar to observations of net autotrophy under ice in the Great Lakes, Arctic, and Antarctic (Tortell and Long, 2009;Twiss et al., 2012;D'souza et al., 2013). The apparent decrease in NOP may be due to a number of factors, such as an increase in heterotrophic respiration of organic carbon and/or the entrainment of waters with a lower O 2 concentration into the mixed layer (Wurgaft et al., 2013;Castro Morales et al., 2013;Nicholson et al., 2014).
Obtaining a time series of O 2 data and obtaining simultaneous gas transfer velocity measurements were both critical for quantifying productivity. The non-steady-state term was a significant contributor to NOP and GOP throughout the time series, and time-series measurements are needed to quantify the non-steady-state O 2 flux. Additionally, because the gas transfer velocity was ∼ 16 times higher at the end of the time series than at the beginning and the values of 17 and (O 2 /Ar) changed with time, the relative importance of the steady-state term versus the non-steady-state term changed substantially during our time series.
In future studies in similar settings, we recommend authors collect measurements of O 2 /Ar and triple oxygen isotopes below the mixed layer in order to better correct for the impact of vertical processes on the O 2 mass balance. For tidally influenced regions, surveys of lateral variability in O 2 /Ar and 17 could also help to constrain the influence of horizontal advection on the O 2 mass balance (Munro et al., 2013;Howard et al., 2017). Measurements of nutrients and dissolved organic carbon could be useful for inferring the causes of changes in GOP and NOP. We also encourage measurement of the triple isotopic composition of H 2 O (Manning et al., 2017a).