Domain Editor-in-Chief: Jody W. Deming; Department of Biological Oceanography, University of Washington, Seattle, Washington, United States
Associate Editor: Stephen F. Ackley; Department of Geological Sciences, University of Texas at San Antonio, San Antonio, Texas, United States

## 1. Introduction

Among the natural systems, the role of sea ice in CO2 cycling is not well constrained. Several studies report air-ice CO2 fluxes that show that sea ice is a permeable medium under certain conditions of temperature and salinity (Semiletov et al., 2004; Nomura et al., 2006, 2010b; Miller et al., 2011b; Delille et al., 2014), thereby refuting the assumption that sea ice impedes the air-ocean gas exchange. However, observational difficulties, in particular the lack of continuous observations covering the entire ice growth and decay cycle, hinder the understanding of carbon exchange processes in ice-covered seas. Some studies suggest that active sea ice processes are significant (Rysgaard et al., 2011; Delille et al., 2014), while others assume them to be negligible (Cross et al., 2014). Given the rapid sea ice changes in the Arctic Ocean in the context of global warming (IPCC, 2013), better constraints on the role of sea ice in CO2 cycling are needed to assess the future capacity of polar oceans to buffer the rise of atmospheric CO2 concentration.

Sea ice is a composite of pure ice, brine and gas inclusions. Various biological, chemical and physical parameters may affect the concentration of gases in sea ice (e.g., Tsurikov, 1979; Zhou et al., 2014a). Temperature appears to be one of the main controls of CO2 concentration in sea ice (Geilfus et al., 2012a; Delille et al., 2014). When temperature decreases, brine inclusions shrink, concentrating salts, gases and other impurities in the brine. In contrast, when temperature increases, the melting of ice with the increase in size of brine inclusions dilutes their content (e.g., Vancoppenolle et al., 2013; Delille et al., 2014). The increasing concentration of gases in brine due to the shrinking of brine inclusions triggers the nucleation of bubbles (Tsurikov, 1979; Killawee et al., 1998; Tison et al., 2002; Light et al., 2003). Gases are therefore found in sea ice in both dissolved (i.e., brines) and gaseous forms (i.e., bubbles). Bubbles are indeed observed at the onset of sea ice growth, when air inclusions are trapped within the ice structure. Crabeck et al. (2014b) and Zhou et al. (2014a) suggested that bubbles further develop during ice growth when the gases concentrates within the brine. Conversely, during ice decay, brine dilution and the related decrease of gas concentration is hypothesized to promote the dissolution of gases in the brines.

Sea ice is considered as permeable when the brine fraction is above 5% (Golden et al., 1998). As sea ice becomes permeable, air-ice gas exchange increases (Delille et al., 2014). Brine convection or gravity drainage (Notz and Worster, 2009; Hunke et al., 2011) is the main process responsible for the rejection of dissolved gases to the ocean during ice growth (Moreau et al., 2014). In the absence of convection, the diffusion of dissolved gases becomes the main pathway to transport dissolved gas across the ice (Gosink et al., 1976; Loose et al., 2011; Shaw et al., 2011). This diffusive flux is driven by dissolved gas gradients (Delille et al., 2007; Nomura et al., 2010b; Miller et al., 2011b). Gas bubbles may provide an alternative gas transport pathway (Zhou et al., 2013; Crabeck et al., 2014a). The efficiency of this pathway depends on whether bubbles can only diffuse molecularly (Loose et al., 2011) or rise due to their buoyancy (Moreau et al., 2014).

Considering the need to better understand carbon cycling in ice-covered seas, the absence of measurements over the whole ice growth and decay cycle and the difficulty of performing these measurements on natural sea ice, we carried out a controlled ice growth and decay experiment during which the air-ice CO2 fluxes were monitored using the chamber method. The CO2 transfer coefficient during both ice growth and ice decay was computed and compared to a sea ice model integrating carbon dynamics.

## 2. Methods

### 2.1 Experimental setting

