Underestimation of oceanic carbon uptake in the Arctic Ocean: Ice melt as predictor of the sea ice carbon pump

,

The sea ice carbon pump is considered to result mostly from three groups of processes: (i) sea ice growth or melt, which implies a freshwater flux (upward or downward) from the ocean to the ice, (ii) brine rejection, which proportionally decreases the uptake of solutes in sea ice, and (iii) active biogeochemical processes, which modify the alkalinity to DIC ratio in sea ice. Most, if not all, Earth System Models (ESMs) lack a representation of biogeochemical processes within sea ice and are 60 therefore unable to account for (ii) and (iii), but encompass (i) by dilution and concentration of tracers, similar to the handling of precipitation and evaporation. In the present study, we do not distinguish between (ii) and (iii) and instead consider that the carbon cycling in sea ice encompasses both aspects. We also consider our reference point (later referred to as CTRL) to be that of current ESMs, i.e. they include processes (i) but not processes (ii) and (iii).
Arctic sea ice extent and thickness have declined rapidly over the past decades at a rate of -83,000 km 2 yr −1 for September ice extent during the 1979-2018 period and with a decline in ice thickness by 65 % from 1975 to 2012 (Meredith et al., 2019).
This decline is expected to continue. Arctic amplification, a combination of positive feedbacks including summer albedo loss and changes in cloudiness, is leading to twice the rate of warming of the atmosphere compared to the global average (Meredith et al., 2019, Box 3.1). Increased "Atlantification" of the Eurasian Arctic Basin, characterized by a progression of Atlantic water masses into the Arctic seas, is contributing to amplified basal ice melt (Polyakov et al., 2017). These dynamic and 70 thermodynamic processes have direct impacts on sea ice seasonality and extent (Perovich and Richter-Menge, 2009) and icefree summers are predicted to happen within the next few decades (Overland and Wang, 2013;Notz and Community, 2020).
Yet, since sea ice extent in winter decreases slower than in summer, the seasonally ice-covered area is expanding. Such an amplified seasonality in sea ice may intensify the sea ice carbon pump, as sea ice forms in open water that had previously been perennially ice-covered. 75 We use two independent and complementary approaches to investigate the supplementary carbon flux in the Arctic Ocean.
We define the supplementary carbon flux ∆F as the fraction of the air-sea CO 2 flux that is solely due to the storage of carbon and alkalinity in ice. This term is quantified here as the difference in air-sea CO 2 flux between a reference situation where there is no ice-ocean carbon flux, i.e including aforementioned processes (i) but not (ii) nor (iii), and situations where ice-ocean carbon flux occurs, i.e including (i), (ii) and (iii). First, we combine a set of mathematical formulations to obtain an equation 80 that provides a theoretical framework for the description of the impact of alkalinity and DIC storage in sea ice on air-sea CO 2 fluxes. These theoretical considerations suggest that sea-ice melt and open-water fraction are the main drivers of an increased oceanic carbon uptake induced by storage of alkalinity and DIC in sea ice. Second, a simple parameterization of the presence of alkalinity and DIC inside the sea ice is implemented in a one-dimensional (1D) ocean model applied to different locations of the Arctic. A large set of sensitivity runs with this 1D model consolidates and expands on the role and importance of melt and 85 open-water-fraction and shows that the alkalinity-to-DIC ratio in sea ice plays a major role in the magnitude of the increased uptake. By forcing the model with a wide range of plausible ice conditions, we obtain a predictive linear relationship between annual ice melt and ice-induced annual supplementary carbon uptake (∆F). This relationship can be used to correct carbon uptake estimates from numerical models that do not account for carbon storage in ice. By applying the relationship to an Earth System Model (ESM) from the sixth phase of the Climate Model Intercomparison Project (CMIP6) ensemble, we show how 90 the impact of sea ice on carbon uptake may evolve under different future emission scenarios.
2 Theoretical Framework for Ice-Sea Carbon Flux and Induced Air-Sea CO 2 Uptake The impact of carbon storage in sea ice on the air-sea CO 2 flux is analyzed using differential equations that account for the impact of freezing and melting on surface water alkalinity and DIC. The air-sea flux is expressed as a function of sea surface pCO 2 , which depends on temperature, salinity, alkalinity and DIC. 95 We assume the flux of alkalinity and DIC between the sea ice and the underlying water to be proportional to the freshwater flux induced by freezing and melting of sea ice, F ice−sea F W (m s −1 ), and the concentration of alkalinity and DIC inside the ice.
The DIC and alkalinity concentrations are assumed to be homogeneous in the ice. The freshwater flux is positive (downward) for melting. The change in sea surface pCO 2 , written ∂pCO ice−sea 2 ∂t , resulting from the freshwater flux can then be expressed as where H 0 is the mixed layer depth (in m), considered constant for ease of interpretation; [Alk] ice and [DIC] ice are the concentrations of alkalinity and DIC inside sea ice (mmol m −3 ) and ∂pCO2 ∂Alk and ∂pCO2 ∂DIC are the fractional change of pCO 2 with alkalinity and DIC, respectively (µatm m 3 mmol −1 ). Note that ∂pCO2 ∂Alk and ∂pCO2 ∂DIC are generally non-linear.

