Underestimation of surface pCO₂ and air-sea CO₂ fluxes due to freshwater stratification in an Arctic shelf sea, Hudson Bay

The Arctic shelves are considered one of the most sensitive locations on our planet to climate change (AMAP, 2017), experiencing fast environmental changes driven by both natural and anthropogenic factors that include sea-ice loss (e.g., Kwok et al., 2009; Notz and Stroeve, 2016), sea-surface warming (e.g., Steele et al., 2008; Shepherd, 2016), intensification of the hydrological cycle (e.g., McClelland et al., 2006; Carmack et al., 2016), and increases in primary production (e.g., Arrigo and van Dijken, 2015). The Arctic continental shelf seas cover more than 50% of the Arctic Ocean area, and they are known to play a significant role in the Arctic marine carbon budget (e.g., Bates and Mathis, 2009; Yasunaka et al., 2016). Despite the Arctic shelves making up only ca. 2% of the global ocean surface area, the estimated CO₂ sink in these waters varies from 70 to 180 TgC yr^–1 (1 TgC = 10¹² g C), accounting for 5%–15% of the global ocean CO₂ sink (Bates and Mathis, 2009; MacGilchrist et al., 2014; Yasunaka et al., 2016; Manizza et al., 2019). Yet, there are large uncertainties when budgeting the carbon sink of these regions due to the spatiotemporal variability of atmospheric CO₂ uptake, including the high seasonal variability associated with sea-ice dynamics and river runoff, and the scarcity of field measurements (e.g., Yasunaka et al., 2016; Ahmed and Else, 2019; Manizza et al., 2019). Also, doubts remain about how the ocean CO₂ sink in the Arctic will change in the future in response to climate change. Therefore, studying the Arctic shelf seas separately and in different seasons to budget the Arctic marine carbon sink accurately and get insights into how it might change due to climate change is very important.