The experiment was carried out at the Arctic Environmental Test Basin facility of the Hamburg Ship Model Basin (http://www.hsva.de) in the framework of the INTERICE V project. Eleven polyethylene bags of 1.2 m3 were filled with about 1000 L of filtered seawater from the North Sea. The experiment reproducing ice growth and ice decay took place over a period of 19 days. The day “0” of the experiment was 30 May 2012. The air temperature above the mesocosms was set to −14 °C the first 14 days (hereafter denoted as the “growth phase”), then to −1 °C until the end of the experiment (the “melting phase”). More information about the experimental settings are provided by Zhou et al. (2014b).

We carried out continuous in situ measurements of ice temperature and air-ice CO2 fluxes: A chain of 10 thermistors was placed at 2 cm intervals through the whole ice thickness. Ice cores were also collected regularly (every 1 to 3 days). On a given day, all samples (ice, under-ice water) came from the same mesocosm. Once the ice in a mesocosm was sampled, it was compromised and not used again in the experiment. So each sampling day corresponded to a different mesocosm, except on day 19, when we sampled two mesocosms. Ice cores were wrapped in polyethylene bags for storage below −25 °C in the dark and for subsequent measurements of bulk ice salinity, total alkalinity (TA) and partial pressure of CO2 (pCO2).

### 2.2 Ice pCO2 at high vertical resolution

We applied the method developed by Geilfus et al. (2012b) and reviewed by Crabeck et al. (2014b) and Geilfus et al. (2015) for the measurement of the bulk pCO2 (denoted as pCO2 bulk) within permeable sea ice. The goal is to equilibrate the sea ice samples with a N2/CO2 gas mixture of known concentration at a temperature as close as possible to the in situ temperature (the temperature of the ice upon sampling). Samples were cut to fit tightly into a square container, 4 x 4 cm, that was 4.4 cm high to minimize the headspace. The container containing the sample was sealed and connected to a vacuum pump for 15 min. A standard gas of known concentration (500 ppm of CO2) was then injected at a pressure of 1013 mbar. Standard gas and ice sample were equilibrated for 20 hours in a constant temperature bath at the in situ temperature. Gas was then recovered and injected in a Varian 3300 gas chromatograph to measure the CO2 concentration. Shortly afterward, the sample temperature was measured to check for experimental drift. Measured pCO2 bulk was corrected for the temperature difference between the sample and the in situ temperature according to the corrections proposed by Copin-Montegut (1988).

### 2.3 Total alkalinity

We measured total alkalinity (TA) on melted sea ice (TAbulk) using 50 g of sea ice melt collected at a 2 cm vertical resolution. TA was measured on all mesocosms except SW 9. We derived TA for all SW 9 ice sections using the strong linear regression between salinity and TA observed for all of the other samples. We also collected seawater for TA measurement. Melted bulk ice and seawater samples were poisoned with a solution of super-saturated HgCl2 and then stored in the dark, until analysis (one year after the sampling). TA was measured by open-cell titration with 0.11 M HCl and the endpoints were determined according to Gran (1952). Routine analyses of Certified Reference Materials (provided by A. G. Dickson, Scripps Institution of Oceanography) ensured that the uncertainty of the TA measurements was less than 4 μmol kg−1.

### 2.4 Dissolved inorganic carbon

We computed the dissolved inorganic carbon of bulk ice from TAbulk and pCO2 bulk using a 2 step computation. We estimated the salinity of brines according to the ice temperature using the relationship of Cox and Weeks (1983) derived from data compiled by Assur (1960). We then computed TAbrines at the salinity of brines from the linear relationship between TA and salinity. We computed the total dissolved inorganic carbon (DIC) of brines (denoted as DICbrines) from TAbrines and pCO2 bulk using the CO2SYS program for the carbonate system (Lewis and Wallace, 1998). We used the CO2 dissociation constants of Mehrbach (1973) refitted by Dickson and Millero (1987) and the other constants advocated by Dickson and Goyet (1994). We then converted DICbrines to DICbulk assuming a linear relationship between DIC and salinity. Nonetheless, there are some limitations with this approach that should be noted: the dissociation constants have been established for the ranges of temperatures and salinities of open ocean waters (i.e., temperatures above 1 °C and salinities of 35). We assumed that the CO2 dissociation constants were applicable at sub-zero temperatures, as suggested by Marion (2001) and Delille et al. (2007). We refer the reader to Brown et al. (2014) for a discussion on the validity of the constants.

### 2.5 Seawater pCO2

To measure the underlying seawater pCO2 a hole was drilled through the sea ice cover. Seawater was pumped from the hole using a peristaltic pump (Masterflex® - Environmental Sampler) and supplied to a sea ice equilibrator system (SIES; Delille et al., 2007) for measurements of the pCO2 and recycled back to the seawater through the same hole. The SIES is based on a membrane contractor equilibrator (Membrana® Liqui-cell) coupled to an infrared gas analyzer (IRGA, Li-Cor® 6262). Seawater flowed into the equilibrator at a maximum rate of 1 L min−1 and a closed air loop ensured circulation through the equilibrator and the IRGA at a rate of 3 L min−1. The IRGA was calibrated before and after the experiment with N2 and CO2:N2 mixtures with mixing ratios of 388 and 813 ppm supplied by Air Liquide Belgium. During the experiment, the drift of the IRGA was corrected with N2. Uncertainty during this experiment was less than 6 µatm.

### 2.6 Air–ice CO2 fluxes

#### 2.6.1 Measurements at the automated chamber

In this paper, positive CO2 flux refers to CO2 flux from the ice to the atmosphere, while negative CO2 flux refers to a flux from the atmosphere to the ice. We measured air-ice CO2 fluxes using an automated chamber placed above the water surface or on top of the ice. The chamber consisted of a mobile cap and a plastic cylinder, or so-called collar, with a diameter of 20 cm and a height of 9.7 cm. A rubber seal surrounded the cylinder and ensured an airtight connection between the ice and the chamber. Each hour, the cap closed the chamber and the pCO2 was measured over 15 min. At the beginning of the experiment one chamber was set above the surface of the water with the collar lowered a few millimetres below the water surface of a dedicated mesocosm apart from the 21 mesocosms used for ice collection. However, ice freezing and consolidation pushed the collar upward so that the collar was not properly sealed. After the fifth day of the experiment, the chamber was therefore moved to mesocosm 11 and was properly sealed. A pump within the LI-COR Multiplexer (LI-8150) circulated the air in the chamber at a flow rate of 2.1 L min−1. When the pCO2 of ice is higher than atmospheric pCO2, CO2 is transferred from the ice to the atmosphere and the automated chamber records a positive flux. A negative flux is observed in the opposite case. Water-corrected CO2 flux was computed automatically with LI-8100 File Viewer 3.1.0 package provided by LI-COR Biosciences. The flux was either calculated with a linear or an empirical exponential regression depending on which method provided the best fit (assessed from the normalized sums of the squares of the residuals).

#### 2.6.2 Computation of a gas transfer coefficient for CO2

As described above, gases are transported through sea ice to the atmosphere by convection, diffusion and/or the ascent of bubbles to the ice surface. In the present study, we calculated an effective gas transfer coefficient for CO2 (K, in mmol m−2 d−1 atm−1), using the equation developed by Liss and Slater (1974) and Sarmiento and Gruber (2004):

$F=K\left({p\mathrm{CO}}_{2\mathrm{bulk}}–{p\mathrm{CO}}_{2\mathrm{air}}\right)$
(1)

where F is the air-ice CO2 flux in mmol m−2 d−1, pCO2 bulk is the pCO2 in the ice and pCO2 air is the pCO2 in the air, both expressed in atm. In this equation we assume that F and K reflect both diffusive flux and bubble buoyancy (i.e., the rise of bubbles to the ice surface).

### 2.7 Assessment of the precision of derived variables

A Bootstrap resampling statistical analysis procedure, using random values of the measured parameters (temperature, salinity, pCO2, DIC and TA) between the mean ± precision over 1000 iterations, was used to estimate the propagation of errors to the computed parameters (ΔDIC and K). This method was used as a way to show the effects of the imprecision of the data set on the calculated parameter.

### 2.8 Modelling air-ice CO2 fluxes

In order to understand the factors that drive air-ice CO2 fluxes during the experiment, we ran a one-dimensional thermodynamic sea ice model representing sea ice halo-thermodynamics and carbon dynamics, including air-ice CO2 fluxes (Vancoppenolle et al., 2010; Moreau et al., 2015). Vertical carbon transport and air-ice carbon fluxes are explicitly separated between dissolved and gaseous form contributions.

The dissolved gas pathway combines the vertical transport of dissolved inorganic carbon in brine and the diffusive air-ice CO2 flux, FCO2 (mmol m−2 d−1), which is assumed proportional to the CO2 partial pressure (pCO2) difference between the surface brine (in the top 5 cm of the ice) and the atmosphere, and is a function of the near-surface brine fraction:

${F}^{{\mathrm{CO}}_{2}}={k}^{{\mathrm{CO}}_{2}}·{e}^{2/3}·\left({\zeta }^{{\mathrm{CO}}_{2}}–{K}_{0}{f}^{{\mathrm{CO}}_{2}}\right)$
(2)

where kCO2 (m s−1) is the piston velocity, e is the brine fraction near the ice surface, e2/3 represents the fractional surface open to brine-air diffusive CO2 fluxes, ζCO2 is the brine CO2 concentration (mmol m−3), f CO2 is the atmospheric CO2 fugacity (atm), and K0 is the Henry’s constant (mmol m−3 atm−1). We consider that fCO2 = patm * rCO2,atm, where patm is the atmospheric pressure (atm) and rCO2,atm is the atmospheric CO2 mixing ratio measured within the chamber. The piston velocity is calculated from the molecular diffusion coefficient of dissolved CO2 (Ddiff in m2 s−1) and the thickness of the diffusive boundary layer (zBL in µm):

${k}^{{\mathrm{CO}}_{2}}={D}_{\mathrm{diff}/{Z}_{\mathrm{BL}}}$
(3)

This piston velocity neglects the effects of snow and wind, absent in these tank experiments, and only includes the contribution of dissolved CO2. In contrast, equation 1 includes both gas bubble and dissolved contributions as F and pCO2 bulk are direct measurements without distinction between dissolved and gaseous forms. The zBL is highly uncertain and therefore has been used as a tuning parameter to adjust the magnitude of air-ice CO2 fluxes, whereas Ddiff is better constrained by observations. For our control simulation (CTRL), we used Ddiff = 0.97 10−9 m2 s−1 (the diffusion coefficient of CO2 in water from Broecker and Peng, 1974) and zBL = 0.5 µm (Moreau et al., 2015; Table 1).

10.12952/journal.elementa.000112.t001

Table 1.

#### Description of the sensitivity runs used to test the sensitivity of air-ice CO2 fluxes to model parameterizationa

Run Name Rbub eTgas Ddiff zBL Ikaite precipitation Number of ice layers Misfit
1 CTRL 10 0.07 0.97 0.5 Yes 10 0.1
2 Loose-Ddiffb 0 0.07 24 0.5 Yes 10 0.08
3 No-bub 0 0.07 0.97 0.5 Yes 10 0.08
4 Low-zBL 0 0.07 0.97 0.05 Yes 10 0.08
5 No-ikaite 0 0.07 0.97 0.5 No 10 0.10
6 20-layers 0 0.07 0.97 0.5 Yes 20 0.06
7 Low-bub 1 0.07 0.97 0.5 Yes 10 0.07
8 High-bub 20 0.07 0.97 0.5 Yes 10 −0.05
9 Moreau-bubc 0.1 0.07 0.97 0.5 Yes 10 0.08

The gas bubble pathway was developed in the model to simulate Argon dynamics (Moreau et al., 2014) and implies explicit gas bubble reservoirs in every layer. The gas concentration in the bubble reservoir changes due to bubble nucleation/dissolution, upward migration of buoyant gas bubbles, and bubble escape to the atmosphere. Nucleation of gas bubbles transfers dissolved CO2 from brine to the bubble compartment as a function of the CO2 super-saturation. At each time step, a fraction Rbub of the CO2 super-saturation is transferred to bubbles. The bubbles migrate upward when the brine network is connected, which is assumed to happen above a given brine fraction threshold (e.g., eTgas = 0.07, Zhou et al., 2013; Moreau et al., 2014). If the fraction of sea ice with e > eTgas includes the ice surface, all gas bubbles escape to the atmosphere and contribute to the air-ice CO2 flux. For our CTRL simulation, we used Rbub = 10% h−1 and eTgas = 0.07 (Table 1).

In order to understand the potential reasons to explain the underestimation of the observed air-ice CO2 fluxes by the model, we performed three series of sensitivity experiments (Table 1). First, we tested the impact of a more intense dissolved CO2 pathway, by changing Ddiff and zBL (runs 2–4). The Ddiff value of Broecker and Peng (1974) is derived from seawater and only includes diffusive effects. We tested the slightly higher Ddiff value of Loose et al. (2011), derived from sea ice experiments, which can potentially include both diffuse and bubble contributions. A lower zBL (0.05 µm) value was also tested. In a second series of experiments (runs 5–6), the impacts of higher resolution (20 layers instead of 10) and ikaite precipitation/dissolution were investigated. Third, we tested the impact of a more intense gas bubble pathway by changing Rbub, the bubble nucleation rate (runs 7–9), which is highly uncertain.

To evaluate which scenario is closest to the observed fluxes, we computed a model-data misfit for each of the model runs. The misfit were computed from the difference between the flux observed and the flux modelled. We report averages for each of the model runs in Table 1.

All other parameters are from Moreau et al. (2015). The simulation spans 31 May – 18 June 2012. The time step is 1 h. Based on the present observations, initial seawater TA and DIC concentrations were set to 2300 mmol m−3 and 2090 mmol m−3, and sea ice TA and DIC concentrations were set to 850 mmol m−3 and 750 mmol m−3.

## 3. Results and discussion

Detailed information about sea ice physical properties can be found in Zhou et al. (2014b). Briefly, the first phase of the experiment – the growth phase – lasted from day 1 to day 15. Air temperature above the tank was set to −15 °C and sea ice grew continuously, reaching a maximum thickness of 24 cm. A strong temperature gradient was observed between the top and the bottom of the ice. Salinity exhibited a typical C-shape profile with a lower salinity at the ice interior compared to the top and the bottom of the ice. Then from day 16 to 19, the air temperature was set to −1 °C, sea ice thickness slightly decreased, temperature exhibited a more homogeneous profile through the whole thickness, and salinity decreased in the top and bottom parts. The brine volume fraction remained above 5% during the whole experiment. On the whole, the sea ice remained thin, warm and permeable during both the growth and the melting phase.

### 3.1 Total alkalinity

TA concentrations in bulk sea ice (i.e., melted sea ice) are consistent with the values reported in the literature (Gleitz et al., 1995; Delille et al., 2007; Nomura et al., 2010b; Fransson et al., 2011; Miller et al., 2011a; Geilfus et al., 2012a; Rysgaard et al., 2013) and were highly correlated with salinity (r2 = 0.95; Figure 1).

doi: 10.12952/journal.elementa.000112.f001.
Figure 1.

#### Linear regression between TA and salinity within bulk sea ice.

The ice growth and melting phases are colored in blue and red, respectively.

The repeated inspection of freshly melted ice using a binocular microscope (Leitz Laborlux® with 125 to 500 x magnification) did not reveal the presence of ikaite, a hydrated calcium carbonate polymorph (CaCO3·6H2O). Taking into account the thermodynamic constraints (Papadimitriou et al., 2013) and the kinetics to simulate the ikaite precipitation (Papadimitriou et al., 2014), the maximum value of ikaite simulated by the model in the ice surface layer was 13 µmol kg−1, whereas the error on TA measurement was 4 µmol kg−1.(parameters from the CTRL simulation of Moreau et al., 2015). This concentration falls at the lower end of the range of ikaite concentrations reported in sea ice (7–93 µmol kg−1 in Dieckmann et al., 2008; 15–19 µmol kg−1 in Geilfus et al., 2013; 100–900 µmol kg−1 in Rysgaard et al., 2013). This evidence indicates that if ikaite precipitated during the experiment, it was not significant.

### 3.2 CO2 exchange at the air-ice interface

#### 3.2.1 Continuous measurements of air-ice CO2 fluxes

During seawater cooling, and before the formation of the ice crystals, CO2 fluxes measured with the automated chamber showed negative values, down to −6 mmol m−2 d−1 (Figure 2). Within hours after the formation of the first ice crystals, CO2 fluxes became mostly positive (ranging between −0.4 mmol m−2 d−1 and 0.75 mmol m−2 d−1 with an average of 0.2 mmol m−2 d−1 for the growth phase), consistent with the observed super-saturation of CO2 in bulk ice (i.e., pCO2 bulk above 400 ppm). During the melting phase, CO2 fluxes turned to negative (ranging between −2.1 mmol m−2 d−1 and 0 mmol m−2 d−1 with an average of −0.24 mmol m−2 d−1) in parallel with the decrease of pCO2 bulk that passed below saturation. The pCO2 of the underlying seawater remained under-saturated during the whole experiment, while the surface (first 5 cm) pCO2 bulk showed values above or below the atmospheric pCO2 during the growth and the melting phase, respectively.

doi: 10.12952/journal.elementa.000112.f002.
Figure 2.

#### Air-ice CO2 fluxes, sea ice surface temperature and pCO2 over time.

a) Evolution of the air-ice CO2 flux (black corresponds to a period with the chamber not properly sealed, blue to a period with a correct sealing of the chamber) and temperature 2 cm above the air-ice interface (in red). The green horizontal step line corresponds to the flux calculated from the DIC anomaly. The shaded areas refer to working hours and the related increase of atmospheric pCO2. b) Atmospheric (black line), ice surface (green triangles down) and seawater (blue triangles up) pCO2. The ice surface pCO2 corresponds to the first 5 cm of the bulk sea ice measured with the method of high vertical resolution.

