A Parameter Model of Gas Exchange for the Seasonal Sea Ice Zone

. Carbon budgets for the polar oceans require better constraint on air–sea gas exchange in the sea ice zone (SIZ). Here, we utilize advances in the theory of turbulence, mixing and air–sea ﬂux in the ice–ocean boundary layer (IOBL) to formulate a simple model for gas exchange when the surface ocean is partially covered by sea ice. The gas transfer velocity ( k) is related to shear-driven and convection-driven turbulence in the aqueous mass boundary layer, and to the mean-squared wave slope at the air–sea interface. We use the model to estimate k along the drift track of ice-tethered proﬁlers (ITPs) in the Arctic. Individual estimates of daily-averaged k from ITP drifts ranged between 1.1 and 22 m d − 1 , and the fraction of open water ( f ) ranged from 0 to 0.83. Converted to area-weighted effective transfer velocities ( k eff ) , the minimum value of k eff was 10 − 5 m d − 1 near f = 0 with values exceeding k eff = 5 m d − 1 at f = 0 . 4. The model indicates that effects from shear and convection in the sea ice zone contribute an additional 40 % to the magnitude of k eff , beyond what would be predicted from an estimate of k eff based solely upon a wind speed parameterization. Although the ul-timate scaling relationship for gas exchange in the sea ice zone will require validation in laboratory and ﬁeld studies, the basic parameter model described here demonstrates that it is feasible to formulate estimates of k based upon properties of the IOBL using data sources that presently exist.


Introduction
The regulation of air-water gas exchange occurs within a small vertical region of ∼ 100 µm at the air-water interface -a region that is too small for direct measurement of the properties governing this exchange, namely partial pressure and diffusive length scale. In addition, the air-water interface is continually being deformed by waves and surfactants and consequently the profile of turbulence and gas abundance are difficult to determine for a given surface condition (Jähne et al., 1987). This limitation at the critical scale has led to a proliferation of rate approximations based on quantities that can be more readily measured. Wind speed is the most common of these quantities, as it is an important mixing mechanism for many surface water bodies, including the ocean (Liss, 1973;Wanninkhof, 1992;Wanninkhof and McGillis, 1999;Nightingale et al., 2000;Ho et al., 2006Ho et al., , 2011Sweeney et al., 2007;Wanninkhof et al., 2009). However, even for the blue water ocean, wind speed parameterizations account for 75 % or less of the observed variability in the gas exchange rate (Frew et al., 2004;Stanley et al., 2006). Wind speed parameterizations are widely used where wind is the first order mixing mechanism, and may be adequate under certain open ocean conditions.
In contrast to the open ocean, there are physical systems such as wetlands, rivers, estuaries, and the coastal surf zone where wind speed derived estimates of gas exchange are incomplete given the other mechanisms that produce aqueous turbulence. These mechanisms include tides, internal waves B. Loose et al.: A parameter model of gas exchange and boundary mixing. In these systems, it is necessary to rely on additional or alternative metrics for relating gas exchange to the underlying forcing conditions (Ho et al., 1997;Zappa et al., 2003;Frew et al., 2004).
For moderately to sparingly soluble (i.e., most) gases, the existence of an aqueous viscous sublayer at the air-sea interface is the dominant restriction to transfer across the interface. This bottleneck ranges in length scale from 20 to 200 µm (Jähne and Haubecker, 1998); turbulent eddies impinging on the free surface can reduce or collapse the viscous sublayer, leading to "bursts" of air-sea gas transfer (Rao et al., 1971). Consequently, measurement of the TKE budget close to the air-sea interface provides the best metric for determining the gas exchange rate when multiple drivers are at work.
We have argued that the sea ice covered polar oceans are another region where wind speed parameterizations are not expected to be adequate (Loose et al., 2009;Loose and Schlosser, 2011), and in this contribution we present an alternative parameter model that is based on the mechanisms that produce turbulence in the ice-ocean boundary layer (IOBL): (1) shear caused by the velocity differential between drifting ice and the water column beneath (McPhee and Martinson, 1992), (2) buoyant convection caused by heat loss and phase changes at the ocean surface (Morison et al., 1992), and (3) short period waves and wave interactions with ice floes (Wadhams et al., 1986;Kohout and Meylan, 2008) (Fig. 1). Our goal is to explore what effect these turbulent production mechanisms may have on the magnitude of k by utilizing the significant advances in theory and observation of turbulence in the IOBL, developed over the past 30 yr (McPhee, 2008). This exercise is meant to provide first-order insight into the processes by which sea ice might simultaneously produce and modulate gas exchange, and to provide a model template upon which laboratory and field studies can begin to establish firm quantitation of these processes.
In Sect. 2 we introduce the basic theory that has been developed to relate differing turbulence production mechanisms to the rate of air-sea gas exchange. Section 3 briefly reviews how to estimate turbulence production from shear in the IOBL as a function of wind speed, ice drift velocity and sea ice cover. Section 4 contains a similar review for the estimate of turbulence production by convection in the IOBL. Section 5 provides a rough estimate of k from capillarygravity waves. Section 6 is used to present estimates of k, calculated from four ice-tethered profiler (ITP) missions. Finally we discuss the results from the ITP data and the parameter model and how they might relate to future investigations of k in the sea ice zone (SIZ).  Fig. 1. A graphic illustration of three mechanisms that can lead to turbulence production and gas exchange in the ice-ocean boundary layer.