105
The relation between the air-sea flux of CO 2 and seawater pCO 2 is where pCO 2 and pCO atm 2 refer to pCO 2 in surface seawater and atmosphere (µatm) resp., k g is the gas transfer velocity (m s −1 ), S CO2 is the CO 2 solubility (mol m −3 µatm −1 ) and λ is the fraction of open water (lead fraction, unitless; Ahmed et al., 2019). Here, the air-sea CO 2 flux is defined as positive downward.

110
The supplementary flux, ∆F t , is calculated as the difference between a case with carbon storage in ice, referred to as ICE, and a control (CTRL) case, where storage is not considered and ice growth or melt only leads to a freshwater exchange, i.e.
with pCO ICE 2 and pCO CT RL 2 the sea surface pCO 2 in the ICE and CTRL cases, respectively. In the rest of this manuscript, We assume that in both CTRL and ICE cases, sea surface pCO 2 experiences the same alterations due to biological processes and changes in temperature and salinity caused by vertical and horizontal mixing and air-sea-ice interactions. This assumption neglects the possibility that non-linearities of the carbonate system lead to differences in the impact of these processes on pCO 2 between the CTRL and ICE cases. Moreover, we assume that ∂pCO2 ∂DIC is constant. Calculations conducted with CO2SYS (Lewis and Wallace, 1998) based on mooring data located in the Beaufort Gyre (DeGrandpre et al., 2019, 78°N, 150°W) and our 120 model data (see Beaufort Gyre setup in Sect. 3.1) yield a coefficient of variation of ∂pCO2 ∂DIC of only 6 % and 5 %, respectively.
This supports the assumption of a constant ∂pCO2 ∂DIC over the range of expected DIC. These two assumptions are only used in this theoretical derivation, not in the numerical analysis.
The change of pCO 2 over time can be written as with the temperature and salinity contributions (the last two terms on the right-hand side) being identical in the ICE and CTRL cases.
The contributions from alkalinity and DIC can come from advection, diffusion, mixing, biological processes (production, respiration, remineralization), and air-sea or ice-ocean carbon fluxes. As already described, the ice-ocean carbon flux modifies the surface seawater pCO 2 , which in turn impacts the air-sea carbon flux. Here, the ice-ocean and air-sea carbon fluxes are the 130 two only processes that are not considered as identical between CTRL and ICE cases and are therefore the only two terms left when subtracting the equations for ∂pCO2 ∂t for the CTRL and ICE cases from each other. The following differential equation governing the evolution of ∆pCO i−c 2 can be derived (see details in the supplement): The solution to Eq. 3 is: where A(t) is a primitive of ∂pCO2 ∂t 1 H0 k g S CO2 λ and s is the variable of integration, with units of seconds. The primitive of a function can be calculated as its time integral plus an unknown constant α This yields a solution for the instantaneous difference in pCO 2 between CTRL and ICE scenarios. To retrieve the previously 140 defined supplementary carbon uptake, i.e., the cumulative air-sea CO 2 flux that is induced by the pCO 2 change, we can insert Eq. 4 into the left-hand side of Eq. 3 and integrate over a period T : A unique derivation to our knowledge, this formulation is composed of three main terms: g(t), which is a function of the concentration of alkalinity and DIC in the ice (Eq. 1); the freezing-melting flux F ice−sea F W ; and the more complicated exponential of the primitive, which contains the lead fraction λ in A(t). A(t) is an integral of the lead fraction and can be interpreted as keeping a memory of the evolution of the ice conditions.
The sign of g(t) determines the sign of ∆F t . Using realistic alkalinity and DIC values for the Arctic Ocean (e.g. and Revelle and alkalinity factors of 14 and -13.3 respectively, as in Takahashi et al. (1993)) yields a negative sign for g(t). It can be shown that the term between parentheses in Eq. 5 is always negative, meaning that for ice melt (F ice−ocean F W > 0), ∆F t is downward (uptake); the opposite is true for ice formation (cf. Supplementary materials). According to this formulation, the dependency of ∆F t on [Alk] ice and [DIC] ice is bi-linear due to the shape of g(t).
It is important to note that the gas transfer velocity and the CO 2 solubility, used in the primitive A(t), require no assumption 155 of shape or value. The gas transfer velocity k g can depend on the wind speed (e.g. Wanninkhof, 2014), on the wave slope (Bogucki et al., 2010) or on turbulence generated by ice drag and convection (Loose et al., 2014). Similarly, the CO 2 solubility could follow Weiss (1974) or any other expression.
One can calculate the solution numerically using the carbonate properties of seawater and sea ice (i.e., their alkalinity and DIC), the sea ice concentration, the ice-ocean freshwater flux, the gas transfer velocity (e.g., using Loose et al. (2014)) 160 and the CO 2 solubility (which depends on temperature and salinity, Weiss, 1974). The product of the Schmidt number and CO 2 solubility can reasonably be considered constant (Etcheto and Merlivat, 1988), therefore removing the dependency on temperature and salinity Wanninkhof (2014, their Eq. 6) and providing an even simpler form than proposed above.
In order to interpret the role of λ, its value can be constrained as follow. We assume that ice formation is associated with full ice cover (λ ≈ 0) and that melting occurs in open waters (λ ≈ 1). We will see that this is supported by ocean model output in This implies that the integrand in Eq. 5 is negligible during freezing and non-negligible during melting. Thus, ice formation has a relatively small contribution to the temporal integral of the supplementary carbon flux, while ice melting significantly increases the CO 2 flux. Since melting leads to uptake, according to the sign examination above, the net outcome of the supplementary carbon flux is uptake.