Although the chamber was not properly sealed between day 0 and day 5, air-ice CO2 fluxes from the whole measurement period were consistent with previous measurements carried out with chambers over artificial sea ice (between 0 and 0.27 mmol m−2 d−1; Nomura et al., 2006) and slightly lower than measurements over natural sea ice. Delille et al. (2014) measured CO2 fluxes ranging from −5.2 mmol m−2 d−1 to 1.9 mmol m−2 d−1 on Antarctic pack ice in spring, and ascribed these fluxes to seasonal pCO2 gradients between the brine and the atmosphere. Geilfus et al. (2012a) measured CO2 fluxes at the sea ice interface ranging from −2.63 mmol m−2 d−1 up to 0.84 mmol m−2 d−1 in the Arctic coastal zone, while Nomura et al. (2010a) measured CO2 fluxes ranging from -1 mmol m−2 d−1 to 0.7 mmol m−2 d−1 over land fast ice in Barrow at the end of spring.

The temperature and the pCO2air measured within the automated chamber exhibited daily variations, rising during “working hours” due to the presence of the researchers in the experimental room and decreasing outside “working hours”. The pattern of air-ice CO2 fluxes is clearly opposite and strongly affected by the rise of atmospheric pCO2 during working hours. These imperfections are inherent to the environmental constraints of this experimental study and could hardly be avoided. Nevertheless, the results can still be interpreted in terms of the main factors driving the CO2 fluxes, and the overall pattern of observed air-ice CO2 fluxes is consistent with in situ observation of a seasonal cycle, with upward CO2 fluxes during ice growth (Nomura et al., 2010a; Geilfus et al., 2012a; Delille et al., 2014) and downward CO2 fluxes during ice melt (Delille et al., 2014; Geilfus et al., 2015).

