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

The objective of this study is to quantify the impact of freshwater stratification on the vertical gradients of partial pressure of CO 2 (pCO 2 ) and estimates of air-sea CO 2 exchange in Hudson Bay during peak sea-ice melt and river runoff. During the spring of 2 0 18, we sampled water in Hudson Bay and Hudson Strait for dissolved inorganic carbon, total alkalinity, salinity, the oxygen stable isotope ratio in the water ( d 18 O), and other ancillary data. The coastal domain and regions close to the ice edge had significant vertical concentration gradients of pCO 2 across the top meters of the ocean due to the presence of a stratified fresh layer at the surface. The pCO 2 and salinity in the central (where sea-ice melt was significant) and the southeast (where river runoff and sea-ice melt were significant) side of the bay generally increased with depth, with average gradients of 4.5 m atm m –1 and 0 .5 m –1 , respectively. Ignoring these gradients causes a significant error in calculating air-sea CO 2 fluxes, especially when using shipboard underway systems that measure pCO 2 at several meters below the sea surface. We found that the oceanic CO 2 sink in Hudson Bay is underestimated by approximately 5 0 % if underway pCO 2 system measurements are used without correction. However, we observed that these gradients do not persist for more than 5 weeks following ice melt. We have derived a linear correction for underway pCO 2 measurements to account for freshwater stratification during periods of 1–5 weeks after ice breakup. Given the lack of measurements in stratified Arctic waters, our results provide a road map to better estimates of the important role of these regions in global carbon cycles. Some


Introduction
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), seasurface 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 2 sink in these waters varies from 70 to 180 TgC yr -1 (1 TgC ¼ 10 12 g C), accounting for 5%-15% of the global ocean CO 2 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 2 uptake, including the high seasonal variability associated with seaice dynamics and river runoff, and the scarcity of field measurements (e.g., Yasunaka et al., 2016;Manizza et al., 2019). Also, doubts remain about how the ocean CO 2 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 2 fluxes are determined by the disequilibrium between CO 2 partial pressure in the atmosphere and in surface seawater (pCO 2atm and pCO 2 , respectively ) . As pCO 2atm is relatively stable over short timescales, changes in pCO 2 will determine whether an area will act as a CO 2 source or sink. Hence, estimating the spatial and temporal variability of pCO 2 in a reliable way to accurately quantify and understand the processes controlling the CO 2 fluxes is important. To date, most studies that estimate regional air-sea CO 2 fluxes from observations are based largely on automated underway systems onboard research vessels that continuously measure pCO 2 at 3-7 m below the sea surface (e.g., Bakker et al., 2016). Determination of air-sea CO 2 fluxes from these continuous underway pCO 2 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 2 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 2 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 2 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 2 fluxes (Calleja et al., 2013;Miller et al., 2019).
These potential errors in air-sea CO 2 fluxes are driven by vertical pCO 2 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 2 (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 2 . Many Arctic rivers, however, are loaded with terrestrial organic carbon which, when degraded, results in pCO 2 supersaturation (e.g., Tank et al., 2012;Semiletov et al., 2016). Hence, a decrease or increase in seawater pCO 2 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., 2011Déry et al., , 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 3 ) 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 2 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 2 , with supersaturated pCO 2 in the coastal waters associated with river inputs and undersaturated pCO 2 in offshore waters. Else et al. (2008b) used remotely sensed sea surface temperature to deduce that the fall season CO 2 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 2 variability and how pCO 2 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 2 , and (3) develop a method to correct for underway pCO 2 measurements in areas of strong stratification. Addressing these objectives leads to a more accurate representation of the surface distribution of pCO 2 , so that we can better understand the biogeochemical processes driving surface pCO 2 variability and air-sea CO 2 fluxes across Hudson Bay and potentially in other strongly stratified regions.

Oceanographic setting and freshwater sources
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 6 km 2 and 0.2 Â 10 6 km 2 , 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(Déry et al., , 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 3 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).

Field sampling
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 2 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 (d 18 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 2 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 2 system (General Oceanics model GO 8050; Pierrot et al., 2009) by directing water flow from a highvolume 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 2 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 2 system are described in . For each location where we collected a surface bottle sample, we calculated a 5-min average of pCO 2 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 2 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 , where g is the gravitational acceleration (m s -1 ), r 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.  showed that the Art. XX, page 4 of 21 Ahmed et al: Freshwater stratification and p CO 2 gradients in Hudson Bay average difference between NARR and ship wind speed measurements in the Canadian Arctic Archipelago was less than 1 m s -1 .

Analytical methods
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 2 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 mmol 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 d 18 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 d 18 O was better than +0.2%.

pCO 2 calculations
In this study, pCO 2 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 2 system. We computed the pCO 2 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 2 from the underway system was calculated using the measured mole fraction of CO 2 (xCO 2 ) and following Dickson et al. (2007): where P equ is the pressure inside the equilibrator and pH 2 O is the saturation water vapor pressure of air (in matm) in the equilibrator calculated after Weiss and Price (1980). To differentiate between the methods for determining pCO 2 , 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 matm) 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 2 calculated from DIC and TA samples in Arctic waters were underestimated by ca. 5 matm relative to underway pCO 2 measurements.