A note about upscaling k in the sea ice zone
The sea ice zone exists as a heterogeneous distribution of ice cover and open water, and this presents a challenge when trying to appropriately represent a gas transfer rate that takes this heterogeneity into account. Because sea ice is a porous medium (Cox and Weeks, 1983), it is possible in principle to ventilate the surface ocean of gas through the ice matrix as across the air-sea interface. To reflect this, the air-sea gas flux can be written as where C is the molar gas differential across the air-sea interface, and k eff is the "effective" gas transfer velocity, where f is the fraction of open water, (1 − f ) is the fraction of sea ice cover, k ice represents the bulk diffusive gas flux through the sea ice D/z ice , and k is the open water gas transfer velocity -what we are most accustomed to working with. This representation of k eff is consistent with the asymptotic homogenization of diffusive transport through a medium with rapidly varying permeability, such as the ice covered upper ocean (Dykaar and Kitandis, 1992;Jähne and Haubecker, 1998;Golden et al., 2006). In practice, the rate of gas diffusion through cold columnar sea ice is within a factor of 10 of molecular diffusion in water (e.g., ∼ 10 −4 cm 2 s −1 ). Thus k ice is very small compared to k -at least for cold, consolidated columnar sea ice (Gosink et al., 1976;Loose et al., 2010). In the presence of thin ice, such as grease ice, nilas or frazil ice, gas exchange may not be so slow, but solute rejection from the ice (Killawee et al., 1998;Loose et al., 2009) makes this circumstance hard to interpret. However, the distinction, between k and k eff is worth clarifying for another reason. There are a variety of methods for estimating the air-sea gas transfer velocity. Based on the scale over which these methods integrate, some yield k eff and others yield k or k ice , and it is important to be able to put the results from these methods on a common scale. For example, tracer-based estimates of gas transfer velocity such as the radon deficit method (Smethie et al., 1985;Bender et al., 2011), and the dual-tracer method (Wanninkhof et al., 1993;Nightingale et al., 2000) will yield estimates of k eff , whereas covariance flux (McGillis et al., 2001), and gradient flux (Orton et al., 2010) methods represent a much smaller spatial footprint and will yield a gas transfer velocity that will represent k ice or k, but not k eff . Examples of past studies that estimate k eff are Fanning and Torres (1991), Loose et al. (2009), andLoose and. The study by Miller et al. (2011) would yield estimates of k ice and the study by Else et al. (2011) would yield estimates of k for a polynya-type water surface.
In this study, the parameter model estimates of gas transfer are based on predictions of the aqueous turbulence and represent values of k. We will present the results as k and k eff , using Eq. (2) and by assuming that k ice is negligible.

TKE dissipation model
We require an estimate of k that accounts for turbulence production mechanisms that are unique to the SIZ. This can be accomplished using the relationship between TKE dissipation (ε) and gas exchange (Lamont and Scott, 1970), where ν is the kinematic viscosity, and Sc is the Schmidt number. For the purposes of this paper, we will assume the constant of proportionality determined by Zappa et al. (2007); 0.419. It should be noted that ε varies vertically in the water column and is likely to assume distinct values at the ice-water and air-water interfaces, therefore unique proportionality constants may exist for the SIZ. Next we assume that over the mixing timescale of the surface ocean TKE production should balance dissipation, leading to a stationary form for the kinetic energy equation (Gaspar et al., 1990), The first term represents TKE production by shear and the second represents TKE production by buoyancy flux. Monin-Obukhov similarity theory (MOST) demonstrates how to relate turbulence produced by convection and turbulence produced by shear (Lombardo and Gregg, 1989). When convection ceases and stratification sets in, MOST also provides a means for estimating the dampening effect on shear driven turbulence based on the Obukhov length scale, where u * 0 is the aqueous friction velocity at the ice-water or air-water interface, w b 0 is the buoyancy flux at the ocean surface, and κ is the Von Karmann constant.
Finally, ε from Eq. (4) can be related to the turbulent friction velocity in the water (u * 0 ) by where z the depth scale, assuming z lies within the boundary layer. By definition, u * 0 is a representation of the shear stress (τ ) that acts on the boundary layer, and it is the shear stress balance that is altered by the presence of ice floes, individually and as part of the ice pack. In the next section we will review and describe the stress balance in partially ice-covered seas and then use the above equations to make estimates of k.

