Air-ice carbon pathways inferred from a sea ice tank experiment

Among the natural systems, the role of sea ice in CO₂ cycling is not well constrained. Several studies report air-ice CO₂ 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 CO₂ cycling are needed to assess the future capacity of polar oceans to buffer the rise of atmospheric CO₂ 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 CO₂ 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 CO₂ fluxes were monitored using the chamber method. The CO₂ transfer coefficient during both ice growth and ice decay was computed and compared to a sea ice model integrating carbon dynamics.

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 m³ 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 CO₂ 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 CO₂ (pCO₂).

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 pCO₂ (denoted as pCO_{2 bulk}) within permeable sea ice. The goal is to equilibrate the sea ice samples with a N₂/CO₂ 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 CO₂) 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 CO₂ concentration. Shortly afterward, the sample temperature was measured to check for experimental drift. Measured pCO_{2 bulk} was corrected for the temperature difference between the sample and the in situ temperature according to the corrections proposed by Copin-Montegut (1988).

We measured total alkalinity (TA) on melted sea ice (TA_bulk) 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 HgCl₂ 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⁻¹.

We computed the dissolved inorganic carbon of bulk ice from TA_bulk and pCO_{2 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 TA_brines 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 DIC_brines) from TA_brines and pCO_{2 bulk} using the CO2SYS program for the carbonate system (Lewis and Wallace, 1998). We used the CO₂ dissociation constants of Mehrbach (1973) refitted by Dickson and Millero (1987) and the other constants advocated by Dickson and Goyet (1994). We then converted DIC_brines to DIC_bulk 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 CO₂ 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.

To measure the underlying seawater pCO₂ 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 pCO₂ 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⁻¹ and a closed air loop ensured circulation through the equilibrator and the IRGA at a rate of 3 L min⁻¹. The IRGA was calibrated before and after the experiment with N₂ and CO₂:N₂ 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 N₂. Uncertainty during this experiment was less than 6 µatm.

In this paper, positive CO₂ flux refers to CO₂ flux from the ice to the atmosphere, while negative CO₂ flux refers to a flux from the atmosphere to the ice. We measured air-ice CO₂ 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 pCO₂ 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⁻¹. When the pCO₂ of ice is higher than atmospheric pCO₂, CO₂ 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 CO₂ 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).

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 CO₂ (K, in mmol m⁻² d⁻¹ atm⁻¹), using the equation developed by Liss and Slater (1974) and Sarmiento and Gruber (2004):

where F is the air-ice CO₂ flux in mmol m⁻² d⁻¹, pCO_{2 bulk} is the pCO₂ in the ice and pCO_{2 air} is the pCO₂ 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).

A Bootstrap resampling statistical analysis procedure, using random values of the measured parameters (temperature, salinity, pCO₂, 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.

In order to understand the factors that drive air-ice CO₂ fluxes during the experiment, we ran a one-dimensional thermodynamic sea ice model representing sea ice halo-thermodynamics and carbon dynamics, including air-ice CO₂ 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 CO₂ flux, F^CO₂ (mmol m⁻² d⁻¹), which is assumed proportional to the CO₂ partial pressure (pCO₂) 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:

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), fCO2 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):

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

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 CO₂ from brine to the bubble compartment as a function of the CO₂ super-saturation. At each time step, a fraction R^bub of the CO₂ 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., e_T^gas = 0.07, Zhou et al., 2013; Moreau et al., 2014). If the fraction of sea ice with e > e_T^gas includes the ice surface, all gas bubbles escape to the atmosphere and contribute to the air-ice CO₂ flux. For our CTRL simulation, we used R^bub = 10% h⁻¹ and e_T^gas = 0.07 (Table 1).