#### 3.2.2 Integrated estimates of air-ice CO2 fluxes

In parallel with direct measurements of air-ice CO2 fluxes, we derived air-ice CO2 fluxes from DIC depletion in the top layers of the ice (Geilfus et al., 2013). DIC was normalized to a salinity of 7 (7 corresponds to the mean bulk ice salinity during the whole experiment, noted as DIC7) to remove the salinity-related changes (brine rejection, concentration and dilution; Figure 3). If no biogeochemical processes occurred (biological activity, CaCO3 precipitation, CO2 transfer to the gas phase and CO2 exchange with the atmosphere), DIC7 profiles should be homogeneous over the ice column. During the growth phase, a clear decrease in the top 12.5 cm of the young ice was measured compared to the bottom ice horizons. Because no significant primary production was possible due to absence of primary producers (Zhou et al., 2014b) and ikaite precipitation was insignificant (section 3.1), we assumed that DIC7 depletion in the top 12.5 cm of the growing ice was due to CO2 release from the super-saturated ice to the atmosphere above. Whereas, during the melting phase we observed a DIC7 increase in the top 12.5 cm. This increase can be linked to the downward air-ice CO2 flux measured during the melting phase.

doi: 10.12952/journal.elementa.000112.f003.
Figure 3.

#### Bulk sea ice DIC and surface DIC depletion during ice growth.