170
Note that if H 0 was assumed to be variable in time, it would remain inside the integrands on both sides of Eq. 5. The integrand is then likely to be small during the melting season when the mixed layer shoals, and larger during the freezing season, when λ is close to 0 and the integrand is already small. A variable mixed layer depth would therefore reinforce the already dominant influence of the melting season in the value of the supplementary carbon uptake.
3 Numerical ocean model 175 We implemented a parametrization of carbon storage and release by sea ice in a 1D ocean model, independent of the theoretical arguments in Sect. 2, to investigate its impact on the air-sea CO 2 flux in different regions of the Arctic Ocean. We do not use any of the results or assumptions from Section 2. By using a wide range of initial and forcing conditions derived from a realistic 3D model, a large ensemble of 1D simulations is generated to account for spatial and temporal variability in forcing conditions.
Analysis of the ensemble provides insights into the main drivers of the supplementary carbon uptake and allows us to derive a 180 formula to estimate the supplementary carbon flux in existing Earth System Model (ESM) simulations. Here we describe the 1D model set-up and forcings, as well as the ESM outputs used to project the evolution of the supplementary carbon flux in different scenarios.

One-dimensional ocean model
The 1D General Ocean Turbulence Model (GOTM, Burchard et al., 1999) is coupled to the Pelagic Interactions Scheme 185 for Carbon and Ecosystem Studies volume 2 (PISCES-v2, Aumont et al., 2015), specifically adapted to the Framework for Aquatic Biogeochemical Models (FABM, Bruggeman and Bolding, 2014). The vertical grid has fixed layer thicknesses, with a resolution of 1 m near the surface and increasing with depth (9 layers in the first 10 m, 24 layers in the first 100 m). Air-sea CO 2 flux is calculated by the model using values from Wanninkhof (2014) with 10 m wind speed. The carbonate chemistry in the model follows the OCMIP protocols (Orr, 1999).