In order to understand the potential reasons to explain the underestimation of the observed air-ice CO₂ fluxes by the model, we performed three series of sensitivity experiments (Table 1). First, we tested the impact of a more intense dissolved CO₂ pathway, by changing D_diff and z_BL (runs 2–4). The D_diff value of Broecker and Peng (1974) is derived from seawater and only includes diffusive effects. We tested the slightly higher D_diff value of Loose et al. (2011), derived from sea ice experiments, which can potentially include both diffuse and bubble contributions. A lower z_BL (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 R^bub, 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⁻³ and 2090 mmol m⁻³, and sea ice TA and DIC concentrations were set to 850 mmol m⁻³ and 750 mmol m⁻³.

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.

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 (r² = 0.95; Figure 1).

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 (CaCO₃·6H₂O). 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⁻¹, whereas the error on TA measurement was 4 µmol kg⁻¹.(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⁻¹ in Dieckmann et al., 2008; 15–19 µmol kg⁻¹ in Geilfus et al., 2013; 100–900 µmol kg⁻¹ in Rysgaard et al., 2013). This evidence indicates that if ikaite precipitated during the experiment, it was not significant.

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

Although the chamber was not properly sealed between day 0 and day 5, air-ice CO₂ 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⁻² d⁻¹; Nomura et al., 2006) and slightly lower than measurements over natural sea ice. Delille et al. (2014) measured CO₂ fluxes ranging from −5.2 mmol m⁻² d⁻¹ to 1.9 mmol m⁻² d⁻¹ on Antarctic pack ice in spring, and ascribed these fluxes to seasonal pCO₂ gradients between the brine and the atmosphere. Geilfus et al. (2012a) measured CO₂ fluxes at the sea ice interface ranging from −2.63 mmol m⁻² d⁻¹ up to 0.84 mmol m⁻² d⁻¹ in the Arctic coastal zone, while Nomura et al. (2010a) measured CO₂ fluxes ranging from -1 mmol m⁻² d⁻¹ to 0.7 mmol m⁻² d⁻¹ over land fast ice in Barrow at the end of spring.

The temperature and the pCO_2air 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 CO₂ fluxes is clearly opposite and strongly affected by the rise of atmospheric pCO₂ 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 CO₂ fluxes, and the overall pattern of observed air-ice CO₂ fluxes is consistent with in situ observation of a seasonal cycle, with upward CO₂ fluxes during ice growth (Nomura et al., 2010a; Geilfus et al., 2012a; Delille et al., 2014) and downward CO₂ fluxes during ice melt (Delille et al., 2014; Geilfus et al., 2015).

In parallel with direct measurements of air-ice CO₂ fluxes, we derived air-ice CO₂ 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 DIC₇) to remove the salinity-related changes (brine rejection, concentration and dilution; Figure 3). If no biogeochemical processes occurred (biological activity, CaCO₃ precipitation, CO₂ transfer to the gas phase and CO₂ exchange with the atmosphere), DIC₇ 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 DIC₇ depletion in the top 12.5 cm of the growing ice was due to CO₂ release from the super-saturated ice to the atmosphere above. Whereas, during the melting phase we observed a DIC₇ increase in the top 12.5 cm. This increase can be linked to the downward air-ice CO₂ flux measured during the melting phase.

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

where DIC_{th, i} and DIC_{12.5 cm} are computed for each sampling event and expressed in µmol kg⁻¹. 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 (ΔDIC_t in mmol m⁻²):

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

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⁻³. Based on the precision of the measured variables (salinity, temperature, TA, pCO₂; see Methods section), the precision on F is assessed at ± 0.03 mmol m⁻² d⁻¹ (which represents an error of 15% for the averaged CO₂ flux of 0.2 mmol m⁻² d⁻¹).

From these DIC anomalies, the calculated air-ice CO₂ fluxes (black triangles in Figure 2a) are in good agreement with the observed air-ice CO₂ fluxes, except on day 16 where the difference represents more than 1 mmol m⁻² d⁻¹; 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 CO₂ fluxes when snow cover is absent, despite the large variation of the atmospheric pCO₂.

Using equation 1 for gas exchange, we calculated a gas transfer coefficient for CO₂ at the air-ice interface for each sampling day of the experiment (see K; Table 2), where F corresponds to the air-ice CO₂ fluxes measured with the automated chamber, pCO_{2 air} to the pCO₂ at the air interface and pCO_{2 bulk} to the pCO₂ 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⁻² d⁻¹ atm⁻¹. A Bootstrap resampling statistical analysis of the propagation of the uncertainties of the measured variables (salinity, temperature, measured CO₂ fluxes, pCO_{2 bulk} and pCO_{2 air}) gives an uncertainty for K of 0.9 mol m⁻² d⁻¹ atm⁻¹ where standard error deviation is only 0.25 mol m⁻² d⁻¹ atm⁻¹. For the melting phase (which corresponds to two measurements on day 19) the mean K value was 0.4 mol m⁻² d⁻¹ atm⁻¹ (the uncertainty derived from a Bootstrap resampling statistical analysis was 0.08 mol m⁻² d⁻¹ atm⁻¹). K was therefore 6 times higher during the growth phase than during the melting phase.

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:

where Sol (mmol m⁻³ atm⁻¹) is the solubility of CO₂ in salt water (here brine), a function of temperature (T) and salinity (S), following Weiss (1974):

We used the T and S of sea ice brine. The resulting value, k = 0.164 (± 10%) cm h⁻¹, 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 CO₂ 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 CO₂ 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.