a) Bulk DIC profiles and b) DIC profiles normalized to a salinity of 7 (DIC7), where the grey area refers to the DIC depletion zone.

Besides the measurements of air-ice CO2 fluxes, the amount of CO2 released to the atmosphere during the sea ice growth was assessed using a method proposed by Geilfus et al. (2013). First, a “theoretical DIC”(DICth) was calculated from the raw DIC concentration at 12.5 cm (DIC12.5 cm), assuming that, if biogeochemical processes are null, DIC and salinity (S) should follow a linear relationship:

${\mathrm{DIC}}_{\mathrm{th},i}={\mathrm{DIC}}_{12.5\mathrm{cm}}·\left({S}_{i}/{S}_{12.5\mathrm{cm}}\right)$
(4)

where DICth, i and DIC12.5 cm are computed for each sampling event and expressed in µmol kg−1. The index, i, refers to the upper two sampling depths, 2.5 cm and 7.5 cm. From the theoretical DIC, we derive the “DIC anomaly” for each sampling day (ΔDICt in mmol m−2):

${\mathrm{\Delta DIC}}_{t}=\Sigma \left(i=2.5\mathrm{to}12.5\mathrm{cm}\right)\left({\mathrm{DIC}}_{\mathrm{th},i}–{\mathrm{DIC}}_{i}\right)$
(5)

As the DIC anomaly was mainly due to the release of CO2 from the ice to the atmosphere during the present experiment (see above), integrating DIC anomalies over time gives the theoretical air-ice CO2 flux between two sampling events:

$F=\left(\left({\mathrm{\Delta DIC}}_{t+\mathrm{\Delta t}}–{\mathrm{\Delta DIC}}_{t}\right)/\mathrm{\Delta t}\right)·\mathrm{dx}·\rho$
(6)

where Δt is the time frame between two sampling events in days, dx is the distance in meters between two sampling depths, and ρ is the sea ice density for each layer defined from temperature and salinity and using relationships given by Cox and Weeks (1982), in kg m−3. Based on the precision of the measured variables (salinity, temperature, TA, pCO2; see Methods section), the precision on F is assessed at ± 0.03 mmol m−2 d−1 (which represents an error of 15% for the averaged CO2 flux of 0.2 mmol m−2 d−1).

From these DIC anomalies, the calculated air-ice CO2 fluxes (black triangles in Figure 2a) are in good agreement with the observed air-ice CO2 fluxes, except on day 16 where the difference represents more than 1 mmol m−2 d−1; integrated measurements cannot accurately capture the rapid transition from the freezing to the melting decay phase. With the exception of this transition phase, the consistency between the continuous and integrated fluxes suggest that in an indoor experiment, the automated chamber provides accurate measurements of air-ice CO2 fluxes when snow cover is absent, despite the large variation of the atmospheric pCO2.

### 3.3 Determination of a gas transfer coefficient for CO2 in artificial sea ice

Using equation 1 for gas exchange, we calculated a gas transfer coefficient for CO2 at the air-ice interface for each sampling day of the experiment (see K; Table 2), where F corresponds to the air-ice CO2 fluxes measured with the automated chamber, pCO2 air to the pCO2 at the air interface and pCO2 bulk to the pCO2 at the ice interface (0 cm to 5 cm). The numbers in light grey in Table 2 refer to the transition phase which was not taken into account in the final value computed for K. The values for K and ancillary data detailed in Table 2 are close during the whole growth period. The mean K for the sea ice growth period, from day 2 to day 14 inclusive, is K = 2.5 mol m−2 d−1 atm−1. A Bootstrap resampling statistical analysis of the propagation of the uncertainties of the measured variables (salinity, temperature, measured CO2 fluxes, pCO2 bulk and pCO2 air) gives an uncertainty for K of 0.9 mol m−2 d−1 atm−1 where standard error deviation is only 0.25 mol m−2 d−1 atm−1. For the melting phase (which corresponds to two measurements on day 19) the mean K value was 0.4 mol m−2 d−1 atm−1 (the uncertainty derived from a Bootstrap resampling statistical analysis was 0.08 mol m−2 d−1 atm−1). K was therefore 6 times higher during the growth phase than during the melting phase.

10.12952/journal.elementa.000112.t002

Table 2.

#### Gas transfer coefficient calculated for each day of the experimenta