Shear in the IOBL
Sea ice floes present additional surfaces against which normal and tangential forces can operate. An individual ice floe is subject to surface drag from wind over the ice, which is transmitted to the IOBL via bottom drag, when the ice and the water beneath move relative to each other (Shirasawa and Ingram, 1991). In addition to surface stresses, an ice floe experiences form drag as it displaces water by its lateral motion (Perrie and Hu, 1997), and by wave radiation stress as gravity waves propagate from open ocean into the pack Hu, 1996, 1997). Steele et al. (1989) observed that as sea ice concentration decreases, the ice pack behaves more like a "flotilla of barges" where individual floes may experience form drag that exceeds surface drag, when the diameter of the floe is less than 300 times the ice draft, based on the ratio of surface and form drag coefficients. The total shear stress exerted on the aqueous boundary layer can be estimated as a function of sea ice cover (1 − f ), wind speed (U ) and ice drift velocity (V 0 ). By these mechanisms the total drag on the water surface boundary layer is as follows (Steele et al., 1989): where τ skin-iw is the ice-water skin drag, τ form is the form drag around ice floes, τ skin-aw is the air-water skin friction and W is the fraction of air-water stress that is used to produce surface gravity waves (Steele et al., 1989). W = 1.6 × 10 −4 ρ w /ρ a is based on the Charnock (1981) relationship.
(τ form ) is dependent on the profile the ice presents to the flow as it forces its way through the water. This effect was represented by Steele et al. (1989) by incorporating the freeboard (D) and length (L) of the ice floe such that form drag decreases as the characteristic dimension of the flow size increases. This behavior is meant to simulate the "wake effect" where the overall pressure gradient from the front to the back www.ocean-sci.net/10/17/2014/ Ocean Sci., 10, 17-28, 2014 of the floe is decreased as L increases.
where V 0 is the relative velocity between the water in the mixed layer and the ice floating on it (Lu et al., 2011), and ρ w is the water density. Tennekes and Lumley (1972) showed that the form drag decreases approximately proportional to L −0.5 . In general, when sea ice concentration is close to 0 %, the characteristic dimension of L is small, e.g., 10 m. As ice cover increases, so does L (Toyota et al., 2006), and the relationship follows a power-law dependence of the form N(L) ∝ L −α , where N is the number density of ice floes per km 2 (Rothrock and Thorndike, 1985). Toyota et al. (2006) observed that this relationship cannot be fit with a single value for α; when L is less than 40 m, the value of α is 1.15 and when L is greater than 40 m, α = 1.87 provides the best fit. N (L) can be related to the fraction of sea ice cover if we assume a characteristic shape for ice floes. The SIC was computed as the average of two "ideal" floes, one squareshaped and one circular, with sides and diameters of length L, multiplied by N (L) (Fig. 2).
(τ skin-iw ), the next term in the boundary layer stress balance (Eq. 8), represents the skin friction at the ice-water interface. The ice-ocean boundary layer is often compared with the atmospheric boundary layer as both are affected by planetary rotation, both experience relatively constant stressdriven turbulence (Blackadar and Tennekes, 1968) and both tend to follow the law of the wall, which permits a complete determination of the velocity profile and turbulence closure in terms of two variables: u * 0 and the roughness length scale z o (McPhee, 2008), where U (z) is the Euclidean norm velocity at depth z in the surface boundary layer. Referring to Eq. (7), the skin friction can alternately be expressed as a velocity or a stress with no loss of generality. Observations in the IOBL show that the surface boundary layer, where velocity varies linearly with the log of depth (∼ 2 m) behaves differently from the outer layer (∼ 70 m) (Shaw et al., 2008). Outside of the surface layer, it is necessary to account for Ekman turning on the shear stress, and this can be accomplished using the Rossby similarity theory (Blackadar and Tennekes, 1968;McPhee, 2008). To determine the surface shear stress in terms of the ice velocity it is necessary to account for the geostrophic drift at the base of IOBL because the magnitude of shear and the magnitude of rotation are interdependent (McPhee, 2008). This can be accomplished using the Rossby similarity theory (Blackadar and Tennekes, 1968;McPhee, 2008), which relates V 0 with u * 0 , the interface friction velocity overlying a neutral geostrophic current (McPhee, 2008; Shaw et al.,  Toyota et al., (2006). N follows a distinct power law distribution when floe size 4 decreases below a threshold (Toyota et al., 2006).  Toyota et al. (2006). N follows a distinct power law distribution when floe size decreases below a threshold (Toyota et al., 2006).
where V 0 and u * 0 are the two-dimensional ice-drift and friction velocities, expressed as imaginary numbers, f is the Coriolis parameter, and A and B are matching coefficients. The term log u * f z 0 is sometimes referred to as the friction Rossby number, Ro. It can be estimated with values from Eq. (10) for the nonrotating law of the wall as an initial guess of the roughness, with z approximated as the sea ice thickness.
When stratification stabilizes the surface layer, shear turbulence is damped and surface turbulence is confined to a shallower region of the IOBL. Under these conditions, the similarity terms become dependent on a stability term (µ * = u * 0 κ/L 0 ); Ro (µ * ), A(µ * ), and B(µ * ), the solution to which is described in Sect. 4.2.3 of McPhee (2008).
The surface roughness z 0 is among the most variable and patchy IOBL parameters, and this variability is caused by the large range in ice thickness that can occur from thin ice in a lead (e.g., ∼ cm) to deep keels underneath sea ice pressure ridges (e.g., ∼ 10 m) (Shaw et al., 2008). Estimates of z 0 span two orders of magnitude from 1 mm to 9 cm, depending on ice type (McPhee, 2008). McPhee (1992) found roughness lengths varying from 2 to 9 cm during three Arctic drift experiments, and (Shaw et al., 2008) found z o roughness values near 10 cm using data from autonomous ocean flux buoys. Younger ice is smoother with z 0 close to 1 mm and in multiyear pack ice z 0 is about 5 cm (McPhee, 2008).
(τ skin-aw ), the final stress term in Eq. (8), refers to the shear stress exerted by the wind on the upper ocean. In the SIZ this stress occurs in the regions where water is exposed by divergent sea ice motions (i.e., leads and polynyas), or between floes in the marginal ice zone. The air-water shear stress that is transmitted into the ocean surface boundary layer can be expressed as where ρ a is the air density, C d is the air-water drag coefficient (1.5 × 10 −3 , Andreas and Murphy, 1986) and U is the 10 m wind speed. It is instructive to visualize how the individual terms of the total stress exerted on the IOBL vary with the sea ice cover, or conversely with the fraction of open water (f ). For this visualization, we use a constant U = 10 m s −1 and a constant V 0 = 0.12 m s −1 , based on the average drift velocity of the autonomous ocean flux buoy #3 (http://www. oc.nps.edu/~stanton/fluxbuoy/index.html). Using these values, we have computed the individual terms in Eq. (8), (1−f )τ form , (1−f )τ skin-iw , and (f )τ skin-aw , for a variable ice cover ranging from f = 0 (100 % ice cover) to f = 1 (0 % ice cover), (Fig. 3). For U = 10 m s −1 , the air-water stress leads to a τ skin-aw = 0.17 N m −2 and V 0 = 0.12 m s −1 results in τ skin-iw = 0.05 N m −2 . Individually, these terms are proportionately modulated by the fraction of ice cover and open water, respectively, and therefore τ skin-aw increases linearly as f increases and τ skin-iw follows the opposite trend, achieving its maximum value when f = 0. In comparison, the relationship between f and τ form is more complicated. We observe two tendencies, the first as τ form increases quadratically from f = 0 to f = 0.8, and then decreases again from f = 0.8 to f = 1. The physical explanation for the two modes results from the competing influences of ice cover and floe dimension on the form drag. At f = 0, the influence of the floe dimension (1/L ∼ 10 −4 m −1 ) strongly diminishes τ form . As the ice cover decreases so does L and we observe a rapid increase in τ form , as the "wake effect" of surrounding floes becomes less and each individual floe is subject to greater form drag. If we compute the gas transfer velocity from τ using Eqs. (8), (6) and then (3), we observe that shear-driven turbulence in the IOBL contributes between 1.6 and 2.0 m d −1 to the overall gas transfer velocity in the sea ice zone, for a wind speed of 10 m s −1 and an ice velocity of 0.12 m s −1 (Fig. 3).

Buoyant convection/stratification
As with the shear-driven turbulence, the goal is to make a determination of ε, which can then be used to estimate k. The first step is to calculate the Obukhov length scale, to determine whether convection (L 0 < 0) or stratification (L 0 > 0) characterizes the surface layer of the ocean (McPhee, 2008). In the SIZ, as contrasted with the lower latitudes, L 0 can be affected by freezing and melting of ice as well as latent and sensible and shortwave heat transfer. Here, its worth pointing out that sea ice freezing (melting) can concentrate (dilute) the gas content of the upper ocean through solute rejection (dissolution) (Killawee et al., 1998;Loose et al., 2009), and change the upper ocean air-sea gas partial pressure differential. For the purposes of this parameter model we will not consider the effect of solute rejection and dissolution.
L 0 from Eq. (5) can be rewritten as (Lombardo and Gregg, 1989) where J b0 is the surface buoyancy flux in W kg −1 or m 2 s −3 equivalently. For the IOBL, the major contributions to the buoyancy flux are sensible heat transfer, sea ice freeze-melt, and the salt balance from evaporation −precipitation, where α is the coefficient of thermal expansion, β is the coefficient of haline contraction, ρ is the density of seawater, E − P is the rate of evaporation less precipitation, S 0 is the surface seawater salinity in kg kg −1 , c p is the specific heat capacity, and g is the gravitational constant. Again, the terms  14) represents the surface freshening or salinification and we estimate it from the net latent heat flux, which is an NCEP Reanalysis product. J SH q and J LH q over open water can be estimated using the bulk approximation for scalar flux in the atmospheric boundary layer flux (Andreas and Murphy, 1986), where T w is the water surface temperature, C H -the surface drag coefficient -is approximated as 1.26 × 10 −3 , and Q s is the specific humidity at the dew point, defined by T w . Q z is the specific humidity at the reference height, z, which should be the same as the height of the wind speed measurement (Andreas et al., 2010).
In the absence of boundary layer flow, the TKE dissipation is directly equivalent to the heat flux: ε = J 0 b . In a mixed condition of convection + shear, a similarity scaling for ε is chosen using L 0 to characterize the relative balance between convection and shear (Lombardo and Gregg, 1989): L 0 z, the boundary layer thickness, indicates free convection and L 0 z indicates a shear-dominated regime.

