Widespread surface water pCO2 undersaturation during ice-melt season in an Arctic continental shelf sea (Hudson Bay, Canada)

Estimating sea–air CO2 fluxes in coastal seas remains a source of uncertainty in global carbon budgets because processes like primary production, upwelling, water mixing, and freshwater inputs produce high spatial and temporal variability of CO2 partial pressure (pCO2). As a result, improving our pCO2 baseline observations in these regions is important, especially in sub-Arctic and Arctic seas that are experiencing strong impacts of climate change. Here, we show the patterns and main processes controlling seawater pCO2 and sea–air CO2 fluxes in Hudson Bay during the 2018 spring and early summer seasons. We observed spatially limited pCO2 supersaturation (relative to the atmosphere) near river mouths and beneath sea ice and widespread undersaturated pCO2 in offshore and ice-melt-influenced waters. pCO2 was highly correlated with salinity and temperature, with a limited but statistically significant relationship with chlorophyll a and fluorescent dissolved organic matter. Hudson Bay on average was undersaturated with respect to atmospheric CO2, which we attribute mainly to the dominance of sea-ice meltwater. We calculated an average net CO2 flux of about –5mmol CO2 m day (–3.3 Tg C) during the spring and early summer seasons (92 days). Combining this result with extrapolated estimates for late summer and fall seasons, we estimate the annual CO2 flux of Hudson Bay during the open water season (184 days) to be –7.2 Tg C. Our findings indicate that the bay on average is a weaker CO2 sink than most other Arctic seas, emphasizing the importance of properly accounting for seasonal variability in the Arctic coastal shelves to obtain reliable sea–air CO2 exchange budgets.