Here we use the model to frame the CO₂ 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 CO₂ flux is consistent between observations and all model simulations. All model runs show an efflux of CO₂ 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 pCO₂ 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 CO₂ fluxes are reproduced. The model gives a good estimation of downward CO₂ fluxes during ice melting although the time of the onset of downward CO₂ fluxes occurs earlier than observed in the model.

By contrast, the magnitude of the air-ice CO₂ 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 D_diff, boundary layer thickness z_BL, vertical resolution, ikaite, etc.). In turn, the total (dissolved + bubble) air-ice CO₂ flux is essentially determined by the gas bubble nucleation rate R^bub, which acts as a bottleneck parameter for the efflux of CO₂ to the atmosphere via the gas bubble pathway. This bottleneck effect explains why R^bub was used as a tuning parameter, set to 10% h⁻¹ 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 z_BL and D_diff could increase the upward CO₂ fluxes during ice growth. None of our trials was successful. For example, increasing D_diff up to 2.4 x 10⁻⁸ m² s⁻¹ (run 2: Loose-D_diff) barely affects the upward CO₂ fluxes during ice growth. Decreasing z_BL down to 0.05 µm did not significantly increase the simulated air-ice CO₂ fluxes either. This model behavior arises because the ultimate limitation for the air-ice CO₂ diffusive flux is the near surface stock of DIC: supply to the near-surface DIC is not rapid enough to sustain CO₂ fluxes to the observed magnitude. This result was already obtained by Moreau et al. (2015; see their Figure 9b where the simulated outward CO₂ fluxes during winter at Barrow reach asymptotically 1.3 mmol m⁻² d⁻¹ when z_BL is decreased down to 10⁻³ µ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 R^bub (10-20% h⁻¹). Yet, the bubble nucleation rate R^bub is far from being well constrained from observations. In Moreau et al. (2014), this value was tuned to 0.1% h⁻¹, 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 CO₂ 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 R^bub values of 10–20% h⁻¹ 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 CO₂ 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 CO₂ fluxes. The precipitation of ikaite could have influenced ice-atmosphere CO₂ 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⁻¹. 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 CO₂ flux (Figure 5b).

We also tested the role of vertical resolution of the simulated CO₂ 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 pCO₂, which slightly intensifies the dissolved gas pathway during the growth phase.

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 CO₂ 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 CO₂ fluxes (Figure 6c) are dominated by the escape of gas bubbles during ice growth. However, only the dissolved pathway affects CO₂ 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.

Based on these considerations, we expand to ice-atmosphere exchanges of CO₂ 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 CO₂ 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 CO₂ 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 CO₂ fluxes during a long period. The presence of frost flowers would be transient and have only short-term effects on CO₂ fluxes. Other processes (e.g., melting/freezing cycle of the ice surface) modulate air-ice CO₂ fluxes, but few of them can actually sustain CO₂ fluxes for extended periods of time as we observed during the experiment. Processes that sustain CO₂ for extended periods of time must involve either the production of CO₂, like bacterial respiration, or the transfer of CO₂ from a large reservoir (like the ice interior or the ocean) to the surface.

The first aim of this study was to determine experimentally the air-ice CO₂ transfer coefficient from continuous CO₂ flux measurements during an ice growth and decay cycle in an ice-tank experiment. Discrete measurements of air-ice pCO₂ gradients and DIC anomalies reflected well the amplitude and the patterns of the air-ice CO₂ 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 CO₂ 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, O₂, CO₂) at a few locations.

There are three key findings in the paper. 1. The observation-based gas transfer coefficient, retrieved by dividing the observed CO₂ flux by the air-brine pCO₂ difference, was ∼6 times higher during growth (K = 2.5 mol m⁻² d⁻¹ atm⁻¹) than during melt (K = 0.4 mol m⁻² d⁻¹ atm⁻¹). 2. The time evolution of the sign of the air-ice CO₂ 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 CO₂ flux is consistent with previous literature but not between the different model simulations. The simulated dissolved CO₂ 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 CO₂ 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 CO₂ 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 CO₂ 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 CO₂ is, to our best knowledge, a unique feature in marine environments. Clarifying the partitioning of CO₂ 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.

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: [email protected], [email protected] and [email protected]).

© 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.

The authors would like to thank the Hamburg Ship Model Basin (HSVA), and the ice tank crew, for the hospitality, technical and scientific support and the professional execution of the test program in the Research Infrastructure ARCTICLAB. We are grateful to Michael Fischer for his help with the deployment of the automated chambers and treatment of flux data, to Nicolas-Xavier Geilfus, one anonymous reviewer, the Editor Jody Deming and the associate editor for their comments that improve the overall quality of the manuscript.