Freshwater and seawater fractions
The fractions of sea-ice meltwater (F SIM ), meteoric water (F MW ), and seawater (F SW ) were estimated using d 18 O and salinity in a three-end-member mixing model.  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 d 18 O in SIM and by 0.5 for salinity and 1.5% for d 18 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 endmember 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 (1s) of the end-member values we derived in this study do not introduce a significant difference to the calculated end-member fractions.

Sea-ice data
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 .

Stratification in Hudson Bay
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 (s y ), 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 s y 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 2 concentrations at several meters below the surface as input for calculating air-sea CO 2 fluxes.

Surface water characteristics
The spatial extent of these shallow surface layers can be visualized by mapping the surface properties we measured Art. XX, page 6 of 21 Ahmed et al: Freshwater stratification and p CO 2 gradients in Hudson Bay 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 s y 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 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 2 is more prominent in our results as we were sampling during the spring season. As evidenced by  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 2 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 seaice 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.

Impacts of stratification on vertical pCO 2 gradients
The existence of a shallow stratified layer in some regions of Hudson Bay suggests potentially strong vertical pCO 2 gradients. To examine this possibility, we used our bottle samples collected at the surface and a depth of 7 m. The vertical pCO 2 gradient (rpCO 2 ) 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 2 increases with depth). Figure 8a and b shows the surface distribution of rpCO 2 and the vertical gradient in salinity (rSalinity) across the study area. Not surprisingly, the rpCO 2 and rSalinity maps show a similar spatial pattern to the surface distributions of salinity and F SIM shown in Figure 6a and e. The rSalinity and rpCO 2 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 2 and salinity within the Bay. This relationship is further demonstrated by plots of rpCO 2 versus rSalinity, shaded with the fractions of SIM and river runoff (Figure 8c and d). As expected, a high vertical gradient in pCO 2 is positively associated with a high vertical gradient in salinity. The observed high rpCO 2 and rSalinity 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 rSalinity and rpCO 2 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 rDIC and rTA versus rSalinity ( 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 2 gradients. On that transect, we observed a strong rpCO 2 (*10.4 matm m -1 ) at station 24 close to the ice edge ( Figure 3A) due to the dilution of upper surface water by SIM. The rpCO 2 decreased to 3.7 matm m -1 and 1.4 matm 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 2 . Station 34, in the river-influenced coastal conduit, shows a strong stratification, with low salinity (*27), low s y (*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 s y (*25.5 kg m -3 ), and low Art. XX, page 10 of 21 Ahmed et al: Freshwater stratification and p CO 2 gradients in Hudson Bay 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 s y over the upper 16 m. The difference between pCO 2 at the surface and 7 m at the stratified versus well-mixed stations were about 28 matm (rpCO 2 ¼ 4 matm m -1 ) and 2 matm (rpCO 2 ¼ 0.3 matm m -1 ), respectively. The pCO 2 difference at the stratified station is larger than the uncertainties (about 6 matm) introduced from using different CO 2 equilibrium constants to calculate pCO 2 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 2 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 rpCO 2 and rSalinity against wind speed and weeks since ice breakup ( Figure 10). The highest impacts on rpCO 2 and rSalinity were observed within the first week after ice breakup (Figure 10b) and under light wind conditions (Figure 10a). Observations of particularly high rpCO 2 and rSalinity 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 rpCO 2 and rSalinity 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 rpCO 2 . This short-lived impact of stratification on the vertical gradient of pCO 2 associated with sea-ice melt might be reinforced by biological CO 2 drawdown if a phytoplankton bloom also occurs at this time . However, chlorophyll a (Chl a) concentrations (average about 1.1 mg L -1 ) were mostly low throughout the bay during our study time, with high values (>3 mg 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 2 compared to sea-ice melt during our study time. That said, the role of biological activity in lowering pCO 2 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 2 in the sea-surface microlayer directly in contact with the atmosphere as a result of CO 2 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 2 in the microlayer in the field (Cunliffe and Wurl, 2014), sampling pCO 2 close to the surface was beyond the capabilities and scope of this study. However, the difference in pCO 2 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.

Surface versus underway pCO 2
As we discussed in the previous sections, using an underway pCO 2 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 2 sink. Therefore, developing a method to correct for underway pCO 2 measurements in our study area is important. Although Miller et al. (2019) found a substantial error in CO 2 flux calculations using the shipboard surface pCO 2 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 rpCO 2 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 2 measurements at these locations challenging. However, the relatively moderate variations in rpCO 2 between Weeks 1 and 5 after ice breakup could possibly be corrected. Furthermore, as the magnitude of rpCO 2 approaches zero beyond Week 5, we argue that no correction is required, as the pCO 2 difference between surface and 7 m is usually lower than the 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 2 measurements to the actual surface pCO 2 conditions. The average values of rpCO 2 and rSalinity for the stations used in Figure 11 are 3.8 matm 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 2 correction decreases (i.e., gets closer to a 1:1 relationship) when underway pCO 2 increases. The rationale for using a subsurface underway pCO 2 measurement to assess the correction factor necessary to infer surface pCO 2 is that the processes producing strong vertical gradients also appear to produce lower pCO 2 at 7 m. This point is clearly made by the colormap of rSalinity in Figure 11, where underway pCO 2 measurements below approximately 400 matm tend to be associated with high rSalinity, while underway pCO 2 measurements above 400 matm are associated with low or nearzero rSalinity gradients. Also, Figure 8 shows that rpCO 2 and rSalinity are strongly correlated. The underway pCO 2 measurements may be lower under strong stratification conditions due to the entrainment of low-pCO 2 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 2 correction (corrected obs. ¼ 1.09 Â uncorrected obs. -49.44) by calculating the average air-sea CO 2 fluxes for the stations used in Figure 11 following the bulk flux equation: where F CO2 is the CO 2 flux and a is the CO 2 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 2 fluxes using the measured and corrected pCO 2 are * -3 + 2.4 and -6 + 3.2 mmol CO 2 m -2 day -1 , respectively. Therefore, using underway pCO 2 measurements without correction might cause an error of *50% less CO 2 uptake when calculating the oceanic CO 2 sink. Such an error for 5 weeks extrapolated over our study area would result in a difference in the annual CO 2 flux by approximately 0.16 TgC yr -1 . An error of 3 mmol CO 2 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 2 uptake by the Arctic Ocean estimated by Yasunaka et al. (2016), indicating that the influence of stratification on air-sea CO 2 fluxes estimated by shipboard underway pCO 2 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 secondlargest 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 2 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 , 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 2 .
Although we cannot apply a universal correction for underway pCO 2 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 2 , 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 Figure 11. Underway pCO 2 correction. Scatterplot of surface pCO 2calc , derived from bottle samples, against pCO 2meas at 7-m depth, measured by an underway pCO 2 system located on the ship. We plotted only stations where the surface samples were collected within 2-km distance from the ship and where 1-5 weeks had elapsed since ice breakup at the time of sampling. The regression relationship is represented by the black solid line, and the dashed black line indicates the 1:1 relationship. The prediction intervals for the fitted regression line at 95% are shown in gray. The regression equation, root mean square error, standard error of the regression slope, and the correlation coefficient (R) are also presented.DOI: https://doi.org/ 10.1525/elementa.084.f11 Figure 12. The spatial area where underway correction would be appropriate at the time of our study. Surface distributions of the number of weeks since ice breakup with the vertical gradient of pCO 2 (rpCO 2 ) are shown as black contour lines across the study area. The red contour lines represent 1 and 5 weeks since sea-ice breakup, with the spatial area between both lines indicating where correction is necessary. The white area represents sea-ice cover (>9/10) as of July 9, 2018, based on weekly ice charts provided by the Canadian Ice Service. DOI: https://doi.org/10. 1525/elementa.084.f12 Ahmed et al: Freshwater stratification and p CO 2 gradients in Hudson Bay Art. XX, page 15 of 21 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 2 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.

Conclusions
We studied surface freshwater stratification and its impacts on the vertical pCO 2 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 2 gradients (especially between 1 and 5 weeks since ice breakup) and introduces substantial uncertainty when calculating air-sea CO 2 fluxes from a shipboard underway pCO 2 system drawing water from several meters below the surface. In this study, we estimated a 50% lower atmospheric CO 2 uptake when using underway pCO 2 measurements at 7-m depth instead of the surface pCO 2 . We derived a correction to align the underway pCO 2 measurements to the surface conditions that is applicable between 1 and 5 weeks after ice breakup. This impact of freshwater stratification on airsea CO 2 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 2 systems are widely used for monitoring surface pCO 2 , 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 2 system and accurately estimate air-sea CO 2 fluxes.

Supplemental files
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 2 values, with isopycnals of density anomaly (s y ) displayed in dashed lines (blue rectangle indicates samples mainly from Hudson Strait), and (b) oxygen water isotope (d 18 O) versus salinity, shaded with pCO 2 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.