Introduction
Global oceans slow the effects of climate change by absorbing anthropogenic carbon dioxide (CO 2 ) emitted by the burning of fossil fuels, land-use change, and cement production. However, the invasion of this anthropogenic CO 2 into the ocean has caused a decline in surface ocean pH leading to ocean acidification (Orr et al., 2005;Gattuso and Hansson, 2011). Friedlingstein et al. (2019) quantified the major components of the global carbon budget (atmosphere, land, and ocean reservoirs) over the previous decade (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) and estimated that the global ocean absorbed about 23% of the total anthropogenic CO 2 emissions. However, this global budget neglected the carbon cycle in coastal ocean waters largely because the spatial and temporal heterogeneity of biogeochemical processes in these waters make it difficult to accurately scale up results from observational studies. Despite great efforts over the last decade to assess the coastal contribution to the global carbon budget (e.g., Chen and Borges, 2009;Bauer et al., 2013;Laruelle et al., 2014;Bourgeois et al., 2016;Roobaert et al., 2019), the CO 2 budgeting of these waters is still subject to uncertainty and lacks consensus. The coastal ocean has been impacted significantly by human activities (hydroelectric damming, eutrophication, overfishing, etc.) and is now considered one of the most sensitive parts of the marine environment to climate change (e.g., Arctic Climate Impact Assessment, 2005; Liu et al., 2010;Liao et al., 2015;Breitburg et al., 2018). Furthermore, the coastal shelves play a disproportionately significant role in the global carbon cycle in comparison with its surface area (about 7%-10% of the oceanic surface area) by linking the oceanic, terrestrial, and atmospheric carbon reservoirs (Wollast, 1998;Liu et al., 2010). Hence, studying the carbon cycle in these waters is essential to understand future shifts in response to climate change.
The coastal ocean consists of diverse but closely connected ecosystems that include estuaries, rivers, tidal wetlands, and the continental shelf. The continental (or coastal) shelves are shallow and have been the focus of many studies due to their key role in the global biogeochemical cycling of carbon and nutrients (e.g., Chen et al., 2013;Bourgeois et al., 2016) as well as their importance to local livelihoods. For instance, the ocean's biological pump operates more efficiently in coastal shelves relative to the open ocean as a result of high terrestrial inputs and efficient use of nutrients (Gattuso et al., 1998;Chavez et al., 1999). The high biological production in these waters usually leads to atmospheric CO 2 drawdown and subsequent CO 2 export to the subsurface waters via the biological pump (Ducklow et al., 2001). The eventual outflow of CO 2 -enriched deep water from coastal shelf regions to the subsurface waters of the adjacent deep oceans constitutes "the continental shelf pump," which is thought to contribute significantly to the global ocean's atmospheric CO 2 uptake (Tsunogai et al., 1999). Furthermore, the production, export, and degradation of organic matter in these waters are several times higher than in the open ocean (Wollast, 1998;Borges et al., 2005). The coastal shelves contribute 10%-30% of global marine primary production, 80% of organic carbon, and 30%-50% of inorganic carbon burial in sediments, and about 50% of the organic carbon supplied to the deep open ocean (Wollast, 1998;Liu et al., 2010). As a result, the sea-air CO 2 exchanges are expected to be more intense in these waters than in the open ocean.
Coastal shelves show a latitudinal pattern in CO 2 fluxes, with the mid-to-high-latitude (30 -90 ) shelves acting as oceanic sinks of CO 2 , whereas the low latitude shelves (equator to 30 ) act as sources of CO 2 to the atmosphere or are nearly neutral (Borges et al., 2005;Cai et al., 2006). This latitudinal pattern can be attributed to a complex interaction among temperature (which alters the solubility of CO 2 ), the amount of biological activity, the magnitude of upwelling and mixing from deeper waters, and the export of terrestrial organic and inorganic carbon from rivers (Takahashi et al., 2009;Dai et al., 2012). Adding further complexity, the processes controlling CO 2 dynamics vary with proximal distance from shore (Wollast, 1998;Bauer et al., 2013). Inner shelf waters close to land tend to be sources of CO 2 due to respiration of terrestrial organic carbon and lateral transport of CO 2 -enriched waters from the adjacent inshore system, whereas the mid-to-outer shelf waters generally act as a sink of CO 2 due to increased primary production as a result of nutrient availability from upwelling and mixing across the shelf break (Walsh, 1991;Jiang et al., 2008). Moreover, coastal shelves show a temporal variability in sea-air CO 2 fluxes as a result of variability in riverine input, phytoplankton blooms, ice formation and melting, and sea surface temperature (e.g., Takahashi et al., 2009;Rysgaard et al., 2012;. For example, Borges et al. (2005) documented seasonal changes of surface seawater CO 2 partial pressure (pCO 2 ) as high as 400 matm in the coastal shelves as a result of riverine input and other coastal processes.
In high-latitude coastal seas, changes in sea-ice coverage (Notz and Stroeve, 2016), water temperature (Shepherd, 2016), riverine input (Carmack et al., 2016), glacial melt input (Pierre et al., 2019), and primary production (Arrigo and van Dijken, 2015) will have a significant impact on the distribution and magnitude of pCO 2 and thus the strength of atmospheric CO 2 uptake. Given these potential changes, as well as the spatial and temporal variability in pCO 2 , studying the spatiotemporal patterns of pCO 2 independently in the Arctic and sub-Arctic coastal shelves is important to improve our knowledge of the main drivers affecting pCO 2 and to increase the reliability of sea-air CO 2 flux budgets. The Hudson Bay system (including James Bay, Foxe Strait, Ungava Bay, and Hudson Strait) is a sub-Arctic/Arctic continental shelf sea that receives nearly one-third of Canada's river discharge from more than 42 rivers (Déry et al., 2011). Hudson Bay and its watershed have been recognized as highly sensitive to climate variability relative to other Arctic shelves (Arctic Climate Impact Assessment, 2005;Markus et al., 2009). For instance, several studies have documented changes in timing and volume of river discharge (Déry et al., 2011(Déry et al., , 2018, the timing of sea-ice formation and melting (Saucier et al., 2004;Hochheim et al., 2011), and biological production (Sakshaug, 2004;Sibert et al., 2011). Permafrost melting due to climate change may also profoundly affect the quality, quantity, and composition of organic matter in Hudson Bay (Gough and Leung, 2002;Mallory et al., 2010;Godin et al., 2017). These potential changes could have a drastic impact on sea-air CO 2 fluxes in Hudson Bay.
Although included in several global coastal ocean CO 2 budgets (e.g., Chen and Borges, 2009;Chen et al., 2013), Hudson Bay has not been included in other global and Arctic-specific budgets (Cai et al., 2006;Bates and Mathis, 2009;Laruelle et al., 2014). This limited inclusion is due to a lack of observations on spatial and temporal scales. Only two studies investigating sea-air CO 2 fluxes in Hudson Bay have been published (Else et al., 2008a(Else et al., , 2008b. Using pCO 2 measurements at discrete stations collected in September and October, Else et al. (2008a) found that the bay acts as a minor sink in October, while Else et al. (2008b) applied a remote sensing approach to determine that Hudson Bay acts as a weak CO 2 source during August and September. The limited temporal and spatial extent of those observations means that we have no reliable estimate of sea-air CO 2 exchange in Hudson Bay over the complete open water period. Therefore, our baseline estimates in this large inland sea need to be extended, especially into the spring and early summer seasons, when the ecosystem is likely subject to rapid and pronounced physical and biological changes (e.g., Else et al., 2012).
Here, we report the first spring and early summer measurements of pCO 2 and estimates of sea-air CO 2 fluxes in Hudson Bay and Hudson Strait. The main objectives of this study were to (1) examine the spatial and temporal patterns of pCO 2 ; (2) identify the key processes driving the variability in the observed pCO 2 ; and (3) estimate the seaair CO 2 fluxes during the 2018 spring and early summer seasons in Hudson Bay and Hudson Strait.

Study area
Hudson Bay is a large shallow inland sea located in central Canada, with an average depth of 125 m and a maximum depth of 250 m, while adjacent areas of Foxe and Hudson Straits can reach depths of 400 m (Prinsenberg, 1986;Saucier et al., 2004). Hudson Bay transitions from complete ice cover in the winter to open water in summer, which adds about 140-160 cm of freshwater when sea-ice melts from May to July (Prinsenberg, 1988). Freshwater from river runoff contributes about 85 cm, distributed around the perimeter of the bay (Prinsenberg, 1986), and largely constrained near the coast (Granskog et al., 2007). Peak river runoff occurs in May, whereas the minimum occurs in March, with a notable change in recent decades toward higher winter discharge for hydroelectric power generation (Déry et al., 2011). This massive seasonal freshwater input to the bay leads to a strong stratification in the upper surface layer Burt et al., 2016;Ahmed et al., 2020), which impacts the mixing of heat and nutrients between the surface and deep waters (Drinkwater and Jones, 1987;Prinsenberg, 1988). Water flowing out of the Canadian Arctic Archipelago enters Hudson Bay at its northwestern end through Foxe Strait, while the northern portion of Hudson Strait allows water inflow of both Atlantic and Arctic origin from the Labrador Sea. Water exits Hudson Bay via the southern portion of Hudson Strait (Ingram and Prinsenberg, 1998). The general circulation in the bay is cyclonic (counterclockwise), which can be attributed to wind-driven transport and estuarine circulation buoyancy resulting from river runoff and ice melt/formation (Prinsenberg, 1986;Saucier et al., 2004). In this study, we focus on Hudson Bay (approximately 0.8 Â 10 6 km 2 ) and Hudson Strait (approximately 0.2 Â 10 6 km 2 ), where we conducted the majority of our field measurements.

Data and methods
We collected our field measurements during a 6-week survey from May 25 to July 13, 2018, onboard the research icebreaker CCGS Amundsen (Figure 1). This study was part of the BaySys project (see other papers in this special feature of Elementa), a comprehensive study examining the relative influences of climate change and hydroelectric regulation on Hudson Bay's physical, chemical, and biological systems. The field program was the first-ever research cruise in Hudson Bay in the spring and early summer, at a time when sea-ice cover was still partially intact, and sea-ice melt and river runoff are near their peaks.

Surface water measurements
Surface water pCO 2 was measured during the cruise using an underway pCO 2 system (General Oceanics model 8050; Pierrot et al., 2009). This automated pCO 2 system continuously sampled water through a shower-type equilibrator system at a nominal flow rate of 2.5 L min -1 from a through-hull inlet located approximately at 7 m beneath the water surface. The precision of the underway pCO 2 system is 0.1 matm, with an overall accuracy estimated at 2 matm (Pierrot et al., 2009). The system was calibrated every 8 h using four certified standard gases (ultra-high purity N 2 as a zero gas and three CO 2 /air mixtures at concentrations of 295.8, 448.6, and 551.7 ppm) traceable to World Meteorological Organization standards. The system was very stable during our sampling time, with drift typically less than 0.2 matm between calibrations, contributing to uncertainty in our measurements less than 0.5%. We calculated the seawater pCO 2 using the mole fraction of CO 2 (xCO 2 ) based on Dickson et al. (2007) and following the recommendation of Pierrot et al. (2009)   calculations and xCO 2 corrections using the measured standards.
Using discrete water samples, Ahmed et al. (2020) found that the sea-ice melt created a shallow, stratified freshwater layer that played a major role in diluting the surface pCO 2 across Hudson Bay. To correct for vertical pCO 2 gradients resulting from freshwater stratification during our study, we applied a linear regression that predicts the pCO 2 gradient as a function of observed pCO 2 (fully described in Ahmed et al., 2020) to about 45% of the total underway measurements (those measured 1-5 weeks after ice breakup with wind speed less than 6 m s -1 ). For pCO 2 in areas measured less than 1 week after ice breakup, stratification was considered to be too strong for confident corrections, and those records (about 10% of the total measurements) were discarded. For pCO 2 in areas measured more than 5 weeks after ice breakup or areas 1-5 weeks after ice breakup with wind speeds greater than 6 m s -1 , no correction was applied (45% of the measurements) as the near-surface vertical pCO 2 and salinity gradients were not likely significant under those conditions (Ahmed et al., 2020).
The underway pCO 2 system has flow-through temperature, salinity, and oxygen sensors (Idronaut model Ocean Seven 315), which were factory-calibrated before deployment. However, we opted to use the temperature and salinity recorded by a separate thermosalinograph (TSG) that draws water from the same intake line but is installed closer to the ship inlet compared to the underway pCO 2 system (distance difference between two instruments in the engine room was about 3 m). A chlorophyll a (Chl a) fluorescence sensor (WestStar fluorometer with an accuracy of 0.03 mg L -1 ) and a fluorescing dissolved organic matter (FDOM) sensor (WETLabs ECOFLD; 370 nm excitation, 460 nm emission; Guéguen et al., 2015) were also installed on the TSG line. Both the FDOM and Chl a sensors were also factory-calibrated before deployment. Measurements of surface nitrate (NO 3 ) concentrations and surface photosynthetic active radiation (SPAR) across Hudson Bay were collected using the Seabird Satlantic MBARI-ISUS V3 nitrate and the Biospherical Instruments SPAR (QCR 2200) sensors, respectively, which were mounted on the ship's rosette. Additional details about sensors, calibrations, and data processing of the TSG and rosette data are reported in Amundsen Science Data Collection (2019).
Despite the proximity of the pCO 2 system to the inlet (<5 m), we observed a slight increase in the recorded water temperature in the equilibrator relative to the onstation measurements made by the ship's conductivitytemperature-depth (CTD) sensor on the rosette. Therefore, we corrected the pCO 2 measurements for this temperature increase following the procedure described by Takahashi et al. (1993). Using the maximum uncertainty in the in situ seawater temperature measurements through the pCO 2 temperature correction, we estimate that the total uncertainty in the final pCO 2 measurements is less than 2%. We processed the underway data to remove any measurements collected when the system was undergoing cleaning, calibration, or restricted water flows from the intake line (i.e., water flow <1.5 L m -1 ) during icebreaking operations. Further details about the data processing of the underway pCO 2 system are described in .

Atmospheric measurements
Atmospheric CO 2 mixing ratio was measured every minute using a closed-path infrared gas analyzer (Li-COR model LI-7000). The analyzer was installed in a container on the foredeck and connected to a gas sampling system that drew air from an inlet located on the ship's foredeck meteorological tower at a height of 17 m above the sea surface. The LI-7000 was calibrated twice a day with two certified gas standards traceable to World Meteorological Organization standards, and the mixing ratio measurements were converted into pCO 2atm following Dickson et al. (2007). Wind speed and air temperature were measured using a conventional propeller anemometer (RM Young Co. model 15106MA) and HMP45C212 temperature sensor, respectively. These sensors were mounted on the meteorological tower at a height of approximately 16 m above sea level. Wind speed was corrected to a height of 10 m assuming a log-linear wind profile and a neutral surface layer (Stull, 1988). We discarded the data when instruments on the tower were cleaned and/or calibrated or when the wind was blowing from the stern of the vessel.
For logistical reasons, the foredeck tower and associated sensors were removed 4 weeks into the cruise (on July 2, 2018). As a result, we opted to utilize wind velocities at 10 m height and surface air temperature provided by the National Centers for Environmental Prediction North American Regional Reanalysis (NARR) to have complete measurements along the ship track. The NARR wind and air temperature data are available eight times daily (i.e., every 3 h) with a horizontal grid spacing of about 32 km Â 32 km (Mesinger et al., 2006). We clipped and resampled the NARR data using a bilinear algorithm to 1-km spatial resolution and then matched it with the ship track based on spatial location and closest estimate to our sampling time. The NARR and in situ measurements show good agreement over the period when the in situ sensors were deployed, especially for the wind speed data (Table  S1). The atmospheric pCO 2 concentrations from May 25 to July 2, 2018, ranged from 400 to 419 matm, with a mean of 408 + 2.8 (1s) matm (n ¼ 17,308).

Sea-ice data
Daily sea-ice observations along the ship track were obtained from the Advanced Microwave Scanning Radiometer 2 (AMSR2) passive microwave ice product (approximately 3-km spatial resolution) provided by the University of Bremen. The AMSR2 characteristics and the ARTIST Sea Ice algorithm used for calculating the daily ice concentrations are described in detail by Spreen et al. (2008). We used the spatial locations of our ship track to index the daily ice maps during our study time and then retrieve the sea-ice concentration from the closest temporally located observation. We compared the AMSR2 ice data with weekly ice charts from the Canadian Ice Service Digital Archive (CISDA) to ensure that the AMSR2 product did not underestimate ice concentrations along the ship track due to melt ponds on ice and/or weather conditions. The CISDA data are based on integration and analysis of data from satellite imagery (i.e., synthetic aperture radar and optical), weather and oceanographic information, and visual observations from ship and aircraft. Past studies have demonstrated the reliability of CIS data to monitor sea-ice cover change in the Canadian Arctic (Tivy et al., 2011) and its application to sea-air CO 2 exchange studies . Figure S1 shows the sea-ice observations obtained from AMSR2, CISDA, and ship observations recorded by the coast guard officers on the ship. Although the AMSR2 product showed relatively better agreement with the ship observations compared to the CISDA ( Figure S1), we noticed that the AMSR2 underestimated ice concentrations in a few locations in Hudson Bay, especially along the southwest and east coasts of Hudson Bay (June 27-30, 2018), where we observed more melt ponds on the ice. For these reasons and to make sure we accounted for any biases in sea-ice concentrations, we decided to use both AMSR2 and CISDA products as upper and lower estimates, respectively, when calculating the sea-air CO 2 fluxes in Hudson Bay. Figure 2 shows the spatial variations of sea-ice concentration, salinity, temperature, and pCO 2 throughout the study area. At the start of our cruise, we observed a large polynya on the northwest side of Hudson Bay, with higher ice concentrations (>50%) encountered mainly in the south of Hudson Strait and the central and eastern side of the bay (Figures 2a and S2). The northwest polynya (hereafter NW polynya) in Hudson Bay has been observed occasionally in winter and spring (e.g., Gough et al., 2004;Landy et al., 2017), resulting from strong offshore westerly winds that advect the existing ice cover eastward and enhance local ice production (Saucier et al., 2004). However, sea-ice melt was actively underway throughout our study period, especially close to the ice edge and coastal areas. The rapid change in the sea-ice cover is evident from maps of sea-ice concentrations throughout the bay at the start, middle, and end of our cruise ( Figure S2). At the end of our cruise, sea ice was mainly limited to the eastern side of the bay ( Figure S2).

Underway observations
The spatial patterns of salinity and temperature were associated with the spatial patterns of sea-ice concentrations ( Figure 2). Salinity ranged from minima of <20 close to the sea-ice edges and the river-dominated coastal region to a high value of 33 in offshore ice-free regions ( Figure 2b and Table 1). We observed higher surface seawater temperatures (>8 C) close to the Churchill and Nelson rivers, and a low value of surface temperature (about -1.7 C) close to the sea-ice edge ( Figure 2c and Table 1). Typically, colder (<1 C) and less saline waters (salinity between 10 and 25) were encountered in the southern and eastern side of the bay near the ice edge, while relatively warmer (>4 C) more saline (salinity from 15 to 29) waters were observed in the river-dominated coastal zone (Figure 2a and c). We also observed high Spatial variability of (a) sea-ice concentration using the AMSR2 ice data and underway measurements of (b) sea surface salinity, (c) sea surface temperature, and (d) surface seawater pCO 2 along the ship track in Hudson Bay from May 25 to July 13, 2018. pCO 2 ¼ CO 2 partial pressure; AMSR2 ¼ Advanced Microwave Scanning Radiometer 2. DOI: https:// doi.org/10.1525/elementa.2020.00130.f2 Ahmed et al: Variability of pCO 2 and sea-air CO 2 fluxes in Hudson Bay Art. 9(1) page 5 of 22 salinity (>30) associated with cold water (<0 C) in Hudson Strait and the southwest side of Baffin Island. The observed variability in salinity matches the general cyclonic (counterclockwise) circulation in Hudson Bay (Prinsenberg, 1986;Ingram and Prinsenberg, 1998), as the marine water entering Hudson Bay from Hudson Strait and Foxe Basin is usually diluted by sea-ice melt and river runoff before exiting the bay to the Labrador Sea via Hudson Strait. We calculated the area-weighted average and standard deviation of these observations based on a 1 km Â 1 km grid to minimize the effect of sampling biases on our calculation of the bay-wide flux resulting from long occupations of some station locations. This binning resulting in decreasing the number of ship observations on the waterside from about 19,500 to 7,817 and on the atmospheric side from about 17,500 to 5,580. Hereafter, all the statistics reported in this study are area-weighted and based on n ¼ 7,817 for the sea and n ¼ 5,580 for the atmosphere. A summary of the range, average, and variability of the observed underway measurements is available in Table  1. The pCO 2 ranged between 125 and 650 matm ( Figure  2d, Table 1), with the lowest measurements observed close to the ice edge on the east side of the bay, the NW side of the bay, and in Hudson Strait (125-280 matm), while the highest pCO 2 values were mainly observed in ice-covered waters (380-550 matm) and close to the Churchill and Nelson Estuaries (390-465 matm). We calculated the average of pCO 2 measurements across Hudson Bay to be 317 + 61 matm, which is significantly lower (p ¼ 0) than the average of atmospheric pCO 2 (408 + 3 matm), indicating that overall Hudson Bay was undersaturated in pCO 2 relative to the atmosphere.
The measured pCO 2 showed no clear evidence of temporal changes during the study period, although it was clearly correlated in various locations with variations in salinity and temperature (Figure 3a, b). This result suggests that the observed variability was linked to the ship's spatial proximity to sea-ice or riverine input rather than the timing of sampling. However, we observed a relatively narrow range of salinity and temperature values (cooler, more saline water) during the first 2 weeks (May 30 to June 15), with higher variability and frequency of warmer, less saline water thereafter. This pattern is consistent with the heavy sampling in the NW polynya at the beginning of the cruise and increasing freshwater abundance as the ship entered the southern coastal domain where most of the runoff enters the bay. We also observed an increase in pCO 2 along the SW coast, associated with a decrease in salinity, an increase in temperature, and low sea-ice concentrations. On the east coast, we observed a decrease in pCO 2 values associated with low temperature and salinity. Low pCO 2 values were also observed in Hudson Strait associated with relatively higher temperatures and salinity. Additionally, the pCO 2 values showed some correlation with variability in sea-ice concentrations (Figure 3c). For instance, we observed that the high spikes of pCO 2 were usually associated with a higher sea-ice concentration (>80%), whereas the lowest pCO 2 observations were mainly associated with lower sea-ice concentration (<40%; Figure 3c).

Processes affecting surface pCO 2 variability 4.2.1. Seawater temperature
The influence of seawater temperature on pCO 2 is governed by the thermodynamic equilibrium of the inorganic carbon system (Takahashi et al., 1993). For example, if pCO 2 is only controlled by temperature change with no other physical or biogeochemical processes playing a role, the pCO 2 in water will double for every 16 C increase in temperature (Takahashi et al., 1993). The weighted-area average of seawater temperature observed during our Calculated sea-air CO 2 fluxes using the AMSR2 sea-ice data and NARR wind data in Hudson Bay.
Art. 9(1) page 6 of 22 Ahmed et al: Variability of pCO 2 and sea-air CO 2 fluxes in Hudson Bay sampling period was 0.8 C + 1.7 C, with values ranging from -1.8 C to 11.1 C ( Table 1). In order to assess the relative importance of the competing effects of thermal (i.e., temperature) and nonthermal (i.e., biological activity and/or physical water mixing) effects on changes in pCO 2 values, we calculated the temperature-normalized pCO 2 (pCO 2nonthermal ) and the effects of temperature on the pCO 2 (pCO 2thermal ) using the decomposition used by Takahashi et al. (2002). More details about the thermal and nonthermal pCO 2 calculations are described in Text S1.
In Figure 4, we show the temporal variability of pCO 2nonthermal and pCO 2thermal . The high temporal variability of pCO 2nonthermal compared to pCO 2thermal in the first 3 weeks of our sampling reflects the relatively stable temperatures and significant contribution of nonthermal processes (i.e., biological activity, water mixing, and sea-air CO 2 exchange) to the observed variations in pCO 2 . For example, the high pCO 2nonthermal (>450 matm) in the NW region could be attributed to the inflow of pCO 2 -rich water from rivers, mixing with subsurface waters, and/or respiratory CO 2 production.
Starting at about the 4th week of our study period (June 27), during the shipping transit along the SW coast of the bay, pCO 2thermal became more variable (Figure 4). This increased variability is not surprising as we measured higher water temperatures during this time (Figure 3b), which were mostly linked to an increase in the air temperature with time ( Figure S3a) and/or the presence of high-temperature riverine waters from Nelson and Churchill Rivers (Figure 2c). The ratio of thermal to nonthermal effects on pCO 2 (DpCO 2ratio ) shown in Figure 4 illustrates that the surface water pCO 2 in Hudson Bay and Hudson Strait was controlled mainly by nonthermal processes at the start of the sampling period, with the thermal effects more dominant at the end of our study.
As pCO 2nonthermal varies substantially through the entire study, we show the spatial variability of pCO 2nonthermal in Figure 5a. The average of pCO 2nonthermal is 319 + 66 matm, with a somewhat smaller range (124-548 matm) than the observed pCO 2 values ( Table 1). The largest differences between pCO 2nonthermal and observed pCO 2 were close to river mouths. For example, the differences between pCO 2nonthermal values and the observed pCO 2 measurements in the Nelson and Churchill Estuaries on average were about 55 matm (Figures 5a and 2d), indicating that mixing with relatively warm water introduced by Nelson and Churchill Rivers was contributing to the observed high pCO 2 .

Biological production and FDOM
The distributions of oxygen saturation, Chl a, and FDOM (Figure 5b-d) indicated that for the most part, biological processes and terrestrial organic matter were not a major driver of pCO 2 distribution in Hudson Bay at the time of sampling. This interpretation is supported by the low Pearson's correlation coefficients (r < .2; p < .05) among pCO 2 , Chl a, and FDOM across the study area (see Section 4.3). However, there were some notable, localized exceptions. For example, the FDOM signal was relatively high close to the coast, especially in the SW region of Hudson Bay (Figures 5b and S5), associated with high terrestrial input in this region (Guéguen et al., 2011), while low FDOM signals (<700 counts) were observed in the offshore icefree water and the NW polynya, likely due to photobleaching in the clearer waters (Dainard et al., 2015). The highest FDOM values (up to 9,000 counts) in our study area were measured east of the Nelson Estuary close to sediment-laden sea ice. Similar to our observations, Granskog et al. (2007) and Guéguen et al. (2011) found high chromophoric dissolved organic matter and  FDOM values near the coast in Hudson Bay due to terrestrial input by rivers. Overall, we did not observe any consistent relationship between FDOM and pCO 2 in our study area, indicating that FDOM played a minor role in pCO 2 variability within Hudson Bay during the study period. The average of Chl a during our sampling time was 1.1 + 1.2 mg L -1 , with minimum and maximum observations of 0 mg L -1 and approximately 10 mg L -1 , respectively ( Table 1). O 2 saturation averaged 111 + 9% with a range varying from 90% to 199% across the study area (Table  1). Typically, the response of primary production to ice melt during the spring varied both spatially and temporally as a result of local differences in nutrient supply (e.g., Tremblay et al., 2015). We measured relatively high surface Chl a concentration (>4 mg L -1 ; Figure 5d) in Hudson Strait (especially in our transit at the end of the cruise) and in many locations in Hudson Bay (particularly close to the Nelson and Churchill Rivers and in the NW Polynya), associated with oxygen supersaturation (>120%; Figure  5c) and depleted nitrate (NO 3 ) concentrations (<1 mmol m -3 ; Figure S4), but variable pCO 2 values (Figure 2d). Because dissolved O 2 equilibrates much more quickly with the atmosphere than dissolved CO 2 , we cannot use O 2 saturation to rule out primary production. However, we can use the observed O 2 supersaturation as evidence of primary production. For example, the strong correlation among low pCO 2 (< 350 matm), high Chl a, and O 2 supersaturation measured in the NW polynya and Hudson Strait, where the measurements were made before significant gas exchange (i.e., the sea ice was actively melting; Figure 2a), suggest that the O 2 supersaturation is related to ice-associated blooms and that primary production played a significant role in lowering the pCO 2 values in these locations.
By contrast, in the Nelson Estuary, we observed a high Chl a concentration associated with high pCO 2 values (>440 matm), suggesting that the relatively high productivity in this area was not large enough to overcome other processes acting to increase pCO 2 . This is not surprising as the Nelson River has the highest annual discharge to Hudson Bay (Déry et al., 2011), with large loads of organic compounds and nutrients (e.g., Granskog et al., 2007), which potentially can sustain both the high biological production and the respiration that generates high pCO 2 values. Indeed, Demarty et al. (2009) observed the annual average of pCO 2 in the Kettle generating station reservoir located on the Nelson River (about 150 km upstream of Hudson Bay) to vary between 550 and 650 matm, showing the highly supersaturated nature of this river that might not be completely offset with high biological productivity. Furthermore, we observed variable Chl a concentrations close to the ice edge and under the ice, with values varying from less than 1-4 mg L -1 , and associated with O 2 saturation varying from saturated to supersaturated values. Similar to our high Chl a observations near the coast, Roff and Legendre (1986) observed higher primary production in these waters compared to offshore waters in Hudson Bay and attributed differences to nutrient replenishment of the upper productive layer through coastal resuspension processes (Griffiths et al., 1981), internal waves (Ingram et al., 1989), or freshwater plume dynamics . As well, Kuzyk et al. (2010) found that river runoff supports new primary production around the margin of the bay, where river flow is constrained. In contrast, the low surface Chl a (<1 mg L -1 ) in many locations in Hudson Bay (Figure 5d), suggests either low productivity or we sampled either before or after the spring bloom. Recently, Matthes et al. (2020) observed that 30% of the annual primary production occurs during the ice melt season (May to June) with stratification likely playing a significant role in delaying the bloom onset into the summer months. That said, our findings of low biological production in many areas of Hudson Bay are consistent with several previous studies that described Hudson Bay as an oligotrophic basin characterized by relatively weak biological production due to a short ice-free season and strong stratification (e.g., Legendre et al., 1996;Ferland et al., 2011;Sibert et al., 2011;Lapoussiere et al., 2013).
Although we cannot quantify whether deep primary production affects our surface pCO 2 measurements, it seems unlikely to have an impact given the depth of the deep Chl a maximum (which varied from 20 to 40 m based on similarities in CTD profiles across the study area) and its position beneath the strong and persistent mixed layer that exists in the region during spring (Ahmed et al., 2020). Similarly, Ferland et al. (2011) observed deep Chl a maxima at 40-60 m during summer and noted a negative correlation between primary production and stratification strength in the upper water column across Hudson Bay. We therefore conclude that biological production had a limited effect on the observed pCO 2 values during our sampling period, except in the NW polynya and Hudson Strait where primary production may have contributed to the observed pCO 2 undersaturation.

Ice melt and riverine input
To study the impact of freshwater from sea-ice melt and river runoff on pCO 2 variability in Hudson Bay, we plotted temperature versus salinity visualized with pCO 2 values throughout the study area ( Figure 6). Past studies have shown that sea-ice meltwater is initially undersaturated in pCO 2 (e.g., Geilfus et al., 2015) with low seawater temperature and salinity, whereas Arctic river waters tend to be relatively warmer and often supersaturated in pCO 2 as a result of degrading terrestrial organic carbon (e.g., Semiletov et al., 2011;Denfeld et al., 2013).
We measured low surface salinity seawater (i.e., <28) associated with colder water (<1 C) and low pCO 2 (<350 matm) in many locations close to the ice edge in Hudson Bay (Figures 6 and 2b and d), indicating that sea-ice melt contributed to pCO 2 undersaturation in these regions. We also measured high temperature (>5 C) associated with relatively low salinity observations (<29) and high pCO 2 (>420 matm) in some of the coastal regions in Hudson Bay (especially near the Nelson and Churchill Estuaries; Figures 6 and 2b and d), suggesting that river water contributed to pCO 2 supersaturation in these regions. Similarly, Capelle et al. (2020) noted that the inputs of terrestrial organic matter into Hudson Bay were a significant driver of CO 2 evasion in the coastal surface waters. Interestingly, the observations from the Churchill Estuary were warmer (0.5-11.2 C) and had relatively lower pCO 2 (averaged 350 matm with a range of 230-480 matm) compared to the temperature (0.9-8.3 C) and pCO 2 (averaged 416 matm with a range of 370-470 matm) observations recorded in the Nelson Estuary ( Figure 6). These differences between both estuaries could be attributed to differences in watershed chemistry, associated with the different physiographic regions, which likely influences the biology and chemistry of the two rivers (Rosenberg et al., 2005); river discharge (Déry et al., 2016), which is expected to influence carbon loading and marine biogeochemistry; and/or atmospheric conditions (e.g., air temperature and solar radiation) prior to sampling. For example, water and air temperatures in the Churchill Estuary were about 4 C and 9 C warmer, respectively, than in the Nelson Estuary during sampling (Figures 6 and S3a). Also, the Churchill River has higher concentrations of dissolved organic carbon than the Nelson River (Mundy et al., 2010), which is consistent with the Churchill watershed having a lower carbonate rock content and a higher permafrost abundance than the Nelson watershed.
Results from past studies in Hudson Bay are consistent with our observations, with low dissolved inorganic carbon associated with sea-ice melt (Burt et al., 2016), which typically results in low pCO 2 values (e.g., Ahmed et al., 2020). Additionally, our observations of high pCO 2 values close to river-dominated nearshore waters match the measurements of Else et al. (2008a). Overall, our results suggest that Hudson Bay was undersaturated in pCO 2 relative to the atmosphere during our study time as a result of the dominance of low-pCO 2 ice-melt waters, with few observations of supersaturated pCO 2 .
To more closely examine the effect of ice melt and riverine input on pCO 2 in Hudson Bay, we visualized the ship transects from ice-covered waters to coastal areas in the NW polynya and the Nelson and Churchill Estuaries (Figure 7). This plot shows typically high pCO 2 (>410 matm) in the coastal domain and low pCO 2 (<330 matm) observed mostly in the offshore waters. Figure 7a confirms that the Nelson and Churchill Rivers were responsible for the supersaturated pCO 2 during our sampling time, but their impact diminishes rapidly offshore toward the ice-melt-influenced waters. Similar to the transects in Nelson and Churchill Estuaries, off Chesterfield Inlet, we observed relatively high pCO 2 values that decreased rapidly as we moved offshore across the NW polynya ( Figure  7b). The undersaturated pCO 2 observations in offshore waters of the NW polynya are probably a result of ice melt and/or biological productivity promoted by the delivery of upwelled and/or riverine nutrients. Many of these processes have been observed in past studies of Arctic polynyas (e.g., Yager et al., 1995;Else et al., 2012).

Under-ice processes
We often observed high pCO 2 (>415 matm) associated with relatively high salinity (>31), cold (< -0.5 C) waters when the ship was navigating through the heavy ice cover (- Figure 7). The observed high pCO 2 measurements were also associated with high nitrate concentration ( Figure S4), suggesting that the under-ice phytoplankton bloom had not started yet when we were sampling at these locations. These high under-ice pCO 2 observations are likely remnants of biogeochemical and physical processes during the winter season. For example, community respiration has been observed to exceed photosynthesis during the winter as a result of low incoming solar radiation attenuated by sea ice and the overlying snowpack (Shadwick et al., 2011;Else et al., 2012). Another potential source of the observed high pCO 2 might be the expulsion of CO 2 -enriched brines into the water column during sea-ice formation and the retention of alkalinity by calcium carbonate crystals within the sea ice (Rysgaard et al., 2007;Dieckmann et al., 2008). Overall, Figure 7 shows both CO 2 supersaturation and undersaturation in the marginal ice zone, suggesting very high spatial variability in the distributions of sea-ice melt at these locations.

Upwelling
The observed relatively high salinity (>32) and low temperature (<0 C) on the east side of Southampton Island in Hudson Bay (Figure 7b) and the southwest side of Baffin Island indicate that upwelling contributed to high pCO 2 values in these regions (>420 matm). Upwelling is a winddriven process that brings the deep, cold, and usually high pCO 2 water to the surface, often resulting in oceanic CO 2 outgassing. Several past studies (e.g., Prinsenberg, 1986;Kuzyk et al., 2010) have observed upwelling events in the NW region of Hudson Bay during the summer and fall seasons, suggesting that upwelling in this area may be a long-lasting event rather than an episodic phenomenon. Therefore, although the pCO 2 supersaturation we observed in association with upwelling during our sampling was limited in spatial extent, it may persist through much of the open water season and possibly occur in other locations within Hudson Bay.

pCO 2 regression analysis
To explain statistically the key variables affecting pCO 2 in Hudson Bay, we performed a linear regression analysis between pCO 2 as the dependent variable and salinity, temperature, Chl a, and FDOM as independent or explanatory variables. We limited our spatial analysis to the ice-free regions (i.e., ice concentrations less than 10%) in Hudson Bay where the majority of our field measurements were collected. As well, we excluded the pCO 2 values impacted by the upwelling event in the NW region to ensure that our regression relationship is not biased by this spatially limited phenomenon. We found that salinity described the major variability in pCO 2 (Pearson's product-moment: r ¼ -.71; p < .05), followed by temperature (r ¼ .65; p < .05; Figure 8 and Table 2). Although Chl a and FDOM measurements showed a significant linear relationship with pCO 2 , the Pearson's correlation coefficients were low ( Table 2), showing that their impact on the pCO 2 variability is minor but statistically significant.
The regression relationship between salinity and pCO 2 , statistically significant at a 95% confidence level, was pCO 2 ¼ À14:16ð+0:17Þ Â Salinityþ753:58ð+5:22Þ: ð1Þ A root mean square error (RMSE) of 36 matm and relatively strong coefficient of determination (R 2 ¼ 0.51) indicate that this algorithm is explaining more than half of the variability in pCO 2 , which is a good result, given that this algorithm operates over a range of about 270 matm (i.e., pCO 2 ranged from 198 to 468 matm in the icefree regions). To examine the significance of other variables on pCO 2 , we added temperature, FDOM, and Chl a, sequentially, and selected the best model based on the lowest RMSE, Akaike information criterion, and Bayesian information criterion. The resulting multivariate relationship combined salinity, temperature, FDOM, and Chl a, with RMSE equal to 35 matm and R 2 of .55:

Sea-air CO 2 fluxes
To estimate CO 2 fluxes, we used a bulk flux equation: where F CO2 is the CO 2 flux (mmol CO 2 m -2 day -1 ), a is the solubility of CO 2 in seawater (mol m -3 atm -1 ) as a function of temperature and salinity (Weiss, 1974), and DpCO 2 is the difference between the measured surface seawater pCO 2 and the average atmospheric pCO 2 (in matm). We calculated sea-air CO 2 flux scaled by sea-ice concentration (C i , reported as a percentage) following the studies of Butterworth and Miller (2016) and Prytherch et al. (2017). The gas transfer velocity k (cm h -1 ) was calculated according to Wanninkhof (2014): where U 2 10n is the square of wind speed (m s -1 ) corrected to a height of 10 m, and Sc is the Schmidt number (unitless). Negative sea-air CO 2 fluxes represent atmospheric CO 2 uptake (oceanic sink), whereas positive fluxes represent oceanic CO 2 outgassing to the atmosphere (oceanic source).
As the shipboard wind speed data were not available after July 2, we estimated the CO 2 fluxes in Hudson Bay using both shipboard wind speed data and NARR wind data (Figure 9). The NARR wind speeds were usually underestimated compared to the in situ measurements (Figure 9a), with the highest wind speeds (>11 m s -1 ) during our study time observed in the NW polynya and the Nelson Estuary (Figure 10a). The shipboard wind speeds ranged up to 15 m s -1 with an average of 5.8 + 2.9 m s -1 (Table S1 and Figure 9a), whereas the NARR wind data  All coefficients are based on 10,000 observations and are significant at a 95% confidence interval.
Art. 9(1) page 12 of 22 Ahmed et al: Variability of pCO 2 and sea-air CO 2 fluxes in Hudson Bay showed a similar range with a lower spatial average of 4.4 + 2.9 m s -1 (Table S1 and Figure 10a). Consequently, CO 2 fluxes calculated using the shipboard wind data had a higher magnitude than those calculated using the NARR wind data (Figure 9b). The average CO 2 fluxes using the NARR and the ship wind data were -5.8 and -7.8 mmol CO 2 m -2 day -1 , respectively (Table S1). Therefore, we consider the CO 2 fluxes calculated based on the NARR wind data as lower-bound estimates of the oceanic CO 2 sink in Hudson Bay. The calculated DpCO 2 varied from about -283 matm to 242 matm, with an average of -91.5 + 60.9 matm (i.e., the surface waters were generally undersaturated, with respect to the atmosphere). During our study period in Hudson Bay, we observed supersaturated values of DpCO 2 (>10 matm) mainly close to the Nelson and Churchill Estuaries and along the southeast of Southampton Island ( Figure  S3b). Whereas, we observed strongly undersaturated values of DpCO 2 (< -100 matm) mostly in Hudson Strait and along the east side of the bay and in the NW polynya.
The high variation in the CO 2 fluxes shown in Figure  10b is attributed mostly to the variability in seawater pCO 2 resulting from the regional and dynamic variability of ice melt and riverine input (as discussed in Section 4.2). However, the high CO 2 uptake observed in the NW polynya was associated with high wind speeds (Figure 10a) over icefree waters (Figure 2a). To examine this spatial variability of sea-air CO 2 fluxes more closely, we visualized the ship transects in the Nelson and Churchill Estuaries and the NW polynya ( Figure S6). Overall, the strong CO 2 outgassing observed in these locations ( Figure S6) was correlated with the supersaturated pCO 2 (Figure 2d), resulting in positive DpCO 2 values ( Figure S3b) calculated for these ice-free areas.
Using all of the observations collected during our study, we found that Hudson Bay overall acted as a moderate to weak net sink for atmospheric CO 2 with an average CO 2 uptake of 5.1 (+9.3) mmol CO 2 m -2 day -1 ( Table  1). We estimated the uncertainty in the CO 2 flux to be 38% (+1.9) mmol CO 2 m -2 day -1 based on the cumulative uncertainties derived from the field observations, the gas exchange parameterization, and sea-ice concentration (see Text S2 for additional details on the uncertainty calculations). The moderate average sink of CO 2 that we calculated is a result of the dominance of uptake versus outgassing over the study area. Using the 1-km spatial bin analysis, 90% of the sampled area acted as a sink for CO 2 , while only 10% of the area acted as a source of CO 2 . Varying the source of sea-ice data did not seem to have a significant effect on these results; using the CISDA ice data in Equation 3, we calculated the average CO 2 uptake to be 4.8 (+9.0) mmol CO 2 m -2 day -1 .

CO 2 budget and comparisons to other Arctic coastal seas
While reporting the instantaneous CO 2 fluxes based on ship observations, understanding the hourly to daily trends and patterns in spring and early summer seasons is important (Figures 9b and 10b); otherwise, estimates might not be representative of the overall source-sink status of the bay. The problem with extrapolating ship observations to calculate CO 2 fluxes over space and time is that the processes influencing seawater pCO 2 are variable, hence biasing estimates to the sampling times and locations. For example, the biological activity and the magnitude of river inputs almost certainly change over time and space beyond our sampling survey. Although this spatial-temporal variability should be kept in mind, we expect that our resulting errors will be relatively low, given that Hudson Bay is known as an oligotrophic basin with low biological activity (e.g., Ferland et al., 2011) and that river discharge to Hudson Bay has relatively consistent seasonal hydrographs (Déry et al., 2016;. Assuming relatively consistent or slow change of freshwater input and biological production, we can extrapolate the average CO 2 fluxes over space and time to estimate the sea-air CO 2 fluxes in the open water areas in Hudson Bay. This approach of extrapolating ship CO 2 fluxes has been adopted by many studies in the Arctic regions to calculate the regional and annual CO 2 fluxes (e.g., Bates and Mathis, 2009, and references therein). Considering the sea-ice coverage during our sampling time, we estimated the area of open water in each month separately (i.e., May, June, and July in 2018) using the CIS-DA ice data. Using the average sink of -5.1 mmol CO 2 m -2 day -1 and the calculated open water areas, we estimated the total CO 2 fluxes for Hudson Bay in May, June, and July to be -0.3, -1.1, and -1.9 Tg C (10 12 g C), respectively, for a total net sea-air CO 2 flux of -3.3 (+ 1.2) Tg C.
Compared to other Arctic shelf regions, Hudson Bay during the spring and early summer seasons acts as a relatively moderate to weak CO 2 sink ( Table 3). Our calculated average uptake rates are similar to rates estimated for the Laptev and East Siberian Seas (Pipko et al., 2017), which are shelf seas similarly affected by strong river runoff. The calculated rates are also close to estimates for the Canadian Arctic Archipelago, which like Hudson Bay overall experiences low primary production . In contrast, the Chukchi Sea (which has high rates of primary and net community production during the icefree period) experiences much stronger CO 2 uptake (Bates, 2006). The high CO 2 uptake observed in the Kara sea (a shelf sea that receives a lot of riverine input) compared to Hudson Bay was more likely due to sampling close to seaice meltwater under high wind speeds .
Our estimates of CO 2 flux through May, June, and July complement past research in Hudson Bay that calculated CO 2 fluxes in late summer and autumn (Else et al., 2008a(Else et al., , 2008b. Else et al. (2008a) found that warm sea surface temperature in late summer (August) causes weak outgassing, which eventually reverts to a weak CO 2 sink as the surface waters cool through the fall. Combined with the moderate uptake we observed during spring and early summer, a fairly complete picture emerges that is consistent with the typical seasonal CO 2 cycle for high-latitude shelf seas described by Bates and Mathis (2009, and references therein). However, the temporal separation of the two studies (13 years) limits our confidence in any calculation of the total sea-air CO 2 flux throughout Hudson Bay over the entire open water season solely from the observational results.
Instead, Else et al. (2008b) found that the variability in seawater temperature is a good predictor for the pCO 2 variability during the summer and fall seasons in Hudson Bay. Therefore, to compute the total CO 2 fluxes during the entire 2018 open water season in Hudson Bay, we can use satellite-derived sea surface temperature and the average pCO 2 during our study (approximately 317 matm), corrected for temperature following the approach of Takahashi et al. (2002) to estimate pCO 2 during the 2018 summer and fall seasons. To estimate the average seawater temperature from August to October in 2018, we used the monthly Level-3 MODIS Aqua (https://oceancolor.gsfc.nasa.gov/), clipped the images to the study area, and resampled them to 1-km spatial resolution. The calculated average water temperatures for August, September, and October months were 5.6 C, 4.3 C, and 1.9 C, respectively. The resulting estimated pCO 2 values for August, September, and October were about 402, 380, and 344 matm, respectively. We calculated k (gas transfer velocity) using the NARR wind data and Equation 4 and estimated the sea-air CO 2 fluxes using Equation 3. The estimated Figure 10. Spatial variability of sea-air CO 2 fluxes in Hudson Bay. Spatial variability of (a) NARR wind speed and (b) seaair CO 2 fluxes based on NARR wind data. NARR ¼ National Centers for Environmental Prediction North American Regional Reanalysis. DOI: https://doi.org/10.1525/elementa.2020.00130.f10 Art. 9(1) page 14 of 22 Ahmed et al: Variability of pCO 2 and sea-air CO 2 fluxes in Hudson Bay CO 2 fluxes during the 2018 August, September, and October were -0.6, -1.2, and -2.2 Tg C, respectively. Adding these values to those we estimated for the spring and early summer seasons, we obtain a rough estimate of the total open water CO 2 uptake in Hudson Bay of -7.2 Tg C. While a complete annual budget should account for the sea-air CO 2 exchange during winter, through-ice fluxes are orders of magnitude smaller than open water fluxes and are not expected to have a large impact on the annual budget (Butterworth and Else, 2018). At an areal extent of about 1 Â 10 6 km 2 , Hudson Bay and Hudson Strait make up approximately 8% of the Arctic Ocean and its associated shelf seas (approximately 13 Â 10 6 km 2 ). Here, we estimated a net sea-air CO 2 flux of -7.2 Tg C during the open water season, suggesting that Hudson Bay and Hudson Strait are likely to account for 4-6% of the Arctic Ocean carbon sink, which previous estimates have shown is on the order of 118-180 Tg C year -1 (e.g., Manizza et al., 2019, and references therein). As a result, Hudson Bay is considered a slightly weaker sink than might be expected based on its surface area likely due to the influence of high-pCO 2 river water along the coastal regions. Yet, this contribution remains significant at the pan-Arctic scale and underscores the importance of including Hudson Bay in carbon budgets for the Arctic Ocean.

Conclusions
This study presents the first spring and early summer measurements of pCO 2 and estimates of sea-air CO 2 fluxes in Hudson Bay and Hudson Strait and thus fills up a major gap in our understanding of sea-air CO 2 fluxes in the Arctic/Sub-Arctic shelf seas. During our study, sea surface pCO 2 exhibited large spatial variability, with localized pCO 2 supersaturation in areas influenced by rivers and upwelling, as well as under the sea ice, and more generalized pCO 2 undersaturation in ice-melt-influenced waters. Overall, Hudson Bay was mostly undersaturated in pCO 2 . Temperature variations, and thereby their effect on pCO 2 , were small at the beginning of our study, when water mixing and biology activity were the primary drivers of pCO 2 variability. Along the southern coast, where most of the relatively warm riverine waters enter the bay, the thermal effect on pCO 2 became more important. The contribution of biological production to the variability in pCO 2 was minor but nonnegligible, especially in the NW polynya and Hudson Strait. In the open water areas, pCO 2 showed the highest correlations with salinity and temperature, with minor but statistically significant correlations with Chl a and FDOM. Our observations in the Nelson and Churchill Estuaries highlight the need to better understand variability in the biogeochemical makeup of rivers entering Hudson Bay and their impacts on the receiving marine system. Hudson Bay as a whole acted as a CO 2 sink with an average flux of about -5 mmol CO 2 m -2 day -1 (-3.3 Tg C) during the study, with an estimated uncertainty of 38% that corresponds to about +1.9 mmol CO 2 m -2 day -1 . The calculated sink in this study is significantly larger than previous estimates for Hudson Bay in late summer and autumn, underlining the significance of properly accounting for seasonal variability in Arctic shelf seas. Nevertheless, on average, the bay was a weaker CO 2 sink than Arctic shelves due to the strong influence of high-pCO 2 river water. We estimated the total CO 2 uptake in Hudson Bay during open water season to be -7.2 Tg C, contributing  about 4%-6% of the annual CO 2 flux estimated for the Arctic Ocean. Future work will attempt to use remotely sensed data and machine learning techniques to extend our observations and generate a more comprehensive estimate of regional sea-air CO 2 fluxes in Hudson Bay.

Data accessibility statement
The underway data presented in this article were uploaded to the Surface Ocean CO 2 Atlas (SOCAT) and will be available to the public in SOCAT Version 2021.

Supplemental files
The supplemental files for this article can be found as follows: Figure S1. Sea-ice concentration from daily AMSR2 and weekly CISDA sea-ice products and ship-recorded observations. Sea-ice concentration (%) from daily AMSR2 products (blue), weekly CISDA sea-ice products (black), and in situ observations recorded by the ship officers (red circles). AMSR2 ¼ Advanced Microwave Scanning Radiometer 2; CISDA ¼ Canadian Ice Service Digital Archive. Figure S2. Regional maps of changes in sea-ice cover over time. Sea-ice concentration (%, color-coded) at approximately the start (May 28), middle (June 18), and end (July 16) of our cruise in Hudson Bay. Sea-ice concentrations were provided by the Canadian Ice Service. Dates are year-month-day. Figure S3. Air temperature and sea-air CO 2 gradient in Hudson Bay. Spatial variability of (a) NARR air temperature and (b) the measured sea-air gradient of CO 2 partial pressure (DpCO 2 ) during our study time (May 25 to July 13, 2018). Figure S4. Nitrate and SPAR variability in Hudson Bay. Surface distribution of (a) nitrate (NO 3 ) concentrations and (b) surface photosynthetic active radiation (SPAR) across Hudson Bay during our study. The dots indicate locations of CTD casts, and the white area represents sea-ice concentration greater than 90%, as of July 9, 2018, based on weekly ice charts provided by the Canadian Ice Service. SPAR ¼ surface photosynthetic active radiation; CTD ¼ conductivity-temperature-depth. Figure S5. pCO 2 versus FDOM. Scatterplot of pCO 2 versus fluorescent dissolved organic matter (FDOM) in Hudson Bay during our study, with the measurements from the SW coast in red. pCO 2 ¼ CO 2 partial pressure. Figure S6. CO 2 fluxes in the Nelson Estuary and NW polynya. Spatial variability of sea-air CO 2 fluxes in (a) the Nelson and Churchill Estuaries with the sea-ice cover as of July 2, 2018, and (b) the NW polynya with the sea-ice cover as of June 16, 2018. NW polynya ¼ northwest polynya. Table S1. Shipboard and NARR match-up measurements of wind speed, air temperature, and sea-air CO 2 fluxes from May 31 to July 2, 2018. a NARR ¼ National Centers for Environmental Prediction North American Regional Reanalysis.
Text S2. Uncertainties of CO 2 flux estimation. their help during sampling; Tonya Burgers and Sebastian Luque for logistics and help during ship mobilization. They also thank Dr Elizabeth Jones and the second anonymous reviewer for their valuable comments and suggestions to improve the article. The ship thermosalinograph and conductivity-temperature-depth data used in this study were collected by systems onboard the Amundsen and made available by the Amundsen Science Program.

Funding
This study was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) Collaborative Research and Development (CRD) Program, Manitoba Hydro, and Fisheries and Oceans Canada.

Competing interests
The authors have no competing interests to declare. LAM serves as Elementa Associate Editor. She was not involved in the review of this article.

Author contributions
Contributed to the conception and design: MMMA, BGTE, TP.
Contributed to the acquisition of data: MMMA, DWC, BB, TP.
Contributed to the analysis and making the figures: MMMA.
Contributed to the interpretation of results: All authors.
Drafted the article: MMMA. Revised the article: All authors. Approved the submitted version for publication: All authors.