Day of the experiment T Bulk S K
2 −3.0 11.5 2.43 (± 0.9)
5 −5.1 8.8 2.59 (± 0.9)
8 −5.0 9.5 2.66 (± 0.9)
14 −5.8 8.0 2.22 (± 0.9)
15b −3.3 3.0 2.1
16b −2.5 6.1 2.9
19 −1.8 6.9 0.5 (± 0.08)
19 −1.8 6.9 0.3 (± 0.08)

In order to get a measure of gas exchange that does not depend on the gas solubility, Wanninkhof (1992) proposed to use k, the gas transfer velocity, defined from:

$K=k.\mathrm{Sol}$
(7)

where Sol (mmol m−3 atm−1) is the solubility of CO2 in salt water (here brine), a function of temperature (T) and salinity (S), following Weiss (1974):

(8)

We used the T and S of sea ice brine. The resulting value, k = 0.164 (± 10%) cm h−1, falls within the range of values given by Liss and Merlivat (1986) for a smooth surface regime and the values given by Crusius and Wanninkhof (2003) for gas exchange over a lake at low wind speed. However, solubility as described by eq. 8 is most likely not adequate for the range of temperature and salinity encountered within sea ice, as suggested by Zhou et al. (2014a). K, the gas transfer coefficient, is solely deduced from measurements and does not depend on an uncertain solubility value.

Thin sections of the ice cores were processed following Tison et al. (2002). Bubbles were observed during the growth phase, particularly near the ice surface (Figure 4). We suggest that these bubbles (trapped at the very beginning of ice growth or newly formed) were likely moving upward due to their own buoyancy. Once bubbles reached the ice surface, they collapsed and CO2 was released to the atmosphere, corresponding to an ebullition flux. K-values calculated during the growth and melting phases indicate that gas transport was about 6 times faster during ice growth than during ice decay. We hypothesize that air-ice CO2 fluxes during the melting phase are driven mainly by molecular diffusion, while during the growth phase, ebullition fluxes add to molecular diffusion. This hypothesis suggests that the ebullition fluxes are a dominant pathway for gas transport and exchange during the growth phase.

doi: 10.12952/journal.elementa.000112.f004.
Figure 4.

#### Thin sections of surface ice cores.

Thin sections of ice cores for day 2, 5, 14 and 15 of the experiment. The contrasted circles and tubes are bubbles while the gray-shaded features are brine inclusions.

### 3.4 Model sensitivity experiments on the CO2 transport pathways through sea ice

Here we use the model to frame the CO2 flux into further theoretical considerations. The various model simulations described in Table 1 are compared to the observations, in order to understand which parameters are most influential in shaping the model response (Figure 5).

First, the sign of the CO2 flux is consistent between observations and all model simulations. All model runs show an efflux of CO2 during the growth period and an influx during the melt period. The direction of the flux depends mostly on the dissolved gas pathway, itself driven by the air-brine pCO2 difference, ultimately set by the near-surface ice temperature evolution, which is quite well understood. In addition, some of the short-term variations of the observed CO2 fluxes are reproduced. The model gives a good estimation of downward CO2 fluxes during ice melting although the time of the onset of downward CO2 fluxes occurs earlier than observed in the model.

doi: 10.12952/journal.elementa.000112.f005.
Figure 5.

#### Observed and simulated air-ice CO2 fluxes.

a) Simulated air-ice CO2 fluxes (mmol m−2 d−1) for runs Loose-Ddiff (grey, where Ddiff = 24 10−9 m2 s−1), No-bub (green, where Rbub = 0% h−1 and Ddiff = 0.97 10−9 m2 s−1) and Low-zBL (dark blue in dash, where zBL = 0.05 µm). The three curves for runs Loose-Ddiff, No-bub and Low-zBL are confounded on the plot. Observed CO2 fluxes are given in blue.

b) Simulated air-ice CO2 fluxes (mmol m−2 d−1) for runs No-ikaite (dark red, where there is no ikaite precipitation) and 20-layers (dark blue, where the number layers of ice was set to 20 instead of 10). Observed CO2 fluxes are given in blue.

c) Simulated air-ice CO2 fluxes (mmol m−2 d−1) for runs CTRL (in black), Low-bub (grey, where Rbub = 1% h−1), High-bub (cyan, where Rbub = 20% h−1) and Moreau-bub (dark grey, Rbub = 0.1% h−1). Observed CO2 fluxes (mmol m−2 d−1) are given in blue.

By contrast, the magnitude of the air-ice CO2 flux is highly variable among the different model runs, in particular during the growth phase. The most influential but uncertain factors clearly lie within the gas bubble pathway (Figure 5c). By contrast, the dissolved pathway is found rather insensitive to parameterization choices (e.g., gas diffusivity Ddiff, boundary layer thickness zBL, vertical resolution, ikaite, etc.). In turn, the total (dissolved + bubble) air-ice CO2 flux is essentially determined by the gas bubble nucleation rate Rbub, which acts as a bottleneck parameter for the efflux of CO2 to the atmosphere via the gas bubble pathway. This bottleneck effect explains why Rbub was used as a tuning parameter, set to 10% h−1 in the CTRL run, the value that gave the best agreement with observations.