Air-sea CO₂ fluxes are determined by the disequilibrium between CO₂ partial pressure in the atmosphere and in surface seawater (pCO_2atm and pCO₂, respectively₎. As pCO_2atm is relatively stable over short timescales, changes in pCO₂ will determine whether an area will act as a CO₂ source or sink. Hence, estimating the spatial and temporal variability of pCO₂ in a reliable way to accurately quantify and understand the processes controlling the CO₂ fluxes is important. To date, most studies that estimate regional air-sea CO₂ fluxes from observations are based largely on automated underway systems onboard research vessels that continuously measure pCO₂ at 3–7 m below the sea surface (e.g., Bakker et al., 2016). Determination of air-sea CO₂ fluxes from these continuous underway pCO₂ measurements is based on the assumption that seawater observations at depths of 3–7 m are representative of waters at the surface due to the homogeneous distribution of pCO₂ in the mixed layer. However, recent studies (e.g., Calleja et al., 2013; Miller et al., 2019) found that this assumption is not applicable everywhere and that vertical pCO₂ gradients can exist within the upper few meters of the ocean. The possible causes of such gradients can be attributed to thermal stratification (mainly caused by solar radiation), heterogeneous biological activity, or freshwater stratification (Murata et al., 2008; Calleja et al., 2013; Miller et al., 2019). In this study, we focus on the impacts of freshwater stratification on vertical CO₂ concentration gradients.

Currently, the Arctic shelves receive over 11% of global river discharge (e.g., Shiklomanov et al., 2000; Dai and Trenberth, 2002) and experience strong seasonal inputs of sea-ice meltwater (e.g., Stroeve et al., 2007; Wang and Overland, 2009). As a result, stratification introduced by freshwater (a shallow and low-salinity surface mixed layer on the top of more saline and deeper mixed layer) is a dominant characteristic of these regions (e.g., Aagaard et al., 1981; Carmack et al., 2016). Overall, freshwater stratification has been identified as a key factor controlling primary production and nutrient fluxes (e.g., Sakshaug, 2004; Ferland et al., 2011), as well as increasing the persistence of ice cover in some regions by constraining the upward diffusive flux of heat from underlying warmer layers (e.g., Polyakov et al., 2013). This stratification has the potential to cause significant errors in calculated air-sea CO₂ fluxes (Calleja et al., 2013; Miller et al., 2019).

These potential errors in air-sea CO₂ fluxes are driven by vertical pCO₂ concentration gradients within the upper few meters of the ocean that arise due to the different chemical composition of seawater, sea-ice meltwater (SIM), and river water. For example, SIM has lower concentrations of total alkalinity (TA) and dissolved inorganic carbon (DIC) than seawater (Rysgaard et al., 2011; Geilfus et al., 2015), so that mixing with SIM decreases not only salinity but also TA and DIC. Due to the typical TA:DIC ratios in sea ice, SIM tends to be undersaturated in pCO₂ (Rysgaard et al., 2011; Geilfus et al., 2015). River discharge typically has higher DIC and TA than SIM (Yamamoto-Kawai et al., 2009) but can span a range of concentrations and TA:DIC ratios, leading to variable pCO₂. Many Arctic rivers, however, are loaded with terrestrial organic carbon which, when degraded, results in pCO₂ supersaturation (e.g., Tank et al., 2012; Semiletov et al., 2016). Hence, a decrease or increase in seawater pCO₂ across a freshwater layer is usually determined by the freshwater source.

Hudson Bay is a shallow subarctic inland sea with a drainage basin covering more than a third of Canada’s land mass (Déry et al., 2011). Recognized as highly sensitive to climate variability, Hudson Bay is already experiencing significant environmental changes (e.g., ACIA, 2005; AMAP, 2017). For example, several studies have documented changes in the timing and volume of river discharge to Hudson Bay (e.g., Déry et al., 2011, 2016), effects on biota and biological activity (e.g., Ferguson et al., 2005), increase in air temperature (e.g., Joly et al., 2011), and changes in timing of sea-ice formation and melting (e.g., Hochheim and Barber, 2014). Hudson Bay receives nearly one third of Canada’s river discharge (around 760 km³) from more than 42 rivers, which is equivalent to 12% of the pan-Arctic river input (Déry et al., 2011). In addition, significant freshwater input comes from seasonal sea-ice meltwater (Granskog et al., 2011).

Freshwater stratification in Hudson Bay has already been identified during the summer season as one of the key factors governing ocean climate, nutrient fluxes, and biological productivity (e.g., Prinsenberg, 1986b; Drinkwater and Jones, 1987). Past studies observed that Hudson Bay during the summer season is characterized by a shallow summer surface mixed layer (15–60 m deep), followed by a cold winter mixed layer (60–125 m deep), and a deep water layer (largely filled with waters of Pacific origin; e.g., Prinsenberg, 1986b; Granskog et al., 2011; Burt et al., 2016). Although stratification in Hudson Bay during summer is reasonably well studied, there are no observations of stratification during the spring season.

Based on discrete surface water pCO₂ measurements (collected approximately at 2-m depth) during the fall season, Else et al. (2008a) found that Hudson Bay acts as a weak sink for atmospheric CO₂, with supersaturated pCO₂ in the coastal waters associated with river inputs and undersaturated pCO₂ in offshore waters. Else et al. (2008b) used remotely sensed sea surface temperature to deduce that the fall season CO₂ uptake is preceded by net outgassing during late summer, with a transition to atmospheric sink coinciding with the seasonal drop in sea surface temperature. Recently, Azetsu-Scott et al. (2014) and Burt et al. (2016) studied the carbonate system during the summer and fall seasons across Hudson Bay, concluding that DIC and TA variability are related strongly to the distribution of river runoff, especially in near-surface waters along the southern coast. Despite the fact that these studies provide a good baseline for understanding the carbonate system during the summer and fall seasons, there is a gap in our knowledge about pCO₂ variability and how pCO₂ varies across the top meters of the ocean as a result of the expected strong stratification during the spring season.

Here, we present recent spring–early summer (May–July 2018) seawater measurements of the carbonate system in Hudson Bay. The main objectives of this study are threefold: (1) characterize stratification in Hudson Bay during the spring season and explain its origins and characteristics, (2) understand the impact of stratification (within the upper few meters of the ocean) on the vertical and spatial variability of pCO₂, and (3) develop a method to correct for underway pCO₂ measurements in areas of strong stratification. Addressing these objectives leads to a more accurate representation of the surface distribution of pCO₂, so that we can better understand the biogeochemical processes driving surface pCO₂ variability and air-sea CO₂ fluxes across Hudson Bay and potentially in other strongly stratified regions.

The Hudson Bay system consists of Hudson Bay, James Bay, Foxe Basin, and Hudson Strait (Figure 1). Hudson Bay and Hudson Strait cover 0.84 × 10⁶ km² and 0.2 × 10⁶ km², respectively, with average depths of less than 150 m in Hudson Bay and about 300 m in Hudson Strait (Prinsenberg, 1986b; Saucier et al., 2004). Water enters Hudson Bay through two openings (Figure 1): one in the northwest called Roes Welcome Sound, where waters of predominantly Pacific origin enter from the Canadian Archipelago, and a second opening at the southeast side of Southampton Island, where waters of both Atlantic and Arctic origin enter via the north side of Hudson Strait (Wang et al., 1994; Ingram and Prinsenberg, 1998). The general circulation in Hudson Bay is mainly wind-driven and cyclonic (counterclockwise), with water outflow to the Labrador Sea along the southern side of Hudson Strait (Prinsenberg, 1986a; Ingram and Prinsenberg, 1998; Saucier et al., 2004). Overall, Hudson Bay acts as a connection basin, where waters of predominantly Pacific origin coming from the Arctic Ocean are modified biogeochemically before moving to the Labrador Sea via Hudson Strait.

Hudson Bay is covered by seasonal sea ice approximately 8–9 months of the year and transitions from complete ice coverage in the winter to complete open water in the summer. The peak of sea-ice melt in Hudson Bay occurs in June to mid-July and delivers more freshwater to the surface waters than river runoff during this time (Prinsenberg, 1988). Unlike most Arctic shelves, river runoff to Hudson Bay does not originate from one predominant river or one spatial location but is instead distributed around the perimeter of the bay. The Nelson River, draining into southeastern Hudson Bay, and La Grande Riviére, flowing to James Bay, are the largest rivers discharging to Hudson Bay, carrying about 60% of the total river input (Déry et al., 2005, 2011). Typically, river discharge to Hudson Bay is low in the winter months, followed by a rapid increase in discharge during May, driven by snowmelt and then a gradual reduction in discharge during summer months (Déry et al., 2016). However, there has been a notable change toward higher winter discharge into this region stemming from the generation of hydroelectric power (e.g., Déry et al., 2016) and a longer open-water period as a result of climate change (e.g., Hochheim and Barber, 2014).

Other sources of freshwater to Hudson Bay include net precipitation and relatively fresh surface water from Davis Strait, which carries Arctic outflow and river runoff along the northern side of Hudson Strait to Hudson Bay (Straneo and Saucier, 2008). However, these sources are very small compared to SIM and river runoff within the system (St-Laurent et al., 2011). According to the annual freshwater budget proposed by Prinsenberg (1988), river runoff adds nearly 1 m of freshwater, precipitation and evaporation together lead to an overall loss of less than half a meter of freshwater, and SIM adds nearly 1.5 m of freshwater in the summer. Overall, the estimated freshwater export from Hudson Bay system to Labrador Sea via Hudson Strait is more than 1,000 km³ yr^–1 (Straneo and Saucier, 2008), suggesting a significant contribution to the Labrador current that influences water mass transformation within the Labrador Sea (Schmidt and Send, 2007).

The data for this study were collected in Hudson Bay and Hudson Strait (Figure 1) between May 31 and July 13, 2018, onboard the research icebreaker CCGS Amundsen as part of the BaySys project. To study vertical pCO₂ concentration gradients and the influence of river runoff and SIM on those gradients, we collected water samples for salinity, DIC, TA, and the stable oxygen isotope ratio in water (δ¹⁸O) within the upper few meters of the ocean using a horizontal Niskin bottle at 47 stations. We collected bottle samples at 7-m depth (to coincide with the depth of the intake line of the underway pCO₂ system on the ship) and at the surface. Due to favorable weather conditions and limited ocean swell, we were able to confidently sample within the upper 0.5 m of the ocean using a marked rope and a horizontal Niskin bottle. Samples were collected from the ship’s foredeck or small boats away from the ship on an opportunistic basis. The number of sampled stations using a small boat was 29 (62% of the total sampled stations) during our study time. When sampling from the foredeck, the ship was usually drifting slowly onto the station to minimize water mixing. Replicate samples were collected at each station and depth by holding the Niskin bottle open for 1 min to allow for adequate flushing before closing. In some cases, we opted not to sample at a location because of seawater mixing associated with greater ship maneuvering (see Table S1 for more details).

During the cruise, we operated an automated underway pCO₂ system (General Oceanics model GO 8050; Pierrot et al., 2009) by directing water flow from a high-volume inlet located at about 7-m depth through a shower-type equilibrator at a flow rate of 2.5 L min^–1. The underway pCO₂ system includes a flow-through temperature and salinity sensor (Idronaut model Ocean Seven 315). The underway measurements were processed to remove any data acquired when the system was being cleaned or calibrated or during icebreaking operations, which often led to a blocked or restricted intake of water to the equilibrator. Additional details about the instrumentation, calibration, and data processing of the underway pCO₂ system are described in Ahmed et al. (2019). For each location where we collected a surface bottle sample, we calculated a 5-min average of pCO₂ measured by the underway system. Station locations and the horizontal distance between the surface-water sampling via Niskin bottle (from ship’s foredeck or small boats) and the closest measurement at 7 m recorded by the underway pCO₂ located on the ship are summarized in Table S1. Underway measurements were sometimes not available at the exact time of water sample collection due to data gaps as described above; when comparing the closest available underway sample, we only considered distances of <2 km from the bottle sample location.

We deployed a YSI CastAway CTD (conductivity-temperature-depth) instrument when sampling from small boats and used both the Sea-Bird (SBE 911) CTD mounted on the ship rosette and the CastAway CTD when sampling from the ship’s foredeck to characterize temperature and salinity structure in the upper meters of the ocean. These data allowed us to identify freshwater stratification spatially and vertically across the study area. Measurements of colored dissolved organic matter (CDOM) and light transmission were collected via WetLabs FLCDRTD-2344 and WetLabs C-star (CST558DR) transmissometer, respectively, which were mounted on the rosette. All sensors were factory-calibrated prior to deployment. Specifications, calibrations, and data processing for the sensors mounted on the ship’s rosette are described in Amundsen Science Data Collection (2019). The accuracy of salinity, depth, and temperature on the YSI CastAway CTD is ±0.1, ±0.01 m, and ±0.05 °C, respectively. More details about the YSI CastAway specifications are described in the user’s manual (YSI Inc., 2019).

We used the built-in derived variables in Ocean Data View (version 5.2) and the CTD profiles obtained during the expedition to calculate the Brunt–Väisälä frequency (BVF, cycles h^–1), which provides a measure of water column stratification. The frequency was calculated as $N=−(g/z)(dρ/dz)$⁠, where g is the gravitational acceleration (m s^–1), ρ is potential density (kg m^–2), z is depth (m), and N is indicative of water column stratification. The shallowest measure of N throughout the vertical CTD profile was reported. More details about the calculations of BVF are described in the Ocean Data View User’s Guide (Schlitzer, 2019). The value of BVF increases with the strength of the stratification.

Wind speeds were measured every minute onboard the ship using a conventional propeller anemometer (RM Young Co. model 15106MA). For cases when wind data from the ship were not available or for samples collected away from the ship via small boat (i.e., >2 km in distance), we used estimates of wind velocities at 10-m height provided by the National Centers for Environmental Prediction North American Regional Reanalysis (NARR; www.esrl.noaa.gov). The NARR product estimates wind every 3 h, and we used the closest estimate to our sampling time. Ahmed and Else (2019) showed that the average difference between NARR and ship wind speed measurements in the Canadian Arctic Archipelago was less than 1 m s^–1.

DIC and TA samples were collected in 125-ml borosilicate glass bottles with isoprene rubber stoppers and aluminum caps. All DIC and TA samples were preserved with a saturated HgCl₂ solution to halt biological activity and stored in the dark at 4 °C. Stored samples were analyzed within 3 months of collection at the University of Calgary, Alberta. DIC was determined using an automated infrared inorganic carbon analyzer (AIRICA, Marianda Company, Kiel, Germany) that is based on gas extraction from an acidified sample and a nondispersive infrared gas analyzer (LI-COR 7000). TA was determined by using a semiautomated open-cell potentiometric titration system (AS-ALK2, Apollo SciTech, Newark, Delaware, USA) based on the modified Gran titration method (Grasshoff et al., 1999). All analytical instruments were standardized twice a day against certified reference materials supplied by A. G. Dickson (Scripps Institution of Oceanography). The precision of our measurements for DIC and TA were ±2 and ±3 μmol kg^–1, respectively.

Salinity samples were analyzed onboard the Amundsen using a Guildline Autosal Salinometer 8400B with a precision better than ±0.002. The salinometer was calibrated with IAPSO Standard Sea Water provided by Ocean Scientific International Ltd before running the samples. As the salinity was measured by conductivity ratio based on the Practical Salinity Scale 1978 (PSS-78), it is expressed with dimensionless values. Samples for δ¹⁸O were collected in 2-ml borosilicate vials, closed tightly and sealed with Parafilm to prevent evaporation. Samples were stored in the dark at 4 °C until analysis at the University of Calgary. We analyzed the samples within 4 months using an integrated off-axis cavity absorption spectrometer (Los Gatos Research, LGR, Triple Liquid Water Isotope Analyzer, model 912–0032). The LGR analyzer was standardized every four samples with in-house standards that were calibrated against VSMOW2 (Vienna Standard Mean Ocean Water 2) and SLAP2 (Standard Light Antarctic Precipitation 2) standards. The analytical precision for δ¹⁸O was better than ±0.2%.

In this study, pCO₂ was determined in two different ways: (1) calculated from DIC and TA measured in discrete water samples and (2) measured directly by the automated underway pCO₂ system. We computed the pCO₂ from the DIC and TA measurements using the CO2SYS program (Lewis and Wallace, 1998), with the dissociation constants of carbonic acid by Lueker et al. (2000), bisulfate ion by Dickson (1990), and the boric acid by Lee et al. (2010). The pCO₂ from the underway system was calculated using the measured mole fraction of CO₂ (xCO₂) and following Dickson et al. (2007):

where P_equ is the pressure inside the equilibrator and pH₂O is the saturation water vapor pressure of air (in μatm) in the equilibrator calculated after Weiss and Price (1980). To differentiate between the methods for determining pCO₂, hereafter we refer to that computed from DIC and TA samples as “pCO_2calc” and that measured by the underway system as “pCO_2meas” To determine whether pCO_2meas and pCO_2calc are comparable, we collected six discrete samples from the intake line for comparison with coincident pCO_2meas. We found that pCO_2calc from the bottle samples tended to be lower on average (6 ± 3.8 μatm) compared to the pCO_2meas, and therefore we applied an offset correction to all pCO_2calc values. Similarly, Chen et al. (2015) found that the pCO₂ calculated from DIC and TA samples in Arctic waters were underestimated by ca. 5 μatm relative to underway pCO₂ measurements.

The fractions of sea-ice meltwater (F_SIM), meteoric water (F_MW), and seawater (F_SW) were estimated using δ¹⁸O and salinity in a three-end-member mixing model. The δ¹⁸O samples used to identify the end-member values were collected on the same cruise and analyzed using a Picarro L2120i cavity ring-down spectrometer at the University of Washington with a precision better than ±0.1%. The end-member values of salinity and δ¹⁸O for SIM were calculated based on average values of 63 ice-core samples collected from 10 different stations during the cruise. End-member values for river runoff were based on samples from 13 rivers near their point of discharge into Hudson Bay. Seawater end-member values were determined from 45 surface samples with the highest salinity from 12 different stations away from the shore with no sea ice around. The average difference in δ¹⁸O samples collected for inter-calibration between the labs at the Universities of Calgary and Washington was 0.15 ± 0.08%.

The end-member values of salinity and δ¹⁸O and their 1σ for river runoff, ice meltwater, and seawater used in this study and reported in Granskog et al. (2011) are summarized in Table 1. Our seawater end-member is similar to that determined by Granskog et al. (2011), whereas our end-member values are lower by 3.6 for salinity and 0.8% for δ¹⁸O in SIM and by 0.5 for salinity and 1.5% for δ¹⁸O for meteoric water. These differences could be attributed to the timing of data collection; we collected our samples during spring, while Granskog et al. (2011) collected their data during the fall season. To evaluate the impact of using different end-member values, we compared the water mass fractions calculated using the values reported by Granskog et al. (2011) versus ours. The different end-member values gave significant differences in F_SIM and F_MW and no significant difference in F_SW, indicating the importance of calculating the end-member properties based on seasonal field sampling. The ranges implied by standard deviations (1σ) of the end-member values we derived in this study do not introduce a significant difference to the calculated end-member fractions.

We used ice charts obtained from the Canadian Ice Service Digital Archive (CISDA) to characterize sea-ice conditions across the study area in 2018. The CISDA archive reports weekly ice conditions based on a combination of remote sensing and ship-borne observations. Several past studies have shown the reliability of the CISDA archive to monitor sea-ice cover in the Canadian Arctic (e.g., Tivy et al., 2011; Howell et al., 2013). We converted the ice charts from the native ArcInfo interchange format to shapefiles and extracted sea-ice concentration and weeks of open water since the sea-ice breakup. The weeks of open water before we sampled were defined for each location by identifying the week when ice concentration transitioned to <90%. For more details about the CISDA data processing, see Ahmed et al. (2019).

When we entered Hudson Strait on May 31, 2018, Hudson Bay and Hudson Strait were still almost entirely covered with sea ice; only northern Hudson Strait and a perennial spring-opening polynya in northwest Hudson Bay were ice-free (Figure 2). Sea ice remained present in central Hudson Bay throughout the cruise, although it was melting, particularly around the coastal margins. The spring increase in river discharge to Hudson Bay usually starts in mid-April and peaks toward the end of May (Déry et al., 2011). Therefore, the overall context for this cruise was essentially at the peak of sea-ice melt and river discharge.

To study the stratification of different water masses in Hudson Bay during the spring season, we examined two CTD transects from the ice edge to the coastal domain (Figures 3 and 4). Based on the salinity and density gradients observed along Transect 1 (June 12–15, 2018), we identified six distinct vertical layers close to the ice edge (Figure 3A, station 24), with an upper mixed layer characterized by relatively low salinity within the upper 20 m. Moving further from the sea-ice edge toward the coast, the depth of the upper mixed layer increased (Figure 3A). In addition, we observed a relatively shallow, warm surface layer between stations 26 and 28 (Figure 3A). The contour plots of salinity and BVF in Figure 3B emphasize a strongly stratified fresh layer close to the ice edge (likely due to ice meltwater) that weakened toward the shore. We did not observe any stratification caused by river runoff at the station closest to shore, which was still about 150 km away from shore.

Figure 4A shows Transect 2 (June 24–30, 2018) from the sea-ice edge to the Nelson River estuary, with data from the rosette-mounted CTD and surface data collected from a small boat at the same stations (40, 41, and 45/45A; for station 45A, the small-boat data are from the river mouth, less than 1 km from shore, but about 30 km from the ship). The CTD profiles from the small boat (especially at stations 40 and 45A; Figure 4A) showed a distinct warm and fresh layer in the upper 2 m, which was not captured by the ship’s CTD (post-processing of the ship CTD data removes the upper 2 m because of mixing caused by ship maneuvers). This fresh layer is also apparent in the contour plots of salinity, temperature, density anomaly (σ_θ), and BVF (Figure 4B) across the transect, suggesting that SIM and river runoff might have significant roles in stratifying the water in this area. However, the stratification diminishes significantly with distance from these freshwater sources (Figure 4A, station 41; Figure 4B, BVF). Although station 45A (close to the Nelson River mouth) had been sea-ice free for about 9 weeks and wind speed was relatively high (∼10 m s^–1) during our sampling, the upper few meters of the ocean were still significantly stratified, indicating that freshwater input from Nelson River maintains stratification under moderate to strong wind conditions and that stratification might be even more pronounced during calm periods.

Based on similarities in 60 CTD profiles across the study area, we can classify the water mass layers within Hudson Bay during the spring season into two domains: (1) open seawater–influenced domain and (2) sea-ice- and river-influenced domain. As an example of these two domains, we plotted the full profile of temperature, salinity, and σ_θ at two stations (stations 24 and 37) in Hudson Bay (Figure 5). The first domain, exemplified by station 24 (Figure 5a), represents areas away from the ice edge and coastal areas. This domain was characterized by a spring mixed layer (SML) with a relatively constant salinity (∼20 m deep). The second domain was mainly observed close to the ice edge and coastal areas and exemplified by station 37 (Figure 5b). This domain was characterized by a shallow stratified layer (∼2 m deep) in the upper layers of the ocean due to the influence of SIM and/or river runoff, followed by an SML.

The water column structures we observed in this study (Figure 5) are consistent with previous studies (e.g., Prinsenberg, 1986b; Granskog et al., 2011). The main difference is that we were able to capture the strong shallow freshwater layers in the upper few meters of the ocean due to strong SIM and river runoff (Figure 5b) and the absence of this layer as we moved away from ice edge and the coastal areas in Hudson Bay. Identifying the presence of these fresh and shallow stratified layers is crucial, as it is expected to cause uncertainties when using the measured CO₂ concentrations at several meters below the surface as input for calculating air-sea CO₂ fluxes.

The spatial extent of these shallow surface layers can be visualized by mapping the surface properties we measured from either small boats or the ship’s foredeck using the Niskin bottle (Figure 6). Sea surface salinity was relatively high (>31) in northern Hudson Strait and the northwest side of Hudson Bay but was relatively low (<25) along the coastal, central, and southeast side of the Bay (Figure 6a). SST values were highest (>3 °C) in the coastal areas (Figure 6b), especially around the Nelson and Churchill Estuaries. The σ_θ shows a similar spatial pattern to salinity, with the lowest values (<20 kg m^–3) observed in the coastal, central, and southeast side of the Bay (Figure 6c). These distributions are associated with variations in the fractions of meteoric water (F_MW) and sea-ice melt (F_SIM) (Figure 6e and f) at the surface, indicating the presence of stratified surface layers due to SIM and river runoff. The influences of sea-ice melt and river runoff on the surface waters are also evident in the temperature-salinity diagram (Figure S1a) and in how our bottle samples fall within the mixing triangle between the three end-members (Figure S1b).

Comparing the distributions of surface pCO_2calc to those of salinity, F_SIM and F_MW (Figure 6), we see that pCO_2calc in our study area is positively correlated with salinity (Figure 6a) but negatively correlated with F_SIM in the central and southeast side of the Bay (Figure 6e), suggesting that SIM is the dominant factor controlling the low pCO_2calc at these locations. For example, our lowest pCO_2calc observations (<150 μatm) were found in the central and southern areas of Hudson Bay, where F_SIM was highest. F_MW and pCO_2calc were both high at one coastal station on the western side of Hudson Bay between Wilson and Thlewaiaza Rivers (station 22, Figure 1), likely a result of river runoff with elevated pCO₂. Else et al. (2008a) also observed higher pCO₂ close to the river-dominated coastal regions than in offshore marine regions during the summer season in Hudson Bay.

The only other regions where we measured significant river water fractions were along the southern and southeast coasts (Figure 6f). Surprisingly, these areas had relatively low pCO_2calc, although SIM fractions were also high in these regions (Figure 6e). The variability of pCO_2calc along these coasts can be attributed most likely to variable mixing of river runoff, SIM, and seawater at each station, leading to different TA/DIC compositions. Our findings in the southern portion of the Bay agree well with Burt et al. (2016) and Azetsu-Scott et al. (2014), as they noticed a variable impact of river runoff and SIM on the carbonate system in this region. However, the magnitude of SIM impacts on pCO₂ is more prominent in our results as we were sampling during the spring season. As evidenced by our description of the two cross-bay transects (Figures 3 and 4) and previous studies of Hudson Bay (Granskog et al., 2009), river water appears to be constrained to a narrow band (100–150 km) along the coast by the prevailing circulation patterns. As a result, SIM appears to have played the dominant role in controlling the pCO₂ variability over the domain we studied (especially in the central and southeast side of the Bay) during spring.

To further examine this spatial variability of freshwater inputs in the Bay, we used all CTD profiles measured by the ship’s CTD (Figure 7). All parameters shown in Figure 7 are based on the shallowest depth recorded from the ship’s CTD, which was usually about 2 m. We observed high BVF values (>150 cycle h^–1; Figure 7a) in the central and southeast of Hudson Bay, associated with a high sea-ice melt fraction (F_SIM > 0.2; Figure 6e), indicating strong stratification at these locations due to SIM. The highest signals of river runoff (F_MW > 0.2) were found mainly along the coastal area (Figure 6f) and were associated with high CDOM (>17 mg m^–3; Figure 7b) and low light transmission (<90%; Figure 7c), as well as high BVF (>90 cycle h^–1; Figure 7a), suggesting once again that stratification due to river runoff is constrained to a narrow coastal band. Overall, we observed shallow (∼1 to 2 m) and low-salinity stratified layers at several locations across the study area, especially close to river mouths, and where the ice had recently retreated.

The existence of a shallow stratified layer in some regions of Hudson Bay suggests potentially strong vertical pCO₂ gradients. To examine this possibility, we used our bottle samples collected at the surface and a depth of 7 m. The vertical pCO₂ gradient (∇pCO₂) is defined as the differences between bottle samples collected at surface and 7-m depth, divided by the depth range (i.e., positive values indicate that pCO₂ increases with depth). Figure 8a and b shows the surface distribution of ∇pCO₂ and the vertical gradient in salinity (∇Salinity) across the study area. Not surprisingly, the ∇pCO₂ and ∇Salinity maps show a similar spatial pattern to the surface distributions of salinity and F_SIM shown in Figure 6a and e. The ∇Salinity and ∇pCO₂ were mainly positive and high in the central and southeast side of the Bay, associated with high F_SIM and low salinity values, suggesting that SIM during our sampling is generally controlling the vertical gradients in pCO₂ and salinity within the Bay. This relationship is further demonstrated by plots of ∇pCO₂ versus ∇Salinity, shaded with the fractions of SIM and river runoff (Figure 8c and d). As expected, a high vertical gradient in pCO₂ is positively associated with a high vertical gradient in salinity. The observed high ∇pCO₂ and ∇Salinity in our study area were mainly linked to high F_SIM (Figure 8c), confirming the major role of SIM on these gradients. The nonlinear relationship between ∇Salinity and ∇pCO₂ in Figure 8c and d was linked to the presence of different slopes of DIC and TA versus salinity at the surface and 7-m depth and also to the relationship between ∇DIC and ∇TA versus ∇Salinity (Figure S2).

Referring back to Figure 3, we can see a useful illustration of how the different vertical domains in Hudson Bay (i.e., Figure 5) affect pCO₂ gradients. On that transect, we observed a strong ∇pCO₂ (∼10.4 μatm m^–1) at station 24 close to the ice edge (Figure 3A) due to the dilution of upper surface water by SIM. The ∇pCO₂ decreased to 3.7 μatm m^–1 and 1.4 μatm m^–1 at stations 26 and 28, respectively, as we moved away from the ice edge. Based on the CIS weekly ice charts, station 24 was open for less than a week at the time of sampling, whereas stations 26 and 28 were open for 3 and 5 weeks, suggesting that the impact of the ice-melt stratification layer diminishes with an increasing number of weeks of since ice breakup.

Figure 9 further demonstrates the impact of surface stratification on pCO₂. Station 34, in the river-influenced coastal conduit, shows a strong stratification, with low salinity (∼27), low σ_θ (∼22.5 kg m^–3), and relatively high temperature (∼–0.6 °C) at the surface, all quite different from waters at 7 m, characterized by high salinity (∼31.5), high σ_θ (∼25.5 kg m^–3), and low temperature (∼–1.5 °C). In contrast, station 21, in the central bay, shows a well-mixed surface layer with little difference in salinity, temperature, and σ_θ over the upper 16 m. The difference between pCO₂ at the surface and 7 m at the stratified versus well-mixed stations were about 28 μatm (∇pCO₂ = 4 μatm m^–1) and 2 μatm (∇pCO₂ = 0.3 μatm m^–1), respectively. The pCO₂ difference at the stratified station is larger than the uncertainties (about 6 μatm) introduced from using different CO₂ equilibrium constants to calculate pCO₂ from DIC and TA (Lueker et al., 2000; Woosley et al., 2017). Based on the weekly CIS charts, we sampled station 34 only 1 week after ice breakup, whereas we sampled station 21 at least 5 weeks after ice breakup, suggesting that the pCO₂ differences between the surface and 7 m decline both with distance from the river mouths and with time after ice melt.

To determine the persistence of these stratified layers, we compared ∇pCO₂ and ∇Salinity against wind speed and weeks since ice breakup (Figure 10). The highest impacts on ∇pCO₂ and ∇Salinity were observed within the first week after ice breakup (Figure 10b) and under light wind conditions (Figure 10a). Observations of particularly high ∇pCO₂ and ∇Salinity were restricted to the central and southeast portions of the Bay (Figure 8a and b), where surface salinity was typically less than 27 (Figure 6a). Both ∇pCO₂ and ∇Salinity gradients diminished with increasing wind speed (>6 m s^–1) and the number of weeks since sea-ice breakup (>5 weeks). Our observations are consistent with past studies identifying a strong bias in vertical dissolved gas concentration at low wind speeds (Fischer et al., 2019) and the transient nature of ice melt lenses (Else et al., 2012; Geilfus et al., 2015). These findings suggest that the effects of stratification associated with sea-ice melt are short-lived and that wind mixing is capable of breaking down the stratification and dramatically lowering ∇pCO₂.

This short-lived impact of stratification on the vertical gradient of pCO₂ associated with sea-ice melt might be reinforced by biological CO₂ drawdown if a phytoplankton bloom also occurs at this time_. However, chlorophyll a (Chl a) concentrations (average about 1.1 µg L^–1) were mostly low throughout the bay during our study time, with high values (>3 µg L^–1) observed only at a few locations. These results suggest that biological activity more likely has a minor role in lowering surface water pCO₂ compared to sea-ice melt during our study time. That said, the role of biological activity in lowering pCO₂ might be underestimated, as we possibly missed the peak of the spring phytoplankton bloom that occurred before or soon after the ice breakup. Although we were sampling within the upper 0.5 m of the ocean, under very calm conditions, the pCO_2calc could be different from the “true” surface pCO₂ in the sea-surface microlayer directly in contact with the atmosphere as a result of CO₂ slow diffusion, near-surface temperature gradients, and/or reactions occurring in the microlayer (Ward et al., 2004; Garbe et al., 2014). As there is no systematic way to measure pCO₂ in the microlayer in the field (Cunliffe and Wurl, 2014), sampling pCO₂ close to the surface was beyond the capabilities and scope of this study. However, the difference in pCO₂ between the microlayer and the bulk water at 0.5-m depth would be substantial only under calm conditions, with little or no wind, limited near-surface turbulence, and thereby low gas transfer velocity (Wanninkhof and Knox, 1996; Wurl et al., 2016). We rarely encountered such calm conditions in Hudson Bay, concluding for the most part that our measurements at or less than 0.5-m depth are likely representative of the waters exchanging with the atmosphere.

As we discussed in the previous sections, using an underway pCO₂ system in a highly stratified region such as Hudson Bay (especially during the spring season) could potentially lead to substantial errors when calculating the atmospheric CO₂ sink. Therefore, developing a method to correct for underway pCO₂ measurements in our study area is important. Although Miller et al. (2019) found a substantial error in CO₂ flux calculations using the shipboard surface pCO₂ measurements in the coastal Canadian Arctic, they were not able to correct for these errors due to variations in the freshwater source from one station to another and vertical variations in primary productivity. Also, their stations were scattered through three different oceanographic settings (Hudson Bay, Canadian Arctic Archipelago, and Baffin Bay). According to our findings in Figure 10b, the large variances of ∇pCO₂ observed during the first week after ice breakup (i.e., very close to the ice edge) are associated with a variable SIM pulse that makes correcting the underway pCO₂ measurements at these locations challenging. However, the relatively moderate variations in ∇pCO₂ between Weeks 1 and 5 after ice breakup could possibly be corrected. Furthermore, as the magnitude of ∇pCO₂ approaches zero beyond Week 5, we argue that no correction is required, as the pCO₂ difference between surface and 7 m is usually lower than the uncertainties introduced from calculating pCO₂ from the bottle samples.

Considering all the stations where we sampled between 1 and 5 weeks after ice breakup (Figure 10b), we found a significant relationship (R = 0.80) between the pCO_2calc at the surface and the pCO_2meas at 7 m (Figure 11), providing a possible correction of the underway pCO₂ measurements to the actual surface pCO₂ conditions. The average values of ∇pCO₂ and ∇Salinity for the stations used in Figure 11 are 3.8 μatm m^–1 and 0.2 m^–1, respectively. Furthermore, the slope and offset of the relationship in Figure 11 imply that the size of the underway pCO₂ correction decreases (i.e., gets closer to a 1:1 relationship) when underway pCO₂ increases. The rationale for using a subsurface underway pCO₂ measurement to assess the correction factor necessary to infer surface pCO₂ is that the processes producing strong vertical gradients also appear to produce lower pCO₂ at 7 m. This point is clearly made by the colormap of ∇Salinity in Figure 11, where underway pCO₂ measurements below approximately 400 μatm tend to be associated with high ∇Salinity, while underway pCO₂ measurements above 400 μatm are associated with low or near-zero ∇Salinity gradients. Also, Figure 8 shows that ∇pCO₂ and ∇Salinity are strongly correlated. The underway pCO₂ measurements may be lower under strong stratification conditions due to the entrainment of low-pCO₂ ice-melt water to 7 m or due to the coincident timing of under-ice (or bottom ice) blooms with ice breakup and melt.

We tested the impact of applying this underway pCO₂ correction (corrected obs. = 1.09 × uncorrected obs. – 49.44) by calculating the average air-sea CO₂ fluxes for the stations used in Figure 11 following the bulk flux equation:

where F_CO2 is the CO₂ flux and α is the CO₂ solubility in seawater calculated according to Weiss (1974), and k is the gas transfer velocity calculated according to Wanninkhof (2014). The average calculated air-sea CO₂ fluxes using the measured and corrected pCO₂ are ∼ –3 ± 2.4 and –6 ± 3.2 mmol CO₂ m^–2 day^–1, respectively. Therefore, using underway pCO₂ measurements without correction might cause an error of ∼50% less CO₂ uptake when calculating the oceanic CO₂ sink.

Such an error for 5 weeks extrapolated over our study area would result in a difference in the annual CO₂ flux by approximately 0.16 TgC yr^–1. An error of 3 mmol CO₂ m^–2 day^–1 for a similar 5-week period over the entire Arctic shelves (excluding the Arctic basin) would result in about a 9.3 TgC yr^–1 annual difference, which would constitute to ∼5% of the annual CO₂ uptake by the Arctic Ocean estimated by Yasunaka et al. (2016), indicating that the influence of stratification on air-sea CO₂ fluxes estimated by shipboard underway pCO₂ systems should not be ignored. This error, however, may be an overestimate as it assumes that all Arctic shelves follow a stratification pattern similar to what we observed in our study region. Hudson Bay is a unique Arctic/sub-Arctic shelf that receives considerable freshwater inputs from rivers and sea-ice melt during the spring season (Prinsenberg, 1988; Déry et al., 2016), which might not be the case for other shelves. In fact, Hudson Bay receives the second-largest volume of river discharge among Arctic shelves (Déry et al., 2016) and is only surpassed by the Kara Sea, which receives most of the freshwater supply from the Ob and Yenisei rivers. Nevertheless, the large vertical gradients of salinity and pCO₂ previously observed in contrasted Arctic environments, including the Beaufort Sea (Murata et al., 2008) as well as Hudson Bay, the Canadian Arctic Archipelago and Baffin Bay (Miller et al., 2019), are comparable to our observations in the central and eastern side of the bay and suggest that freshwater stratification is strong and ubiquitous over Arctic shelves. As a result, we believe that the magnitude of the error should not be ignored and might be exacerbated in the future in response to increases in the length of the melt season (e.g., Markus et al., 2009) and the freshwater stratification associated with the predicted rise in river discharge across the Arctic domain (e.g., Nummelin et al., 2015). Furthermore, we assume that freshwater stratification will cause vertical concentration gradients of other dissolved gases, particularly for those with similar solubility in seawater as CO₂.

Although we cannot apply a universal correction for underway pCO₂ measurements in our study region due to the variable and seasonal impact of SIM (especially close to the ice edge) and river waters (particularly close to shore) on surface pCO₂, our results indicate that at least regional and seasonal corrections are possible. Based on our previous discussion that the correction we have derived for the underway measurements is relevant to the period between 1 and 5 weeks after the sea-ice breakup, we visualized the spatial area over which our correction would be appropriate during our study (Figure 12). We found that the correction would be mostly appropriate in the north of Hudson Strait and along the western and southeastern sides of Hudson Bay (Figure 12). In order to determine air-sea CO₂ exchanges in ice-covered oceans accurately, we recommend not using underway observations in locations where the ice is currently melting (especially if the ship’s intake line is below 2-m depth) but rather using only discrete surface samples as close to the surface as possible.

We studied surface freshwater stratification and its impacts on the vertical pCO₂ concentration gradient in Hudson Bay during the spring season. We found that sea-ice meltwater and river runoff create a distinct shallow stratified layer in the upper few meters of the ocean, although sea-ice melt was the dominant process leading to stratification across the region during our study. This stratified layer leads to strong near-surface CO₂ gradients (especially between 1 and 5 weeks since ice breakup) and introduces substantial uncertainty when calculating air-sea CO₂ fluxes from a shipboard underway pCO₂ system drawing water from several meters below the surface. In this study, we estimated a 50% lower atmospheric CO₂ uptake when using underway pCO₂ measurements at 7-m depth instead of the surface pCO₂. We derived a correction to align the underway pCO₂ measurements to the surface conditions that is applicable between 1 and 5 weeks after ice breakup. This impact of freshwater stratification on air-sea CO₂ fluxes may not be limited to polar waters, as the coastal areas at all latitudes host many freshwater sources. Additionally, we assume the air-sea gas exchange of other gases might be affected in different ways by near-surface freshwater stratification. As automated underway pCO₂ systems are widely used for monitoring surface pCO₂, we argue that collecting coincident discrete water samples as close to the surface as possible, especially in stratified environments such as Arctic and river-influenced shelves, is important to correct the shipboard underway pCO₂ system and accurately estimate air-sea CO₂ fluxes.

The data used in this study are available here:

Ahmed, M, Else, B. 2019. Discrete measurements of dissolved inorganic carbon (DIC), total alkalinity (TA), temperature and salinity during the CCGS Amundsen cruises in Hudson Bay and Hudson Strait from 2018–05-31 to 2018-06-24 (NCEI Accession 0206288). NOAA National Centers for Environmental Information, Dataset, https://doi.org/10.25921/hmfe-nj95.

The supplemental files for this article can be found as follows:

• Figure S1. Water mass characteristics at each station based on our bottle samples. (a) Surface-water potential temperature versus salinity, color-shaded with pCO₂ values, with isopycnals of density anomaly (σ_θ) displayed in dashed lines (blue rectangle indicates samples mainly from Hudson Strait), and (b) oxygen water isotope (δ¹⁸O) versus salinity, shaded with pCO₂ values (the dashed lines represent the mixing lines between the end-members). SIM, SW, and MW refer to sea-ice melt, seawater, and meteoric water end-members, respectively.

• Figure S2. Vertical gradients of DIC and TA in Hudson Bay. Distributions of DIC (a) and TA (b) as a function of salinity at the surface (red) and at 7-m depth (blue), with the distributions of vertical gradients of DIC (c) and TA (d) versus the vertical gradients of salinity during the study time. The solid black lines in (a, b) show linear regression relationships.

• Table S1. Sampling stations with surface and 7-m water samples collected across the study area.

We wish to thank the captain and the crew of the CCGS Amundsen for their hard work and help during the BaySys cruise in 2018, and chief scientist, David Barber, for his leadership and support. We also acknowledge Rachel Mandryk and Yekaterina (Kate) Yezhova for help during sampling; Brian Butterworth, Tonya Burgers, and Sebastian Luque for logistics and help during ship mobilization; Patrick Duke, Stephen Gonski, and Araleigh Cranch for their help during lab analyses. The ship TSG and CTD data presented herein were collected by systems onboard the Canadian research icebreaker CCGS Amundsen and made available by the Amundsen Science program. Some figures in this article were made using Ocean Data View (https://odv.awi.de/).

The work was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) Collaborative Research and Development (CRD) program and Manitoba Hydro.

The authors have no competing interests to declare. Lisa Miller is an editor at Elementa. She was not involved in the review of the article.

• Contributed to conception and design: MA.

• Contributed to the acquisition of data: MA, DC, TP.

• Contributed to analysis and making the figures: MA.

• Contributed to interpretation of results: MA, BE, DC, LM, TP.

• Drafted and/or revised the article: MA, BE, DC, LM, TP.

• Approved the submitted version for publication: MA, BE, DC, LM, TP.

How to cite this article: Ahmed, MMM, Else, BGT, Capelle, D, Miller, LA, Papakyriakou, T. 2020. Underestimation of surface pCO₂ and air-sea CO₂ fluxes due to freshwater stratification in an Arctic shelf sea, Hudson Bay. Elem Sci Anth, 8: xx. DOI: https://doi.org/10.1525/elementa.084

Domain Editor-in-Chief: Jody W. Deming, University of Washington, WA, USA

Associate Editor: Jean-Eric Tremblay, Université Laval, Canada

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: BaySys