Surface waves and interactions with sea ice
In addition to shear and convection, surface waves and their interaction with ice floes will also produce turbulence and gas exchange in the SIZ. As discussed in Sect. 2, the actions of wind-driven surface waves are the dominant mechanism by which aqueous turbulence leads to gas exchange (Jähne et al., 1987;Wanninkhof, 1992), and among these, the highest frequency gravity waves appear most effective at disturbing the viscous mass boundary layer and enhancing k (Bock et al., 1999;Frew et al., 2004Frew et al., , 2007. The mean squared wave slope s 2 of these short period gravity waves correlates well with k (Bock et al., 1999;Frew et al., 2004Frew et al., , 2007, however the gravity wave frequency distribution is variable and interacts strongly with sea ice (Wadhams et al., 1986;Hayes and Jenkins, 2007). To date, most studies of wave interactions with sea ice focus on much longer period waves, because they can be responsible for the breakup and melt of ice floes (Squire et al., 1995;Meylan et al., 1997;Kohout and Meylan, 2008), and relatively less attention has been paid to waves in the capillary-gravity range (40-800 rad m −1 ). Energy attenuation by wave interaction with sea ice is strongest for waves with a wave number of 1 rad m −1 or less (Polnikov and Lavrenov, 2007). Therefore, it might be expected that short period capillary-gravity waves are reflected with near 100 % efficiency, based on wave reflection studies (Squire and Williams, 2008).
Direct measurements of the capillary-gravity wave spectrum in the presence of sea ice are necessary to constrain the contribution of surface wave forcing of gas exchange, and we are not aware of any study that has explicitly examined the wave field in this range. In the mean time, we can apply the relationship between s 2 and 10 m wind speed from Frew et al. (2007) (their Fig. 1) to give an idea of how surface forcing contributes to the overall gas exchange. The choice of a parameterization based on s 2 instead of a relationship between wind speed and gas exchange can be justified by recognizing that Jähne et al. (1987) have shown that waves in the capillary-gravity range arise and decay within seconds, upon excitation by wind, indicating minimal fetch dependence. In contrast, it is not well known how the wind speed-gas exchange relationship varies with fetch.
In lieu of direct measurements of s 2 in the SIZ, we have used data from Frew et al. (2007) that relates the 10 m wind speed to s 2 100 40 , the mean squared wave slope from 40 to 100 rad m −1 . To capture the tendency between data s 2 100 40 and U 10 we fit the Frew et al. (2007) data with a third order polynomial with fit coefficients [7.7814×10 −6 , −0.0001624, 0.0015064, −0.0012018] and an R 2 value of 0.91 (Fig. 4). As above, U 10 from NCEP data was interpolated onto the ITP drift track. Finally, the values of s 2 100 40 were used to estimate k using k mss = 0.336 + 1.82 × 10 5 s 2 100 40 (16) (Frew et al., 2004). The relationship between s 2 100 40 and gas exchange appears most robust under short fetch conditions and where bubbles are not a major contributor to the exchange rate (Frew et al., 2004(Frew et al., , 2007. Although we do not have measurements to support this, anecdotal observations indicate that breaking waves and bubble penetration are minimal in the sea ice zone.

Computation of k from ITP data
We have calculated the gas transfer velocity using the parameter model for shear and convective turbulence in the IOBL, using data from ITPs together with satellite data products for the Arctic Ocean. The inputs to this parameter model (f , U , V 0 , Q) and profiles of T and S, were taken from post-processed ITP missions 3, 6, 7 and 8, (http://www.whoi. edu/itp). Daily estimates of sea ice cover were derived from the Bootstrap algorithm for NIMBUS-7 and SSM/I satellite radiometers (Comiso, 2007). Daily estimates of wind speed and surface air temperature were taken from NCEP Reanalysis data provided by NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, (http://www.esrl.noaa.gov/psd/). The inputs from gridded remotely sensed data were extracted by interpolating onto the recorded daily positions of the ITPs.
The calculation for k proceeded as follows.
1. To start, the equation for buoyancy flux (Eq. 14) is evaluated to determine whether convection or stratification is affecting water column turbulence. The 26 1 Figure 4. Modified from Frew et al., (2007), the relationship between the 10-m wind speed and 2 the mean-squared wave slope has been represented in this study by a cubic (solid black line).  Frew et al. (2007), the relationship between the 10 m wind speed and the mean-squared wave slope has been represented in this study as cubic (solid black line). threshold for buoyant convection is determined using the Obukhov length (L 0 ) criteria: H /L 0 > 1 (Lombardo and Gregg, 1989). H is half the mixed-layer depth, and the mixed-layer depth was found by locating the maximum value of N 2 in the top 50 m of the water column -typically H was 10-15 m.

Shear-driven turbulence in the IOBL is determined by
summing the individual terms in the stress balance (Eq. 8). To solve Eq. (9) for τ form , the floe dimension, L, is estimated from the empirical relationship in Fig. 1, as follows. The floe number density, N, is first computed for the range of L from 10 to 10 4 m, using the parameterizations presented by Toyota et al. (2006). Next, the discrete values of N and L are used to compute the sea ice cover, (1 − f ), by taking the average of square (L × L) and circular π 4 L 2 floe geometry. For a discrete value of sea ice cover from the ITP drift, a value for L is computed as a table lookup. τ skin-iw is computed using Eqs. (11) and (7) as described in Sect. 3: here, the similarity terms, Ro, A and B, in Eq. (11) are evaluated in terms of the buoyancy flux (J 0 b , Eq. 14). Finally, τ skin-aw is computed using the Reanalysis wind speed in Eq. (12). The total shear-stress τ from Eq. (8), and resulting friction velocity (u * 0 ) is then used to compute the shear-driven TKE dissipation (ε s , Eq. 6).
3. Finally, the aggregate total dissipation (ε) was computed as ε = 0.84 0.58J 0 b + 1.76ε s , and k is computed from ε or ε s using Eq. (3). In Figs. 5 and 6, k (labeled k shear ) was computed from ε s , and in Figs Red "+" symbols indicate estimates of k shear that have been affected by stratification. The speckled color indicates the vertical distribution of N 2 > 10 −4 s −2 , giving a rough approximation of the mixed-layer depth. Some data gaps exist in the ITP profiles.
We present estimates of k from shear-driven turbulence separately from convection-driven turbulence to compare their relative contributions to the overall gas transfer velocity. These separate terms are presented for two of the four ITP missions (ITP 3 and 8). The duration of ITP 3 was 786 days, beginning 24 August 2005 with f , the open water fraction, less than 0.1 until the last 100 days, when ITP 3 drifted into waters with f of nearly 0.3, near the center of the Beaufort Gyre. The maximum daily NCEP wind speed, U , was 17.2 m s −1 and the mean U was 5.8 m s −1 , whereas the maximum ice velocity, V 0 was 0.24 m s −1 and the average V 0 was 0.05 m s −1 (Fig. 5). Throughout the drift, the ice-water skin friction dominated the stress balance, because of the large fraction of ice cover. On average, the ice-water skin friction was 10 times greater than the form drag and more than 100 times greater than the air-water skin friction.
It is worth noting that between days 900 and 1000, decreased ice cover (increased f ) corresponded to an increase in τ form as dictated by Eq. (9). However, this effect was small compared to the coincident increase in τ skin-iw that resulted from an increase in stratification. McPhee (2008) describes this effect in Sect. 4.2.3 -assuming the same boundary forcing conditions, the total shear stress in the stratified water column is the same as the neutral water column. However, stratification confines the profile of shear stress to a shallower surface layer. This effectively "concentrates" the shear stress in the surface layer, leading to an enhancement in the apparent TKE dissipation and thus gas exchange. The maximum observed values during ITP 3 all occurred during the ice retreat period and can be associated with the coincident  The k shear that would result from shear-driven turbulence. Red "+" symbols indicate estimates of k shear that have been altered by stratification through the stability function (Eq. 11). The speckled color indicates the vertical distribution of N 2 > 10 −4 s −2 , giving a rough approximation of the mixed-layer depth. Some data gaps exist in the ITP profiles.
increase in both τ form and τ skin-iw (Fig. 5). The mean value daily k shear for ITP 3 was 1.3 m d −1 and the maximum was 7.8 m d −1 .
A somewhat similar picture emerges for ITP 8; the ITP 8 data extended over 784 days, beginning 12 August 2007. The mean of the daily NCEP wind speed during the ITP 8 drift was also 5.7 m s −1 , with a maximum daily average wind speed of 14.5 m s −1 . However, the average ice velocity was 0.075 m s −1 , 25 % larger than the average drift velocity of ITP 3, and the maximum was 0.35 m s −1 . Overall, the variance in the drift velocity was greater for ITP 8, as compared to ITP 3. Near day 1600 (ca. 275 days after starting its drift), ITP 8 drifted into progressively opening sea ice cover, until f increased to a maximum of 0.83, around day 1700 (Fig. 6). Similar to what was observed in ITP 3, the decrease in sea ice cover produced a decrease in L, based on the parameterization in Fig. 2 The parameter model estimates of k conv are on average 3 times smaller as compared to the estimates of k shear . These estimates depict a seasonal cycle that predictably mirrors the air-water temperature differential. This is the case for both ITP 3 and 8 (Figs. 7, 8). During the period of peak heat flux, k conv reached its peak above 1.2 m d −1 for ITP 3 and 1.5 m d −1 for ITP 8. These peaks in k conv tended to occur when f > 0 and when the air-water temperature differential strongly favored heat loss to the atmosphere, i.e., when fall cooling was leading to sea ice formation, but open water area yet remained (e.g., days 1700-1800, Fig. 8). Conversely, the periods of apparent decrease in sea ice cover and increase in air temperature during spring (e.g., days 1600-1700, Fig. 8), caused L 0 to increase and the convection threshold to exceed the Obukhov length criteria H /L 0 > 1. This causes k conv to be "turned off" in the parameter model.
Having combined the TKE production from both shear and convection and using those values of ε to estimate k from ε = 0.84 0.58J 0 b + 1.76ε s , we find that the mean value of k from shear + convection was 1.63 m d −1 and the maximum was 10.6 m d −1 . In general, k is a scale independent or intensive property so we would expect no strong trend between k and f . This lack of tendency is borne out in Fig. 9 (top panel), despite the areal dependence that has been introduced through Eqs. (8) and (14) for the shear and buoyancy budgets. However, there is one interesting feature between f = 0.15 and f = 0.30, with a discernable peak in k, which exceeds 4 m d −1 (Fig. 9, top panel). This increase in k is the result of shear and convection enhancements during the seasonal transition periods when sea ice is in advance or retreat. As reported above, stratification during ice melt enhances k from shear by "concentrating" the net stress closer to the air-sea interface. During periods of decreasing f (e.g., fall), the effect of convection in water that is still open also leads to an enhancement in k. If k from shear and convection is converted to k eff using Eq. (2), we find that the effective gas exchange across the heterogeneous ocean surface is less than 1 m d −1 from f = 0 to f = 0.4, and increases almost linearly across the range from f = 0 to f = 0.8, with a value of ∼ 2 m d −1 at f = 0.8 (Fig. 9,   To estimate the total magnitude of gas transfer in the sea ice zone, we have stacked k from the mean-squared slope relationship in Eq. (9) with the estimates of k from shear and convection (described in steps 1-6 above), displayed in Fig. 10. The k ∝ f line determined by the mean NCEP wind speed for the Arctic, and prior values of k eff from Fanning and Torres (1991), and Loose and Schlosser (2011) have been overlaid to compare how this parameter model compares with those values. The parameter model exceeds the linear proportionality, but no values meet or exceed the median estimates of Fanning and Torres (1991). On average, the magnitude of k eff from ice processes is approximately 42 % of the total magnitude of k eff as can be observed in Fig. 10.

Discussion and conclusions
The objective of this study was to utilize the existing physical understanding of turbulence in the ice-ocean boundary layer (IOBL) to formulate a plausible parameter model for air-sea gas transfer in the seasonal sea ice zone. The model requires only six input parameters: sea ice concentration, wind speed, sea ice-drift velocity, air temperature and profiles of water temperature and salinity. These inputs are available from model output Timmermans et al., 2011), from satellite data products for the polar oceans, and from Arctic Observing Network nodes, such as ITPs and autonomous ocean flux buoys. Here, we estimated the magnitude of k that results from shear and buoyancy production processes that are unique to the sea ice covered ocean. We used input data from four ITP buoy missions: ITP 3, 6, 7, and 8. The parameter model calculation is not intended to provide a definitive map of k for the SIZ as we still lack sufficient constraint to definitively produce this map. Instead it is intended to emphasize the unique mechanisms that produce ocean surface turbulence, and thus gas exchange, in the IOBL.
Several interesting observations have emerged from this simple analysis, which are worth keeping in mind as we move toward developing an accepted method to calculate k in the polar oceans. The estimates of k from ITP data predominantly occurred at low values of f ; 71 % of the estimates of k were made for f < 0.1. This bias exists because the ITP's are purposefully situated in consolidated sea ice to ensure a long duration of the ITP mission. As ice cover becomes excessively low, the likelihood that the ITP will become dislodged and inoperable increases. This implies that marginal ice zone conditions (with f > 0.5) are under-represented in the data and calculations presented here, both in their likelihood of occurrence and in the gas exchange conditions that exist.
The model indicates that k becomes enhanced during the transition periods of spring and fall. In fall, convection coincides with a nonnegligible fraction of open water area resulting in enhanced gas exchange. In spring, ice melt leads to the formation a shallower mixed layer, which confines the surface stress closer to the air-sea interface, enhancing k. However, a shallower mixed layer is a smaller reservoir for gas storage; while gas transfer may be enhanced in spring, this effect on the net flux of a gas is likely more than offset by equilibration with a much smaller portion of the water column.
The characteristics of sea ice: underside roughness (z o ), floe dimension (L) and fractional ice cover (1 − f ) all play a direct role modulated the magnitude of k. Indeed, (1 − f ) modulates both the stress balance as well as the k eff ; consequently high resolution, high accuracy estimates of these characteristics are valuable to the parameterization of gas transfer.  (2011), including estimates of the gas transfer velocity from 4 the radon deficit method from Fanning and Torres (1991). The stacked bars indicate estimates of 5 keff from sea ice zone effects (sea ice) and from surface wind waves (MSS), during the drift 6 tracks of ITP 3, 6, 7, and 8 computed using equation (2) (2011), including estimates of the gas transfer velocity from the radon deficit method from Fanning and Torres (1991). The stacked bars indicate estimates of k eff from sea ice zone effects (sea ice) and from surface wind waves (MSS), during the drift tracks of ITP 3, 6, 7, and 8 computed using Eq.
It is known from sea ice kinematics that divergent motions cause an ephemeral fraction of the ocean surface to remain exposed, even when satellite sea ice cover is labeled as 100 % (Geiger and Drinkwater, 2005), and some estimates put this ephemeral fraction at nearly 0.1 or 10 % open water . When the fraction of ice cover exceeds 90 %, (f < 0.1), the fraction of open water that we need to resolve is of the same magnitude as the error in satellite ice cover estimates (Knuth and Ackley, 2006), resulting in large uncertainties in the estimate of f . In the absence of a more precise downscaled constraint on f , it would be useful to have some constraint on what the real lower limit is on open water area, if indeed it is not 0 %. The constraints on k and f in the winter when f < 0.1 may be particularly important as they relate to polynyas where much of the ocean's deep water is formed and properties are set.
Finally, the parameter model indicates the effect of shear and convection processes, some of which are unique to the sea ice zone, appear to contribute an additional ca. 40 % to the magnitude of air-sea gas transfer, beyond what would be predicted if a wind speed or mean-square slope parameterization is simply scaled by the fraction of open water, which has been a practice of early attempts to estimate net flux in the sea ice zone (e.g., Stephens and Keeling, 2000;Takahashi et al., 2009).
The parameter model that has been presented here can serve as a means for comparison with laboratory and fieldscale investigations into the connection between ice cover, aqueous turbulence and gas exchange. At a minimum it will be necessary to investigate the effects of shear, convection, stratification and short period wind waves as individual drivers of the effective gas transfer velocity. It is unclear how turbulence and gas exchange in the SIZ transitions to turbulence and gas exchange in the open ocean. However, it is likely that a threshold exists between an upper ocean that is significantly affected by the presence of sea ice, and an ocean that is dominated by direct wind-driven turbulence. When f exceeds some critical value, wind, breaking waves and bubble entrainment are likely to exceed the shear-driven effects of ice floes that are present, but the value of this threshold in ice cover is as yet unknown.

Supplement
A copy of the Matlab routines used to calculate k eff in the sea ice zone has been included as Supplement to this manuscript. A more up to date version of this code may be available at http://geotracerkitchen.org. Navigate to "Recipes".