Why is the dissolved pathway rather insensitive between the different runs? As in Moreau et al. (2015), we investigated whether the dissolved gas pathway could have been more intense than simulated. To answer this question, we tested to what extent modifying the uncertain parameters zBL and Ddiff could increase the upward CO2 fluxes during ice growth. None of our trials was successful. For example, increasing Ddiff up to 2.4 x 10−8 m2 s−1 (run 2: Loose-Ddiff) barely affects the upward CO2 fluxes during ice growth. Decreasing zBL down to 0.05 µm did not significantly increase the simulated air-ice CO2 fluxes either. This model behavior arises because the ultimate limitation for the air-ice CO2 diffusive flux is the near surface stock of DIC: supply to the near-surface DIC is not rapid enough to sustain CO2 fluxes to the observed magnitude. This result was already obtained by Moreau et al. (2015; see their Figure 9b where the simulated outward CO2 fluxes during winter at Barrow reach asymptotically 1.3 mmol m−2 d−1 when zBL is decreased down to 10−3 µm). By contrast, the gas bubble pathway is much more variable between the different runs and characterized by much larger uncertainties, because gas bubble flux can originate further down in the ice and is not limited by near-surface DIC stocks. The model-data misfit presented in Table 1 shows that the best model-observation agreement was achieved for a strong gas bubble pathway, which corresponds to high values of Rbub (10-20% h−1). Yet, the bubble nucleation rate Rbub is far from being well constrained from observations. In Moreau et al. (2014), this value was tuned to 0.1% h−1, corresponding to a characteristic time scale of 40 days, in order to match observed gas bubble concentrations over simulations spanning several months. In the present context (curve Moreau-bub in Figure 5c), such value gives a weak gas bubble pathway and an underestimation of CO2 fluxes by the model. Forty days is rather long for bubble nucleation, much larger than expected, typically less than one hour (see, e.g., Brennen, 1995). The most obvious reason why such an unrealistically large 40 day value had to be used by Moreau et al. (2014) is that pressure effects are missing in the model, which should quickly inhibit gas bubble nucleation once a large enough stock of gas bubbles is built up, which would typically occur for impervious sea ice over a few days. Larger Rbub values of 10–20% h−1 strongly intensify the gas bubble pathway. They also correspond to more realistic time scales for gas bubble nucleation (4–8 hours).

The influence of several other processes, known to occur in natural conditions and contribute to the ice-air CO2 fluxes, was also considered. As far as frost flowers are concerned, they were not observed in the tank experiment, and hence are not relevant to the present discussion of air-ice CO2 fluxes. The precipitation of ikaite could have influenced ice-atmosphere CO2 fluxes (e.g., Geilfus et al., 2013; Rysgaard et al., 2014); however, as discussed above, ikaite precipitation was not observed under the microscope and the model CTRL simulation predicts a rather low ikaite concentration of 13 µmol kg−1. The lack of ikaite influence is further corroborated by the model run with no ikaite precipitation (and no bubbles; Run 6, Table 1) which barely differs from the no-BUB run in terms of air-ice CO2 flux (Figure 5b).

We also tested the role of vertical resolution of the simulated CO2 fluxes by running the model with 20 vertical ice layers instead of 10, and with no bubbles (Figure 5b). More layers better resolve and increase the near-surface pCO2, which slightly intensifies the dissolved gas pathway during the growth phase.

### 3.5 Synthesis

As detailed in Zhou et al. (2014b), the artificial sea ice was likely permeable to gas exchange during the whole experiment. Convection was likely only present during the growth period and limited to near the ice-water interface (see Rayleigh numbers in Figure 3 of Zhou et al., 2014b) and was nearly negligible during the whole period of decay. We therefore consider that only diffusion and/or buoyancy processes occurred in the upper parts of the ice. Both the computation of the gas transfer coefficient and the model sensitivity analysis support this assumption. As described previously, bubbles were observed during the growth phase, particularly near the ice surface (Figure 4).

The values of K were about 6 times lower during the melting phase as compared with the growth phase (Table 2). At this stage, based on model simulations, we cannot explain this difference without assuming a significant contribution of the gas bubble pathway during the growth of sea ice. However, this explanation is not based on direct observation of gas bubble processes but rather on inference from model simulations characterized by large uncertainties, in particular in terms of the gas bubble nucleation rate. The latter is the key tuning bottleneck parameter sourcing gas bubbles to the model gas bubble ice-atmosphere pathway. All the other surface processes that we could envision were, to the best of our present investigation capability, not able to supply enough carbon to the atmosphere. During ice melt, the model is in agreement with the much lower K values, suggesting that air-ice CO2 fluxes are dominated by the dissolved pathway, apart from the short under-sampled transition period when other processes such as episodic convection events might have taken place.

Using our best model run (CTRL), we compared the contribution of all processes based on the computation of the model DIC budget, following Moreau et al. (2015; see Figure 6). Sea ice is gaining DIC throughout its growth (Figure 6a). Large quantities of DIC are incorporated into sea ice by growth, although 80% of this DIC is rejected through brine drainage (Figure 6b). Ikaite precipitation does not significantly contribute to the budget of DIC. Ice-atmosphere CO2 fluxes (Figure 6c) are dominated by the escape of gas bubbles during ice growth. However, only the dissolved pathway affects CO2 fluxes during ice melt. As suggested by Moreau et al. (2015) the budget of DIC in sea ice is driven mainly by physics while biogeochemical processes (only chemical here), although significant, are secondary.

doi: 10.12952/journal.elementa.000112.f006.
Figure 6.

#### CO2 budget.

a) Daily budget of vertically integrated DIC (TDIC, mmol m−2) dissolved into CO2 (black), CO32- (light grey), and HCO3- (grey) for the CTRL run. b) Corresponding daily mean changes (mmol m−2 day−1) in TDIC due to total sea ice growth and melt (dark blue), brine drainage (light blue), and ikaite precipitation/dissolution (red). c) Corresponding daily mean changes (mmol m−2 day−1) in TDIC due to diffusive CO2 fluxes (light blue) and gas bubble fluxes (dark blue).

Based on these considerations, we expand to ice-atmosphere exchanges of CO2 the suggestion that gas bubbles provide important contributions to the stocks of gaseous compounds in sea ice (Mock, 2002; Tison et al., 2002; Rysgaard and Glud, 2004; Zhou et al., 2013; Moreau et al., 2014). The conceptual view that we propose is based on three key arguments. 1. Gas bubbles easily form in sea ice because brine shrinking with temperature induces drastic increases in gas concentrations, and because there is a net decrease in gas solubility with decreasing temperature in brine inclusions. 2. Gas bubbles in the liquid brine rise upward if the connectivity of the brine network is sufficient. 3. The dilution of brines during warming brings CO2 concentration in brine below saturation, which brings trapped gas bubbles back into dissolved form. Because there is no analytical method to assess the partitioning of gases between dissolved and gaseous forms in sea ice brine, it is difficult for the time being to confirm this conceptual view directly. In natural sea ice, other surface processes likely affect air-ice CO2 fluxes as well (e.g., the formation of frost flowers, brine skim; Geilfus et al., 2013; Barber et al., 2014). In the present experiments, however, they were not observed.