190
The evolution of alkalinity and DIC in surface waters is parameterized by includes the advective and dispersive transport terms as well as dilution and concentration due to sea ice melting and freezing or due to precipitation and evaporation and Bio represents the biological sources and sinks of Alk and DIC, and F air−sea CO2 is the air-sea CO 2 flux. Preliminary runs showed that the biological terms have a similar impact on carbon uptake regardless of whether the carbonate system inside sea ice is represented or not, and thus yield a negligible impact on supplementary carbon 200 uptake (less than 1 % normalized difference. They were therefore deactivated for the ensemble runs to save computational effort. Surface forcings were prescribed from a 3D physical-biogeochemical-ice-ocean model based on NEMO-LIM-PISCES (Madec et al., 2017;Rousset et al., 2015;Aumont et al., 2015) for the North Atlantic, North Pacific, and Arctic Oceans, hereafter referred to as the NAPA model. The NAPA model, including the validation with observational data, is documented 205 in Zhang et al. (2020) and Zheng et al. (2021). In our application of this 3D model, the atmospheric forcing was obtained from ERA-5 reanalysis product (Hersbach et al., 2020) from 2014 to 2019. Outputs were written out daily, providing the necessary temporal resolution to capture sub-seasonal variability. We used the simulated ice concentration, latent and sensible heat fluxes, longwave and shortwave radiative fluxes, freshwater fluxes (due to ice melt-freeze and evaporation-precipitation), and momentum fluxes (due to wind and ice drift) at the top of the surface layer of the model, calculated as a weighted average between 210 open water and under ice conditions to force the 1D model. This methodology allows us to simulate the impact of sea ice in our 1D model without having to resort to a full ice component. Other inputs necessary for air-sea CO 2 flux include the wind speed and mean sea-level pressure from ERA5, as well as atmospheric pCO 2 from the Alert Station, Northwest Territories of Canada (Keeling et al., 2001).
In generating the ensemble of 1D simulations, every 10th horizontal grid cell of the NAPA domain was selected with the 215 following exceptions. Since our focus is on open-ocean conditions with a significant presence of sea ice, coastal locations with water depths shallower than 100 m as well as the Canadian Arctic Archipelago, Hudson Bay and the Baltic Sea were excluded.
Also excluded were grid cells with both ice melt and freezing rates of less than 0.1 m yr −1 . Given NAPA's average grid spacing of 12 km in the Arctic, every 10th grid cell leads to roughly one cell every 120 km for a total of 732 cells covering a wide range of ice conditions. For each of these locations, we ran the 1D model for six 1-year simulations starting on January 1st for the 220 years 2014 to 2019, with initial conditions from the NAPA model. Since the 1D model cannot explicitly represent horizontal advection, its solutions were nudged toward the properties simulated by the NAPA model with a timescale of 4 months for temperature and salinity and 1 year for alkalinity and DIC.
Based on the above setup, we systematically ran the 1D model in two configurations CTRL (no carbon in sea ice) and ICE (storage of carbon in sea ice). In both configurations, sea ice growth and melt generates a freshwater flux that concentrates or 225 dilutes tracers at the surface ocean. The runs are listed in Table 1. The supplementary carbon uptake ∆F m is calculated as the difference in annual air-sea CO 2 flux between the ICE (or ICE2) and the CTRL runs. We consider the air-sea CO2 flux in the CTRL run as the baseline, since it corresponds to the values reported by numerical models that do not account for the sea ice carbon pump. Potential predictors of the supplementary carbon flux are investigated including the net freezing-melting flux (the integral over a year of the freshwater flux between ice and ocean), the gross melting (freezing) flux which only accounts 230 for ice melt (formation), and the yearly integrated ice concentration (which ranges between 0 and 365). We bin the latter metric into 9 bins of equal size and applied a linear regression between gross annual ice melt and the supplementary carbon uptake for each of these bins, which can be considered different ice regimes.

Application to an Earth System Model
ESM output from the CMIP6 suite of models can be used to estimate the supplementary carbon flux in projected future 235 climate scenarios. We chose the ACCESS-ESM1.5 model (Ziehn et al., 2020) because it has a plausible simulation of sea ice (according to Notz and Community, 2020) and its monthly-averaged freshwater ice-ocean flux due to ice thermodynamics    How this amplification of the seasonal cycle of pCO 2 affects the seasonal air-sea CO 2 flux depends on the ice cover shown in Fig. 1a. According to the formulation in Eq. 2, almost complete ice cover (λ = 0) in winter results in an air-sea CO 2 flux close to 0 when pCO 2 is highest. Lower sea ice cover in summer allows for some air-sea gas exchange directly proportional to 265 the air-sea pCO 2 gradient. Integrated over a full seasonal cycle, the amplification of the pCO 2 cycle results in net oceanic CO 2 uptake added to the baseline. In the case of the Beaufort Gyre station location, averaged over both years, this supplementary uptake ∆F m amounts to 45.5 mmol C m −2 yr −1 for an alkalinity-to-DIC ratio of 1.80 (ICE) and 13.6 mmol C m −2 yr −1 for a ratio of 1.26 (ICE2), over a 3-fold difference. Note that these are low flux values relative to other oceans (usually higher than 1 mol C m −2 yr −1 ), mainly because of the ice cover.

270
The effect of carbon storage on the annual net CO 2 flux is explored more thoroughly for the Beaufort Gyre location by varying [Alk] ice from 340 to 700 mmol eq m −3 and [DIC] ice from 260 to 600 mmol m −3 (Fig. 2). The net CO 2 flux (Fig. 2a) and supplementary carbon uptake ∆F m (Fig. 2b) are strongly dependent on the alkalinity-to-DIC ratio in ice (white contours).
Notably, the net CO 2 flux varies by a factor of 2 to 3 for realistic values of the alkalinity-to-DIC ratio. Thus, alkalinity and carbon storage in ice has a significant impact on the net air-sea CO 2 flux in the model.

275
Next, we investigate the role of ice conditions, including the freezing-melting rate and ice concentration on the air-sea CO 2 flux. The NAPA model simulates a wide range of ice melt rates over the Arctic Ocean, spanning from 0 to over 7 m yr −1 and with areas of high ice melt in the Labrador and East Greenland Currents and the southern edge of the Beaufort Gyre (Fig. 3a).
The NAPA model also simulates freezing conditions that mostly occur when the lead fraction is close to 0 (Supplementary   Fig. 1 and 2. materials, Fig. S2). Indeed, over 88 % of the freezing days occur when the ice concentration is above 0.9. This supports the 280 assumption made in Sect. 2, where we considered freezing to mostly occur when the lead fraction is close to 0. Gross freezing rates and yearly integrated ice coverage are poorly correlated to ∆F m (r 2 =0.12 and r 2 =0.15 respectively).
Yearly net freezing-melting is more strongly correlated with ∆F m (r 2 =0.39). However, a better predictor of ∆F m is ice melt, excluding any freezing, hereafter called gross annual melt (r 2 =0.86; Fig. 4). An explanation for this strong relation is that winter ice cover prevents air-sea flux during the freezing period. This is an independent confirmation of the interpretation of 285 the mathematical derivation made in Sect. 2.
The high correlation between the gross annual ice melt (F M elt ) and ∆F m gives confidence in a linear model relating those two metrics: Another driver for ∆F m is the yearly integrated ice concentration (Fig. 4, colors), which is largest where full ice cover 290 persists for most of the year (Fig. 3b). While model experiments with lower ice coverage (dark blue) follow the regression well (solid black line), runs with higher ice coverage (light blue) have a steeper slope.
The 1D simulation ensemble can be used to calculate a yearly Arctic-wide increase due to ice-ocean carbon flux, for the 2014-2019 period. The ICE runs represent an increase of 30.0 ± 9.1 % (mean ± standard deviation calculated over the 6 years period), compared to the CTRL runs. Equation 6 depends on the parameterization of carbon ice-ocean flux and of air-sea CO2 295 flux, but is not otherwise model-specific. Therefore, it can be applied to other model outputs.

Application to an Earth System Model
The amplification of the air-sea CO 2 exchange due to the storage of carbon and alkalinity in ice is sensitive to the gross annual ice melt and the seasonality of the ice concentration. Both parameters are rapidly changing due to global warming. To investigate the impact of these changes on the supplementary carbon uptake, we turned to outputs from the ACCESS-ESM1.5 300 (Ziehn et al., 2020). This model, as any ESM, does not include any carbon storage in sea ice, although the freshwater flux between the ocean and sea ice is accounted for. We applied the linear relation in Eq. 6 to estimate the missing carbon uptake of CO 2 and, by adding it to the modelled carbon uptake, provide a corrected estimate of the oceanic carbon uptake in polar regions. While subduction processes are simulated in the initial outputs, our offline methodology does not correct mixing and advective carbon transport for the supplementary carbon due to the sea ice carbon pump. Therefore, an inherent assumption to 305 our methodology is the subduction of all the added carbon.
Although the linear relation between gross annual melt and ∆F m (Eq. 6) was derived from daily data, a very similar relationship is obtained when monthly data is used instead (RMSE between daily vs. monthly calculated gross annual ice melt <0.06 m yr −1 , not shown), giving us confidence that Eq. 6 can be used. The linear relation was applied to the extracted gross ice melt, resulting in a yearly supplementary carbon uptake for each grid cell. Spatially integrated over the area of interest, 310 this yields an Arctic-wide supplementary carbon uptake due to ice-ocean carbon flux, which can then be added to the modelderived carbon flux over the same area to yield a corrected carbon flux. The ratio between ∆F m and the model-derived carbon flux, expressed as a percentage, allows for easier interpretation of the magnitude of the process. This ratio can be interpreted as a measure of how much the ESM underestimates the Arctic Ocean carbon uptake. Those metrics were integrated over the different periods, yielding cumulative carbon uptake estimates over the historical and projection periods.

315
Due to the CO 2 undersaturation of the Arctic Ocean, the net carbon flux is positive (into the ocean) for all periods and scenarios. During the historical run, the modelled uptake slowly increases from 180 Tg C yr −1 in 1850 to 200 Tg C yr −1 in 1995 (a linear regression gives a slope of 0.26 Tg C yr −2 with r 2 = 0.5, p-value < 0.001), then stagnates during the last 20 years (r 2 = 0.0, p-value = 0.4) (Fig. 5a). The supplementary carbon flux, on the other hand, remains relatively constant over the whole period (Fig. 5a), meaning the corrected carbon uptake (Fig. 5a) follows a similar pattern as the model estimate. It 320 also leads to a slow decrease in the ratio of ∆F m over the model estimate (Fig. 5b). The increase in uptake may be driving increasing pCO 2 levels in the Arctic Ocean (Ouyang et al., 2020;DeGrandpre et al., 2020).
Projecting into the future, all three climate scenarios show a decrease in modelled and corrected carbon uptakes although interannual variability is high. In scenario SSP5-8.5 (Fig. 5a), and SSP1-2.6 to a lesser extent (Fig. 5a), carbon uptake increases until the 2040s, before dropping rapidly during the remainder of the century. The severe sea ice decline in SSP5-8.5 leads to a 325 similar decrease in ∆F m , while the two other scenarios show a relatively constant ∆F m over the 21st century.
Those scenarios differ in how large the fraction of ∆F m is compared to the total carbon uptake (Fig. 5b). Over the historical period, it slowly decreases starting above 15 % to arrive at around 12.5 % in 2015. It keeps decreasing in SSP5-8.5 to reach  Integrated over 1850-2100, the modelled carbon uptake sums up to 41.6, 40.2 and 42.3 Pg C for scenarios SSP1-2.6, SSP2-4.5 and SSP5-8.5 respectively, and the supplementary carbon uptake adds another 5.7, 5.6 and 5.3 Pg C respectively (Table 2).
Those cumulative supplementary carbon fluxes represent 12.5 to 14.1 % of the model-derived cumulative flux (Table 2).
Therefore, discarding the storage of carbon in sea ice in ESMs can lead to a significant underestimation of the carbon uptake in the Arctic Ocean, with varying impacts depending on the scenario considered, as described in the Discussion.

Discussion
In this study, the link between ice-ocean and air-sea carbon fluxes was investigated using two independent methods: a theoretical framework and numerical modelling. The methods provide consistent, complementary results, both pointing to a linear relationship between ∆F and ice melt and an exponential relation with the open-water fraction (Eq. 5 and Fig. 4).
Only three assumptions were made during the theoretical derivation. The assumption of a constant ∂pCO2 ∂DIC was addressed 340 in Sect. 2. The second assumption was a constant value of the mixed layer depth H 0 , also discussed in Sect. 2. The third assumption is the negligible effect of non-linearities in the carbonate system. Here, it is worth noting that the 1D numerical model does not rely on those assumptions and accounts for the varying ∂pCO2 ∂DIC and H 0 , and for the non-linearities of the carbonate system. Therefore, the good agreement between the theoretical framework and the model ensemble results builds confidence that these assumptions are justified. Back-of-the-envelope calculations using typical orders of magnitudes (pCO 2 345 = 350 µatm, changes in pCO 2 = 20 µatm) also show that non-linearities would represent less than 10 % of the total changes induced by the temperature, salinity, DIC and alkalinity variations, further supporting our assumptions.
A simplified version of the theoretical equation 5 can be evaluated with g and ∂pCO2 ∂t 1 H0 k g S CO2 considered as constant (cf. Supplementary Materials, Section S1.4), to better compare both methods. The ice concentration and freezing-melting flux used to force the 1D model (described in Section 3.1) can then be applied to this simplified version to calculate ∆F t (Fig. 6). While To interpret the relatively complex equation obtained in the theoretical framework (Eq. 5), we considered that ice formation is associated with ice-covered waters, related to the exponential term. Again, this simplification is supported by results from 355 the NAPA model (cf. Sect. 4.1 and Supplementary Fig. S2). The functional form may not apply to some regions with distinct ice regimes, including ice-exporting polynyas and ice-importing marginal ice zones. In those regimes, the exponential term and therefore the slope of the relation between ∆F and ice melt would be different. However, our solution is applicable to most of the Arctic Ocean.
We presented an approach for how Arctic carbon uptake estimates from ESMs can be corrected using our linear relation 360 between ∆F m and sea ice melt (Sect. 4.2). In doing so, past and potential future impacts of the sea-ice carbon pump in the Arctic can be analyzed. Our analysis suggests that uptake due to the sea-ice carbon pump increased during the historical period (Fig. 5a) due to longer open-water seasons and increased atmospheric pCO 2 . This is consistent with observations in the Canadian Arctic where higher pCO 2 levels are correlated with low ice extent (DeGrandpre et al., 2020). Because the sea ice carbon pump only applies to the seasonally ice-covered areas, the decline in ice extent translates into a stagnation of the 365 supplementary carbon uptake toward the end of the historical period and decreases during all SSP scenarios (Fig. 5a). In the SSP5-8.5 projection, the inhibition of the impact of carbon storage in sea ice is linked to drastic ice loss and therefore to less ice melt. In SSP1-2.6 and SSP2-4.5, the ice seasonal cycle remains significant, leading to a larger importance of ∆F m .
In all scenarios, except SSP5-8.5, we deem the current and future role of carbon storage and release by sea ice as nonnegligible. Without it, the ACCESS-ESM-1.5 model could be underestimating carbon uptake over seasonally ice-covered areas 370 by 5 to 15 %, or 10 to 15 % if we exclude SSP5-8.5. Note that this range differs from our calculation using the shorter NAPA model run from 2014 to 2019, in which the supplementary carbon uptake increases the yearly carbon uptake by 30.0 ± 9.1 % (mean ± standard deviation over 6 years). The discrepancy is mostly due to a lower ice melt simulated by the ACCESS model compared to the NAPA model (∼18 % lower), though both models have a reasonable agreement with satellite observations in terms of sea ice extent and concentration. We note that ACCESS-ESM-1.5 is the only CMIP6 model that provided the ice-ocean 375 freshwater flux and air-sea CO 2 flux, which are necessary inputs for our parameterization. An extension of this calculation to other ESMs would be possible if suitable output was available for more models.
Our estimated supplementary carbon flux is consistent with numbers given by Rysgaard et al. (2011) who suggested that the sea-ice carbon pump could represent 20 % of the air-sea CO 2 flux in open Arctic waters at high latitudes. Our estimates are higher than those from two other modelling studies. Grimm et al. (2016) reported that 7 % of simulated net polar oceanic CO 2 uptake is due to the sea ice carbon pump. Moreau et al. (2016) found a weakened Arctic carbon sink when including the sea-ice effect. While a direct comparison with those studies is difficult, we suggest that the vertical resolution is crucial for properly resolving the mechanisms at play. The coarse resolution used by Grimm et al. (2016) and Moreau et al. (2016) (9 and 10 layers in the first 100 m, respectively, compared to 9 layers in the first 10 m in our configuration) prevent them from capturing the shallow summer mixed layer observed in the Arctic. Using the same resolution as Moreau et al. (2016) in our 385 1D model leads to significant changes in the magnitude of the air-sea flux, either positive or negative depending on whether the mixed-layer depth is under-or over-estimated. The importance of high vertical resolution, capable of properly representing the shallow mixed layers in Arctic regions, is not surprising. On top of that, a proper representation of subduction, included in the Grimm et al. (2016) and Moreau et al. (2016) studies but beyond the scope of the present one-dimensional study, would be important to more fully understand the long-term fate of carbon in the global ocean. Yet, in an undersaturated ocean, the 390 amplification of the pCO 2 seasonal cycle can in itself explain an increased seasonal carbon uptake. Without any subduction, this would then lead the Arctic Ocean to reach equilibrium with the atmosphere faster than without accounting for the sea ice carbon pump, eventually saturating the surface ocean and reducing the carbon uptake. The output from the ACCESS-ESM1.5 model account for subduction, but the fate of supplementary carbon estimated here cannot be determined without a proper coupling of a sea ice biogeochemical component. It is therefore unknown whether carbon flux driven by advection and mixing 395 would proportionally increase and export the supplementary carbon or whether the latter would saturate the surface mixed layer, leading seawater pCO 2 to catch-up with atmospheric values faster than without accounting for the sea ice carbon pump.
While the amplified seasonal cycle of carbonate properties found in our study agrees well with Mortenson et al. (2020), they suggest a negligible impact of ice-ocean carbon flux on annual oceanic CO 2 uptake. A potential source of this disagreement could be their lower alkalinity-to-DIC ratio in sea ice (1.25 in their study, 1.8 for this study's reference case). We have shown 400 that the resulting supplementary carbon uptake is sensitive to this ratio (Sect. 4.1 and Fig. 2).
Our parameterization of the alkalinity-to-DIC ratio may be overly simplistic. First, the vertical profiles of alkalinity and DIC in sea ice, assumed homogeneous here, might be C-shaped to follow salinity profiles, though observations do not necessarily support a vertical heterogeneity (e.g. Miller et al., 2011;Rysgaard et al., 2009). As long as the parametrized values are representative of the freezing and melting ice over a seasonal cycle, we believe that the vertical homogeneity assumption is 405 reasonable. Second, the alkalinity-to-DIC ratio is known to increase over time. The ratio can change for several reasons: (1) CO 2 outgassing from ice to the atmosphere when brine is expelled at the surface  or when permeability is restored by increasing temperatures in early spring (Delille et al., 2014;Nomura et al., 2010) decreases DIC, although uncertainties in these fluxes are high (Watts et al., 2022), (2) primary production from ice algae consumes CO 2 and nitrate, therefore reducing DIC while increasing alkalinity (Delille et al., 2007;Rysgaard et al., 2007), (3) formation of ikaite crystals trapped in 410 sea ice retains alkalinity while CO 2 -enriched brine is exchanged with seawater (Rysgaard et al., 2007(Rysgaard et al., , 2009(Rysgaard et al., , 2011. However, the main driver of supplementary carbon uptake is sea ice melt, occurring towards the end of the seasonal cycle, when the alkalinity-to-DIC ratio is expected to be highest (Sect. 4.1). Therefore, applying a constant, high ratio is likely to best match real conditions while keeping the parameterization in its simplest possible form. Moreover, while the value of 1.8 might seem high, it is within the range of observed values (1 to 2, Miller et al., 2011;Rysgaard et al., 2009Rysgaard et al., , 2011. Nonetheless, a better constraint on this ratio is needed, which requires a proper understanding of the conditions of ikaite formation. The empirical linear relation determined in Sect. 4.1 (Eq. 6) involves annual ice melt only, to the exclusion of ice formation.
Outputs from the 3D numerical ice model show that whenever the freeze-melt rate is negative (i.e., ice is forming), the ice concentration is close to 1 preventing gas exchange. While this might be due to artifacts inherent to numerical models (e.g., lack of resolution of small leads), our linear relation is derived for application on the latter and therefore stands in this context.

420
It should be noted that we excluded shallow shelves from our runs, such as the Laptev Sea shelves. Those areas are highly productive with regard to ice formation in polynyas (exceeding 7 meters per year, Dmitrenko et al., 2009) and subject to active leads in winter. Therefore, in those regions, during ice formation, carbon storage in sea ice could yield anomalous outgassing, though intense ice formation has also been linked to enhanced CO 2 uptakes (Else et al., 2011). Brine sinking in those areas is also significant enough to form deep water masses and is therefore likely to provide a carbon export mechanism over multiyear 425 time scales. Investigating this mechanism would require a fully coupled 3D model.
In this study, a 1D model was used preferentially for computational reasons. This provided more flexibility for parameterization and sensitivity tests and allowed us to generate a large ensemble of simulations which would be computationally prohibitive with a full 3D model. For the same reason, we disabled the biological processes in our 1D model. It could be hypothesized that respiration will increase pCO 2 in winter when ice is acting as a lid and primary production will lower it in 430 summer, in phase with the chemical process described here, thus further amplifying the sea ice carbon pump. The storage and release of carbon by sea ice complete the picture drawn by the rectification hypothesis (Yager et al., 1995) which assumes that half of the air-sea CO 2 exchange that would be occurring in the typically ice-free ocean is cancelled by the presence of sea ice.
While this rectification hypothesis is fully applicable in areas of local ice formation and melt, the southern-most areas of our domain of interest (e.g., Labrador Current and East Greenland Current, Fig. 3) only involve melting of advected ice, usually in 435 winter and are therefore out of phase with the previously described seasonal cycle of pCO 2 . Melting of advected sea ice would then decrease pCO 2 and increase carbon uptake in winter without modifying it in summer. Deep convection events frequently happening in those areas could then have important consequences for the carbon export at depth, but this is beyond the scope of this study.

440
In this study, we used two independent but consistent approaches, a theoretical framework and numerical models, to explore the effect of storage and release of alkalinity and DIC by sea ice on air-sea CO 2 fluxes. Our theoretical derivation and numerical results show that the ice-ocean carbon flux amplifies the seasonal cycle of surface pCO 2 in phase with the seasonal cycle of sea ice concentration. This leads to a significant increase of oceanic carbon uptake in seasonally ice-covered areas in the Northern Hemisphere. One of the key findings of this study is that ice melt is a direct driver of the supplementary carbon 445 uptake and can therefore be used to correct carbon uptake estimates. This supplementary carbon uptake accounts for 30 % of Arctic Ocean carbon uptake according to our regional, high-resolution model and for 5 to 15 % in the global, lower-resolution ACCESS-ESM1.5 model, depending on the chosen scenario.
We also provide two novel relations to estimate the impact of sea ice carbonate on air-sea carbon flux. The first (cf. Eq. 5 for the full expression of ∆F t ), derived from a theoretical framework, can be useful for analyzing observational datasets 450 and decomposing sources of pCO 2 variability. The second, ∆F m = 113.6 · F M elt − 10.1, derived from a linear regression on numerical data, can be used to estimate the missing supplementary carbon uptake in numerical models that do not account for the sea ice carbon pump. An important strength of our theoretical framework is that no geographical assumption was made in its derivation. Eq. 5 can therefore be applied to both the Northern and Southern Hemispheres, keeping in mind that alkalinity and DIC values in sea ice may be different between both regions, due to environmental conditions (Delille et al., 2014;Fransson 455 et al., 2011;Rysgaard et al., 2011).
While the results presented here offer a straightforward way for estimating the missing carbon uptake in ESMs, additional sea ice and under-ice observations will help to better constrain the impact of carbon storage in sea ice onto air-sea fluxes. Furthermore, it seems prudent to add sea ice biogeochemistry to numerical models to reduce uncertainties, especially in regional studies.

460
This study emphasizes the importance of accounting for carbon storage in sea ice in numerical models for an accurate simulation of carbon fluxes in polar regions. Further model runs explicitly simulating the sea ice carbon pump in projection scenarios would help validate our results and would provide useful insights into the future carbon cycle in the Arctic and Southern Oceans, including the role of mixing and advective processes on the fate of the added carbon. A high vertical resolution would be crucial to properly resolve the shallow Arctic summer surface mixed layer and the carbon subduction.

465
Modelling studies dedicated to leads and polynyas would also help to qualify and quantify the sea ice carbon pump in those areas of intense mixing, as well as providing guidelines on how to parametrize those mesoscale ice features in low resolution ESMs. Observational constraints on the temporal and spatial variability of the alkalinity-to-DIC ratio in sea ice and a better mechanistic understanding of the fate of brine during ice formation season are crucial for properly simulating those processes.
The importance of the sea ice carbon pump should also be kept in mind when estimating fluxes from observations. A better 470 accounting of the sea ice carbon pump will also facilitate the global effort to better constrain the carbon cycle in the oceans and to understand its changes under climate change.
Data availability. The 1D model outputs are available at https://doi.org/10.5281/zenodo.7038942. The ACCESS-ESM1.5 data can be accessed at https://esgf-node.llnl.gov/search/cmip6/ (last access: 31 August 2022; Ziehn et al., 2020). The BGOS Mooring In Situ pCO2 and pH time-series are provided by Michael DeGrandpre and hosted by the Arctic Data Center, https://doi.org/10.18739/A28C9R46N (last access: Author contributions. BR, KF and EO conceived the study. BR carried out 1D model simulations, validation, and analyses. TB, XH and YL assisted with setup and validation of the NAPA model. TB assisted with ACCESS-ESM1.5 data access and analysis. BR, KF, EO and MDD discussed the results. BR, KF and EO wrote the paper with contributions from the coauthors. Competing interests. The authors declare that they have no conflict of interest. the work contains. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP. We thank the CSIRO for producing and making available the ACCESS-ESM1.5 model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP and ESGF. The BGOS mooring in-situ pCO2 time-series data used in this study for validation are available through the U.S. National Science Foundation (NSF) Arctic Data Center(https://arcticdata.io). BR would like to acknowledge NSF support for participating in the BGOS/JOIS