Changes in the Arctic Ocean Carbon Cycle With Diminishing Ice Cover

Abstract Less than three decades ago only a small fraction of the Arctic Ocean (AO) was ice free and then only for short periods. The ice cover kept sea surface pCO2 at levels lower relative to other ocean basins that have been exposed year round to ever increasing atmospheric levels. In this study, we evaluate sea surface pCO2 measurements collected over a 6‐year period along a fixed cruise track in the Canada Basin. The measurements show that mean pCO2 levels are significantly higher during low ice years. The pCO2 increase is likely driven by ocean surface heating and uptake of atmospheric CO2 with large interannual variability in the contributions of these processes. These findings suggest that increased ice‐free periods will further increase sea surface pCO2, reducing the Canada Basin's current role as a net sink of atmospheric CO2.


Introduction
The rapid loss of sea ice in the Arctic Ocean (AO) is well documented (Meier et al., 2014). Other changes in the AO are also becoming evident. Freshwater content is increasing due to sea ice melt and river runoff (e.g., Krishfield et al., 2014;Proshutinsky et al., 2009;Yamamoto-Kawai et al., 2009). Sea surface temperature has also increased (e.g., Perovich et al., 2011;Steele et al., 2008;Timmermans, 2015;Toole et al., 2010). This evolving physical environment is altering biological production (Arrigo & van Dijken, 2015;Bergeron & Tremblay, 2014) and food web structure (Hunt et al., 2014;Li et al., 2009;Søreide et al., 2010). The carbon cycle in the AO is intimately connected to these processes (Anderson & Macdonald, 2015), but it is not clear how carbon sources and sinks are changing in the AO and if they could affect CO 2 accumulation in the atmosphere and sea surface.
Sea surface pCO 2 is a key carbon cycle parameter because it is used to determine air-sea CO 2 fluxes for global carbon budgets and for understanding the rate of ocean acidification. Despite this, sea surface pCO 2 measurements in the AO are spatially and temporally sparse. While pCO 2 measurements date back decades (Kelley, 1970) and have continued on sporadic research cruises (Ahmed et al., 2019;Bates et al., 2006;Cai et al., 2010;Evans et al., 2015;Jutterström & Anderson, 2010;Robbins et al., 2013), measurements are mostly from AO shelf regions in the summer and fall due to limited access to more heavily ice-covered regions. Few studies have included repeat shipboard pCO 2 measurements from the AO's deep basins. The interior basins comprise~50% of the surface area of the AO (Bates et al., 2011) and, because of their reduced seasonal variability compared to AO coastal margins, might provide an earlier indicator of changes in pCO 2 as the ice-free sea surface is exposed to present-day atmospheric CO 2 levels.
While the AO is known to be a sink for atmospheric CO 2 (Arrigo et al., 2010;Bates et al., 2011;Islam et al., 2016Islam et al., , 2017Yasunaka et al., 2016Yasunaka et al., , 2018, its contribution to global air-sea CO 2 fluxes remains highly uncertain (Anderson & Macdonald, 2015). Estimated to uptake between 70 and 200 teragrams (Tg) carbon per year or 5-14% of the global uptake (Bates et al., 2011;Yasunaka et al., 2018), continued ice loss will make these estimates even more uncertain. Previous studies have found evidence that pCO 2 levels are changing in the oligotrophic AO basins and along the Chukchi shelf (Cai et al., 2010;Miller et al., 2014;Yasunaka et al., 2016Yasunaka et al., , 2018, and a new analysis by Ouyang et al. (2020) more clearly shows rising levels of pCO 2 in the Canada Basin since 1994. Increased CO 2 uptake by the AO will also accelerate ocean acidification, that is, driving a commensurate decrease in pH that increases calcium carbonate solubility (Robbins et al., 2013;Yamamoto-Kawai et al., 2009). In this manuscript, a shipboard pCO 2 time series collected from 2012-2017 in the AO's Canada Basin reveals that pCO 2 increases with decreasing ice concentration. We use a temporal reference, days since ice retreat (DSR), and a mass balance model to examine to what extent ice-dependent processes such as air-sea CO 2 fluxes, surface ocean warming, and biological production drive changes in sea surface pCO 2 after ice retreats.

Observations
Underway sea surface pCO 2 was measured on the Beaufort Gyre Observing System/Joint Ocean Ice Study (BGOS/JOIS) cruises on the CCGS Louis S. St-Laurent during 2012-2017. No pCO 2 measurements were made in 2015. The starting dates for the five~4 week cruises were 6 August 2012, 3 August 2013, 25 September 2014, 24 September 2016, and 8 September 2017. The pCO 2 was recorded using an infrared-gas equilibrator system (SUPER-CO 2 , Sunburst Sensors, LLC) located in the ship's lab. The instrument uses an infrared analyzer (LI-COR, LI-840A) and gas phase equilibrator (Liqui-Cel membrane contactor, Model #G453) as described in Hales et al. (2004). The equilibrator was connected directly to the ship's seawater line. Calibrations were automated using CO 2 gas standards and a zero CO 2 gas sample. Temperature was measured in the equilibrator and at the seawater intake (9 m depth), assumed equal to sea surface temperature (SST), as discussed below. The infrared analyzer CO 2 mole fraction was corrected to SST and converted to pCO 2 using 100% humidity at SST and local barometric pressure (Dickson et al., 2007). Some warming, usually <0.5°C, can occur enroute to the equilibrator so it is essential to correct the pCO 2 for this temperature change. If there were greater than~2°C differences between the equilibrator inlet temperature and SST, it was assumed seawater flow had stopped (e.g., due to ice clogging) and the pCO 2 was discarded during these periods. A flow meter was installed in 2016 to detect periods of low flow rate. The pCO 2 uncertainty is estimated to be ±5 μatm based on the reproducibility of the standards and baseline zero. CTD stations showed that the seawater intake was sometimes within the halocline (i.e., below the mixed layer), and this was also evident from the salinity and temperature variability recorded by the ship's thermosalinograph. These conditions were found mostly on the continental shelf during 2012 and 2013 due to high Mackenzie River outflow. This analysis focuses on the Canada Basin bounded by 155-130°W and 72-82°N where CTD stations consistently found that mixed layer depths were greater than the ship intake depth.
Air temperature, wind speed, wind direction, and barometric pressure were recorded by the ship's weather system. Mixed layer depths, defined as the depth where the density difference from the surface first exceeds 0.25 kg m −3 (Timmermans et al., 2012), were determined using temperature and salinity from~50 CTD casts occupied annually as part of the BGOS/JOIS cruises . Atmospheric pCO 2 was computed from the mole fraction of CO 2 measured at Alert, Nunavut, Canada, using data from the National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory (ESRL) (https:// www.esrl.noaa.gov). Sea ice concentration with daily, 12 km resolution was obtained from the French Research Institute for Exploration of the Sea (IFREMER) (http://cersat.ifremer.fr/oceanography-fromspace/our-domains-of-research/sea-ice) that provides data collected by the satellite-based Special Sensor Microwave Imager (SSM/I) and processed by the National Snow and Ice Data Center (https://nsidc.org).

Data Analysis and Modeling
In this study our goals are to determine if sea surface pCO 2 levels are related to interannual variability in ice concentration and to evaluate processes that might control pCO 2 under low (or no) ice conditions. To facilitate analysis of the spatially and temporally disparate shipboard pCO 2 and other data, Canada Basin data were gridded by identifying 20 × 20 km grid areas that contain data, and those data were then averaged as described in Evans et al. (2015). Underway pCO 2 measurements collected at different times within each grid cell were averaged, and an average time of the measurements was computed. The ship might spend as much as a day within a 20 × 20 km area because of variable activities, for example, mooring deployments, so these underway data were averaged together to avoid a spatial bias. The gridding routine also computed the standard deviation of data found within each grid cell as well as the number of observations. The average number of observations for each grid cell ranged from 17-23 for each cruise. As a temporal reference relative to the beginning of the open water period, we define "days since ice retreat" or DSR, as the difference between the pCO 2 measurement date and the day of ice retreat (DOR) for each grid cell. DOR is the day when ice concentration dropped below 15% in any grid cell (Steele & Dickinson, 2016), corresponding to the approximate uncertainty in ice concentration (Ivanova et al., 2015). The gridded data were also averaged together for each cruise (i.e., each year) to allow interannual comparisons of the mean physical and biogeochemical conditions. Further, we used these yearly values as input to the mass balance model described below. The mass balance model was used to examine how sea surface pCO 2 might change after ice retreats. Many processes contribute to sea surface pCO 2 variability in polar regions including biological production, heating and cooling, physical mixing and upwelling, ice melt and formation, and air-sea gas exchange. Sea surface warming and air-sea uptake are likely the most important factors for increasing pCO 2 in low ice areas in the Canada Basin (Cai et al., 2010;Else et al., 2013). The combined contributions of these two processes to sea surface pCO 2 variability were estimated using a dissolved inorganic carbon (DIC) mixed layer mass balance (Islam et al., 2017;Martz et al., 2009) as follows: where ΔDIC is the change in DIC for a time step Δt (1 hr in this study), F gasex is the air-sea CO 2 flux (e.g., in mmol m −2 day −1 ), MLD is the mixed layer depth, and ρ is seawater density. Warming (increasing SST) is accounted for in the equilibrium calculation as described below. F gasex was calculated using where k is the gas transfer velocity, K 0 is the CO 2 solubility, ΔpCO 2 is the pCO 2 difference between the sea surface and atmosphere, and f is the fraction of open water (Butterworth & Miller, 2016;Prytherch et al., 2017). In this case, f is set equal to 1 because the pCO 2 was modeled only after the day of ice retreat (DOR, defined above). F gasex is negative when there is a net uptake of CO 2 by the ocean from the atmosphere. We used the wind speed relationship in Wanninkhof (2014) to compute k, where ship wind speed was corrected to 10 m height. The average of second moments of wind speed (i.e., wind speed 2 ) was calculated as opposed to wind speed averages because short-term (<daily) variability in the winds is retained leading to higher gas transfer rates during periods of greater wind speed variability (Evans et al., 2015;Wanninkhof, 2014). A simple modification of Equation 1 makes it possible to examine potential contributions from net community production (NCP), where F NCP is the net uptake of CO 2 due to biological production (e.g., in mmol m −2 day −1 ). We used values from Ji et al. (2019) who measured NCP using O 2 isotope and O 2 /argon methods during the same BGOS cruises, excluding 2017. In the mass balance models, DIC was incremented for each time step with ΔDIC from Equation 1 (gas exchange only) or Equation 3 (gas exchange and biological production). The pCO 2 was then recalculated using the new DIC, SST, and salinity and a constant total alkalinity (A T ) in the CO 2 equilibrium program CO2sys (Pierrot et al., 2006). The A T was estimated using an A T -salinity relationship derived from bottle samples during the BGOS cruises (DeGrandpre et al., 2019;Yamamoto-Kawai et al., 2005). A T is a conservative property of seawater that does not change with temperature or air-sea CO 2 exchange (Millero, 2007). The time period of the model calculation is based on the maximum DSR for each cruise. For surface warming, SST was incremented equally for each time step from −1.5°C to the mean SST (Table 1) from DSR = 0 days until the end of the DSR period similar to Else et al. (2013). The initial under-ice condition was chosen to be the freezing point of seawater (−1.5°C) at a salinity of 27 and a seawater pCO 2 of 300 μatm all derived from extrapolation of the mean gridded data versus mean ice concentration to 100% ice concentration. All equilibrium (CO2sys) calculations used the Mehrbach et al. (1973) constants refit by Dickson and Millero (1987). Lastly, mean gridded values were used in the model calculations (Table 1 and Equations 1-3). A model that uses the evolution of ice concentration at each grid cell and includes other changing physical conditions, e.g. MLD or wind, was beyond the scope of this study.

Results and Discussion
The shipboard pCO 2 data collected on the Beaufort Shelf, Canada Basin, and eastern Chukchi Sea are shown in Figure 1 overlaid on % sea ice concentration. A large range of interannual variability in sea ice cover was observed during these cruises. In 2012, ice cover dropped to 3.4 million km 2 , the lowest level observed since the satellite record began in 1978. The minimum ice extent rebounded in 2013 and 2014 to~5.0 million km 2 . All 5 years ranked in the top 10 lowest minima (National Snow and Ice Data Center, Arctic Sea Ice News and Analysis, http://nsidc.org/arcticseaicenews). Sea surface pCO 2 was also highly variable spatially and interannually. Open water in the Canada Basin typically had higher pCO 2 levels than pCO 2 recorded in ice-covered areas; for example, compare 2012 (low ice) and 2014 (more ice) in Figure 1.
Using the mea levels for each cruise reveals a significant correlation with ice concentration (Figure 2). Sea surface pCO 2 is higher and closer to atmospheric saturation during years of low ice concentration and exceeded atmospheric pCO 2 (Table 1)  Anderson, 2010) and a new compilation of data since 1994 shows an increase in pCO 2 in the Canada Basin (Ouyang et al., 2020), no previous studies have documented an interannual connection with ice concentration like that shown in Figure 2. The mean pCO 2 values from each cruise reveal this relationship with ice concentration by consolidating the variability-the gridded pCO 2 data for each cruise (shown in supporting information Figure S1) cover the full range shown in Figure 2 obscuring interannual differences.
The absence of sea ice exposes the surface ocean to direct solar radiation and atmospheric heating and the subsequent warming increases pCO 2 . Air-sea exchange will also increase pCO 2 because, when the pCO 2 is lower than atmospheric levels, the AO will absorb CO 2 from the atmosphere (Figure 1 and Table 1). These mechanisms for increasing pCO 2 in the surface ocean suggest that pCO 2 is not only dependent upon the ice cover but also the duration of open water (Arrigo & van Dijken, 2015). Thus, the days since ice retreat (DSR) are used as a temporal reference, as defined above. We examine specific variables that might be important in driving the relationship between pCO 2 and sea ice concentration, including DSR, sea surface temperature (SST), MLD, wind speed, and NCP (Equations 1-3 and Table 1).
The gridded sea surface pCO 2 shows mostly increasing values after ice retreat with large intra and interannual variability (Figure 3a). Each year's observations appear to have a relatively consistent upward trajectory, except for 2016 and 2017 (discussed below), suggesting similar processes were at work in the open water area of the Canada Basin during each cruise. The pCO 2 observations have more scatter with increasing DSR, possibly due to different physical conditions during each cruise and each year. Some of the interannual variability may be due to differences in cruise timing, but there are significant contrasts among cruises conducted over similar periods. For example, in Figure 2, the mean pCO 2 from the August cruises (2012 and 2013) differs significantly, as do data from the late September cruises (2014 and 2017). These differences are also evident in the pCO 2 data in Figure 3a with lower values in 2013 and 2014 compared to 2012 and 2016, respectively.
The results of the mass balance model (Equation 1) assuming constant heating are shown in Figure 3a. The model encompasses the range of observed variability using the mean conditions in Table 1. Although there is  Modeled pCO 2 (dashed curves, labeled with each year and with the color matching the pCO 2 symbol data) were computed from the predicted change in pCO 2 due to air-sea exchange and increase in SST using Equation 1 and values in Table 1. Models were run to the maximum DSR recorded for each cruise period with panel (a) using a constant heating rate and panel (b) heating only until DSR = 60, then temperature is held constant. Initial pCO 2 before loss of ice was assumed to be 300 μatm (see section 2). The cruise start dates are indicated in parentheses in the legend. disagreement between observations and model curves during some years and some parts of the records, these results suggest that heating and/or gas exchange significantly contribute to the observed increase in pCO 2 with increasing DSR and that these processes are highly variable from year to year (broken down in Figure S3). In 2012, the model predicts that heating and gas exchange increased pCO 2 by~60 and 30 μatm, respectively; whereas in 2016, heating and gas exchange contributions were~15 and~65 μatm, respectively ( Figure S3). Also note that a smaller range of variability is observed and predicted for the years with shorter DSR periods (2013 and 2014) (Figures 3 and S3), where surface ocean warming and air-sea gas exchange were limited by the time of exposure to the atmosphere.
There are several other possible sources of variability and incorrect model assumptions that could contribute to differences between the model and observations. The sensitivity analysis in the Supporting Information shows that much of the variability within each cruise can be explained by varying input values over the observed ranges ( Figure S2 and Table 1). There are some other notable deviations, however. For example, during 2016 and 2017, the 2 years with extended DSRs, the pCO 2 levels decreased after~60 days (Figure 3). In these years, data were collected into October, at which point sea surface cooling is possible which would decrease the pCO 2 . Because the model employs a linear warming trend, it can only predict the mean change not the time-varying rate of warming or cooling. In Figure 3b, the model was modified so that all of the heating occurred before DSR = 60, slightly steepening the pCO 2 increase but not markedly improving the prediction. We did not introduce an arbitrary period of cooling because this would essentially be fitting to the data. There is also evidence that movement of the ship into different water masses over the course of each cruise may have contributed to the variability. For example, the drop in pCO 2 after 60 days in 2017 corresponds to a period of lower salinity (not shown). Also, employing DSR as a time variable in the model assumes that no air-sea exchange or warming occurred when sea ice concentration was greater than 15% which is clearly not true ( Figure S1). While we could not readily model gas exchange prior to DOR in this study, the heating contribution was small because no SST data exceeded 0.0°C from 0 to 10 days after ice retreat. An increase of −1.5°C, the approximate freezing point, to 0.0°C would increase pCO 2 bỹ 22 μatm. Steele and Dickinson (2016) also found that SST is typically <0°C at DOR. Ice melt and formation could change pCO 2 by diluting and concentrating DIC, respectively, and by altering the ratio of DIC and A T in the ice or brine (Cai et al., 2010;DeGrandpre et al., 2019;Rysgaard et al., 2007Rysgaard et al., , 2009. Ice melt occurs before DOR and could decrease pCO 2 levels (Cai et al., 2010). In fact, heating, gas exchange, and ice melt may have all played a role in determining the pre-DOR pCO 2 levels, contributing to the deviations between model and observations evident in Figure 3.
Lastly, we consider possible biological contributions to the observed variability. Biological drawdown of sea surface pCO 2 in the Canada Basin is predicted to be small because of the lack of nutrients in the stratified surface layer. The study by Cai et al. (2010) found no evidence of a biological DIC drawdown in the Canada Basin for waters >72°N. In a previous study in the Canada Basin, an NCP of~4.9 mmol O 2 m −2 day −1 offset a~15 μatm pCO 2 increase expected from atmospheric CO 2 uptake in low ice conditions (Islam et al., 2017). The Ji et al. (2019) NCP values, which ranged from 1.3 to 2.9 mmol O 2 m −2 day −1 from 2011-2016, confirm the very low productivity for the Canada Basin (Bates et al., 2005). For comparison, in the same paper, NCP in the eastern Chukchi Sea was estimated to range from 30 to 240 mmol m −2 day −1 . The Canada Basin rates have a relatively small effect on the modeled pCO 2 levels, estimated using Equation 3 ( Figure S2). In 2012, the pCO 2 would be reduced by 8-25 μatm by the end of the DSR period (60 days) over the range of NCP reported in Ji et al. (2019). These estimates assume the rate of NCP was constant over the DSR period. The O 2 /argon method integrates NCP over the residence time of O 2 in the mixed layer (10-30 days) (Kaiser et al., 2005) and the rates in Ji et al. (2019) varied by only 15-21% during the cruises. This range of variability would not significantly alter the modeled pCO 2 trajectory ( Figure S2 and Table S1).

Conclusions
This study reveals that loss of sea ice leads to large interannual changes in sea surface pCO 2 levels in the Canada Basin. Using the DSR as a temporal reference facilitated implementation of a time-dependent mass balance model to explore the underlying mechanisms that might control sea surface pCO 2 in open water conditions. Results from the model suggest that warming and air-sea uptake drive the pCO 2 toward atmospheric equilibrium to a varying extent (Figures 3 and S3). NCP is persistently low and can only account for a small portion of the observed interannual variability ( Figure S2 and Table S1). These results suggest that neither of two previously proposed scenarios, that is, that increased SST will largely negate uptake of atmospheric CO 2 (Else et al., 2013) or that air-sea gas exchange, will consistently dominate the increase in pCO 2 (Bates et al., 2006;Cai et al., 2010) (Figure S2). An important implication, as discussed by others (Cai et al., 2010), is that increased SSTs can reduce the uptake of atmospheric CO 2 by increasing the rate at which pCO 2 rises toward atmospheric equilibrium, decreasing the air-sea CO 2 gradient. The model simulation suggests that warming reduced the possible air-sea CO 2 flux by 55% in 2012, the year with the most warming ( Figure S3). As Arctic warming and open water periods increase, it is possible that, as observed in 2012, the lowest ice concentration year on record, pCO 2 will more frequently exceed atmospheric saturation. Over the period covered in this study any accumulation of CO 2 is hidden by the interannual variability of ice cover, but the study by Ouyang et al. (2020) over a 23 year period clearly reveals this expected CO 2 increase in the Canada Basin. Even with the insights provided in these new data sets, the future response of the AO carbon cycle to decreased seasonal ice cover remains highly uncertain and can only be understood by continued observations and modeling.