Finally let us stress that there are only a few processes able to sustain enhanced CO2 fluxes during a long period. The presence of frost flowers would be transient and have only short-term effects on CO2 fluxes. Other processes (e.g., melting/freezing cycle of the ice surface) modulate air-ice CO2 fluxes, but few of them can actually sustain CO2 fluxes for extended periods of time as we observed during the experiment. Processes that sustain CO2 for extended periods of time must involve either the production of CO2, like bacterial respiration, or the transfer of CO2 from a large reservoir (like the ice interior or the ocean) to the surface.

## 4. Conclusions

The first aim of this study was to determine experimentally the air-ice CO2 transfer coefficient from continuous CO2 flux measurements during an ice growth and decay cycle in an ice-tank experiment. Discrete measurements of air-ice pCO2 gradients and DIC anomalies reflected well the amplitude and the patterns of the air-ice CO2 fluxes, supporting the reliability of our methods and results, including the calculation of a bulk gas transfer coefficient. The second aim was to discriminate the different drivers of air-ice CO2 fluxes using a 1D sea ice carbon cycle model (Vancoppenolle et al., 2010; Moreau et al., 2015), including explicit empirical representations of dissolved gas and gas bubble ice-atmosphere pathways and testing for several gases (Ar, O2, CO2) at a few locations.

There are three key findings in the paper. 1. The observation-based gas transfer coefficient, retrieved by dividing the observed CO2 flux by the air-brine pCO2 difference, was ∼6 times higher during growth (K = 2.5 mol m−2 d−1 atm−1) than during melt (K = 0.4 mol m−2 d−1 atm−1). 2. The time evolution of the sign of the air-ice CO2 flux, characterized by an efflux from the ice during growth and an influx during melting, was consistent between observations and the nine model simulations, as well as with previous literature (Nomura et al., 2010a; Geilfus et al., 2012a; Delille et al., 2014; Geilfus et al., 2015). Such evolution must therefore be seen as typical and robust. 3. The magnitude of the observed CO2 flux is consistent with previous literature but not between the different model simulations. The simulated dissolved CO2 flux clearly underestimates the observed value because of an intrinsic limitation by DIC stocks in the model already identified by Moreau et al. (2015). The observed magnitude of the ice-air CO2 flux can only be reached in the model by invoking an intense escape of gas bubbles through an ice-atmosphere pathway. Such an intense pathway was achieved by tuning the gas bubble nucleation time scale down to a few hours, which contrasts with the value of a few weeks used in a previous 6-month simulation of Argon dynamics in natural sea ice (Moreau et al., 2014).

Based on these concordant findings, we infer that the gas bubble ice-atmosphere pathway is likely a significant contributor to air-ice flux in young growing sea ice, next to diffusion. This intense gas bubble pathway in growing sea ice proposed here is plausible at this stage but subject to caution, because it does not rely on direct process observations, but rather on inference from air-ice CO2 flux and a single sea ice model, characterized by a number of assumptions. Further evaluation of the proposed scenario and reduction of uncertainties rely on improvements of both observation and modelling techniques. Observations of gas bubbles from X-ray tomography (Crabeck et al., 2016) provide new insights on gas bubble size and number distribution, and might be used to document the upward migration of gas bubbles. Future models should include at least pressure effects on gas bubble nucleation and perhaps tortuosity effects on gas bubble rise, as well as improved representations of brine dynamics (e.g., Griewank and Notz, 2013; Turner et al., 2013; Rees Jones and Worster, 2014).

The multiphase nature of CO2 in sea ice is well established, with contributions of dissolved species, gas bubbles (Tison et al., 2002; Zhou et al., 2013; Crabeck et al., 2014a; Moreau et al., 2014, 2015) and ikaite crystals (Papadimitriou et al., 2004; Dieckmann et al., 2008), a solid inorganic compound that contains inorganic carbon. This multiphase nature of CO2 is, to our best knowledge, a unique feature in marine environments. Clarifying the partitioning of CO2 between the gaseous and dissolved phases and the gas transport pathways, as well as the contribution of surface ice processes, is a critical challenge for future research related to gas dynamics and the carbon cycle in ice-covered seas.

## Data accessibility statement

The observation data and the model code will be made available to any scientist upon request and which is encouraged in order to reproduce the computations and test the model results in other locations (contacts: marie.kotovitch@ulg.ac.be, s.moreau@uclouvain.be and bruno.delille@ulg.ac.be).

© 2016 Kotovitch et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

## Contributions

Lead of the INTERICE V experiment: DT

Contributed to planning and experimental design: DT, GD, K-UE, J-LT, BD

Contributed to carry out the experiment and acquisition of data: JZ, GD, K-UE, DT, J-LT, BD

Contributed to analysis and interpretation of data: MK, SM, JZ, FVdL, J-LT, BD, MV

Contributed to the model simulations and interpretation: SM, MV

MK wrote the paper, with the help of SM, JZ, BD, MV and comments and inputs from all coauthors.

## Competing interests

The authors have declared that no competing interests exist.

## Funding information

This study was supported by the European Community’s 7th Framework Programme through the grant to the budget of the Integrated Infrastructure Initiative HYDRALAB-IV, Contract no. 261520. The work was supported by a FiDiPro award by the Academy of Finland, the Walter and Andree Nottbeck Foundation, and the BIGSOUTH project funded by the Belgian Science Federal Policy Office. MK and JZ, SM and BD are respectively research fellows, postdoctoral researcher and research associate of the Fonds de la Recherche Scientique – FNRS. MV acknowledges support from the EU, through BISICLO (FP7 CIG grant no 321938). This is a MARE contribution.