Modeling Gas Budgets in Marginal Sea Ice Zones

Biogeochemical gas budgets at high-latitude regions and sea ice zones are a source of uncertainty in climate models. The four main processes that regulate these budgets include advection, ventilation, mixing, and accumulation/release from sea ice. Considering the scarcity of data in sea ice zones, specifically during winter time, the environment is too poorly sampled to constrain these processes through direct measurements; hence we proposed models to systematically investigate these processes. The models proposed in this dissertation consist of regional numerical iceocean models, 1D forward and inversion numerical models, and analytical models. Manuscript I of this dissertation focuses on a 3D regional Arctic ice-ocean models. The models are based on MIT general circulation model (MITgcm) code. We used 36 km, 9 km and 2 km horizontal resolution of regional MITgcm configuration with fine vertical spacing to evaluate the capability of the model to reproduce the physical parameters that affect the budget. The model outputs of interest from these simulations are sea ice concentration, sea ice speed, water velocities, and mixed layer depth. From gas budget point of view, sea ice concentration and speed effect ventilation, changes in mixed layer depth lead to mixing, and resolving water velocities quantifies the effects of advection. To assess the accuracy of model, we compared the model outputs to existing field data. We found model sea ice concentration and speed follow data with good fidelity. The model demonstrated the capacity to capture the broad trends in the mixed layer although with a significant bias. We saw improvement in mixed layer depth accuracy with reducing the horizontal resolution of the model. Finally we showed modeled water velocities have low correlation with point-wise in situ data. This correlation remained low in all three model resolution simulations and we argued that is largely due to the quality of the input atmospheric forcing. Manuscript II of this dissertation focuses on 1D forward and inversion modeling of gas budgets. Following our results from the first manuscript, we approximated the effects of advection analytically and utilized a 1D model and its inversion code. We applied the model with combination of numerical passive tracers to reproduce the 53 radon profiles gathered in the Arctic. The optimization based on inversion model reduced the uncertainties in initial conditions and supported the 1D model. We showed mixing, if not resolved, can introduce up to 50% error in estimated budgets. When effects of mixing, melt/freeze and advection taken into the account, we show current estimates of gas exchanges under predicts surface flux in almost cover sea ice areas. Manuscript III presents a new approach in modeling gas exchange in sea ice zones. In this study a sea state dependent gas exchange parametric model is developed based on the turbulent kinetic energy dissipation rate. After comparing this model results with data in the Open Ocean, lakes and marginal ice zones, we applied it to a numerical ice-ocean model of Arctic Ocean. Finally, it is shown that, under the present conditions, gas flux into the Arctic Ocean may be overestimated by 10% if a conventional parameterization is used. In summary, the work presented in this dissertation evaluates and quantifies the effects of environmental forcing on gas budgets in marginal ice zones and offered insight into main factors regulating near surface gas budgets in marginal sea ice zone.


Introduction
The ocean surface is a dynamic region where momentum, heat and salt, as well as biogeochemical compounds are exchanged with the atmosphere and with the deep ocean. At the sea-air interface, gases of biogenic origin and geochemical significance are exchanged with the atmosphere. Theory indicates that the aqueous viscous sublayer, which has a length scale of 20 to 200 μm (Jähne and Haubecker 1998), is the primary bottleneck for air-water exchange. Limitations in measurement at this critical scale have led to approximations of sea-air gas exchange based on indirect measurements. Four approaches involving data are typically used , 1) Parametrization of the turbulent kinetic energy (TKE) at the base of the viscous sublayer 2) Tracing purposefully injected gases Nightingale et al. 2000a) 3) Micro Meteorological methods (H. J. Zemmelink et al. 2006a;Zemmelink et al. 2008;Blomquist et al. 2010;Salter et al. 2011), and 4) Radon-deficit Method.
Here, we examine the radon-deficit method (4), together with a parameterization of the TKE forcing (1) that theoretically leads to the observed deficit in mixed-layer radon.
In the polar oceans wind energy and atmospheric forcing are transferred in a more complex manner as a result of sea ice cover (Loose et al. 2009(Loose et al. , 2014Legge et al. 2015). Sea ice drift due to Ekman flow (McPhee and Martinson 1992), freezing and melting of ice leads on the surface ocean (Morison et al. 1992) and short period waves (Wadhams et al. 1986;Kohout and Meylan 2008) all constitute important sources of momentum transfer. Considering the scarcity of data on marginally covered sea-ice zones (Johnson et al. 2007;Gerdes and KöBerle 2007), especially during Arctic winter time, the environment is too poorly sampled to constrain these processes through direct measurement or empirical relationships.
Lacking sufficient data to constrain these processes, we wonder whether it is possible for a numerical model to adequately capture forcing of air-sea gas exchange in the sea ice zone and consequently improve predictions of air-sea flux. The parameters of interest are sea ice concentration (or fraction of open water), sea ice velocity, mixed layer depth, and water current speed and direction in the ice-ocean boundary layer (IOBL) (Loose et al. 2014). Here we use the budget of 222 Rn gas in the IOBL as an example, because the radon-deficit method has emerged as one of the principle methods to estimate gas exchange velocity in ice-covered waters (Rutgers Van Der Loeff et al. 2014;Loose et al. 2016).
The Radon deficit method involves sampling 222 Rn and 226 Ra in the mixed layer to examine any difference in the concentration or (radio) activity of the two species.
Radon is a gas, radium is a cation; in absence of gas exchange 222 Rn and 226 Ra enter secular equilibrium meaning the amount of 222 Rn produced is equal to decay rate of 5 226 Ra. Any missing 222 Rn in the mixed layer is attributed to exchange with atmosphere (Peng et al. 1979).
Since the 222 Rn concentration in air is very low, less than 5% (Smethie et al. 1985) and considering that concentration is proportional to activity/decay rate A, we can use Eq.
(1) to determine gas exchange. Where k gas transfer velocity in (m d -1 ), A E is the activity or decay rate of 222 Rn which in secular equilibrium is equal to 226 Ra activity, A M is 222 Rn measured decay rate in mixed layer, λ is decay constant of 222 Rn (0.181 d -1 ) and h is the mixed layer depth The mixed layer depth, h, is calculated from the measurements performed at the hydrographic stations during 222 Rn sampling process. Gas transfer velocities from equation (1) reflect the memory of 222 Rn for a period of two to four weeks , which is four to eight times the half-life of 222 Rn (3.8 days).
This memory integrates the physical oceanography properties of the IOBL, including sea ice cover, mixed layer depth and water current speed. These processes are likely to vary significantly during this period and it is important to consider them as a source of uncertainty in Eq. (1). To illustrate this uncertainty, consider a mixed layer that rapidly changes by a factor of 2 just prior to sampling for radon. If the mixed-layer becomes shallower by stratification, h will be smaller by factor of 0.5 while A E /A M in the mixed layer remains the same. Based on equation 1, this causes k to be half of its true value.
That is, prior to stratification TKE forcing was sufficient to ventilate the ocean to a depth greater than the apparent h (Bender et al. 2011). 6 Conversely, if the mixed layer deepens due to mixing, h increases and a new parcel of water with A E /A M = 1 is added to mixed layer, causing the activity ratio to come closer to unity. These two influences on equation 1 (increasing H and A E /A M approaching unity) work against each other, but the net effect is to cause k to appear larger. The change of factor of 2 or higher (in case of convection) in mixed layer depth in less than two weeks has been observed during several studies (Acreman and Jeffery 2007;Ohno et al. 2008;Kara 2003).
The "memory" of gas exchange forcing that radon experiences is further complicated by the presence of sea ice. Consider two alternate water parcel drift paths that lead to the 222 Rn sampling station in sea-ice zone ( Figure 1). Path B demonstrates a history in which water column spends most its back trajectory under sea-ice. Path "A" shows a water column which experiences stratification and shoaling of mixed layer depth equal to δh when drifting through a region that is completely uncovered by ice. During most of Path "B" gas transfer happens in form of diffusion through sea-ice and it will have a very low k (Crabeck et al. 2014;, in contrast Path "A" will have a greater radon deficit, but a smaller h because of stratification. In either case, it is critical to take into account the time history of gas exchange forcing, including changes in the mixed-layer and ice cover, which has led to the apparent radon deficit at the time of measurement.
This observation about drift paths in the sea ice zone strongly implies that we must consider both time and space in estimating the forcing conditions that are recorded in the radon deficit. In other words, we require a Lagrangian back trajectory of water 7 parcels to track the evolution of mixed layer and its relative velocity 4 weeks prior to sampling.
Although satellite data, Ice tethered drifters (Krishfield et al. 2008) and moorings Proshutinsky et al. 2009) have provided valuable seasonal and spatial information about the sea ice zone, they do not track individual water parcels and tend to convolve space and time variations. The spatial limitation of these data poses a challenge to producing a back trajectory of the water parcel.
To address the above mentioned challenges, we use a suite of the Estimation of the Circulation and Climate of the Ocean (ECCO) project's Arctic regional configurations to test the if a numerical model can be used to follow the back trajectory of a radonlabelled water parcel and the gas exchange forcing acting upon it and yield the missing information required for the Radon deficit method.
The variables and derived quantities of interest from the numerical model include mixed layer depth (MLD), sea ice concentration and speed (Loose et al. 2014) and the water velocity in the MLD. We note that as part of the Arctic Ocean Model Intercomparison Project (AOMIP), a number of Arctic ocean-ice models' capability to represent the main ice-ocean dynamics have been assessed (Proshutinsky et al. 2001;Lindsay and Rothrock 1995, p. 995;Proshutinsky et al. 2008). Our reasons for choosing ECCO over other Arctic models stem from the higher correlation between the ECCO's regional Arctic simulated outputs to satellite derived sea ice data (Johnson et al. 2012) and the feasibility in the MITgcm to adapt a high near surface vertical resolution to existing configurations.

8
The remainder of the article is organized as follows: In section 2 we provide the details of the ECCO ice-ocean models. Section 3.1 and 3.2 focus on model outputs of sea-ice concentration and velocity and comparison with observations from satellite and Ice tethered profilers. Section 3.3 investigates the modeled output salinity and temperature structure and the resulting upper ocean density structure and mixed layer.
Section 3.4 evaluates the correlation in near surface water velocity. In section 4 we discuss the results and sources of error and their impact on estimated gas exchange and lastly, section 5 provides the summary of our results.

ECCO model configurations
Three ECCO configurations are used, at horizontal grid spacings of 36 km, 9 km, and 2 km, respectively. The models are based on the Massachusetts Institute of Technology general circulation model (MITgcm) code and employ the z coordinate system described in Adcroft and Campin (2004 (Large et al. 1994) and 36 and 9 km runs utilize salt plume parameterization of Nguyen et. al(2009). The horizontal boundary condition for the 36 km and 9 km configurations comes from existing global ECCO2 model outputs (Marshall et al. 1997;Menemenlis et al. 2008;Losch et al. 2010;Heimbach et al. 2010).
We introduced a set of new vertical grid spacings to allow us to capture near surface small details which cannot be represented with the coarser grid system. In the 36 km (hereafter referred to as A1) and 9 km (called A2) models, the spacing is 2 m in the upper 50 m of the water column and gradually increases to a maximum of 650 meters.
In contrast, the 2 km model (called A3) has 25 layers in the top 100 meters of water column, starting from 1 meter and increasing to 15 meters step. All the boundary conditions from ECCO2 have been interpolated to match the new vertical grid system.

Observations
Satellite-derived estimation of sea-ice cover at 25km horizontal resolution (Comiso 2000) is interpolated to a horizontal grid system to facilitate model-data comparison.

Sea ice concentration
For Sea ice concentration analysis we introduced a grid system covering the Beaufort Gyre and interpolated the data from satellite (Comiso 2008) and A1 on to the grid. The analysis grid extends from 70° to 80° north and 130° to 170° west, covering most of Beaufort Gyre (Figure 3). Grid points can be divided into two main geographic zones that are marked out based on sea ice cover. The first zone contains grid points where the annual average sea ice cover is greater than 80%. These sets of points are fully covered by sea-ice most of the year. The second zone can be described as "marginally ice covered" wherein the ocean surface is free of ice for some fraction of the year. We The RMSE error of 0.2 is the maximum value away from land, this same level of error can also be found near land which is caused by fast-ice generation. Fast Ice in the model is replaced with pack of drifting sea ice; this error is common among numerical models and has been brought to attention during AOMIP (Johnson et al. 2012).
If we compare the monthly climatology for sea ice cover over the 1992-2013 period, the RMS error between model and satellite data is least during the early winter months (e.g. Jan-Mar) when sea ice is close to its maximum extent. Comparing Data and A1 An important source of errors in the model ice concentration comes from the reanalysis surface forcing. Fenty and Heimbach (2012) showed that adjustments in the air temperature that are within the uncertainties of this reanalysis field can help bring the model ice edge into agreement with the observations. Of note also is that the uncertainty in satellite-derived ice cover can be the highest in the marginal ice zone due to tracking algorithms that are sensitive to cloud liquid water or cannot distinguish thin ice from open water (Ivanova et al. 2015), this error also manifests itself in quantification of model-data misfits.

Sea ice velocity
13 Ekman turning causes ice and water to move at divergent angles with respect to each other. Ice moves the fastest, with mean values of 0.09 m s -1 , and the water column progressively winds down in velocity, along the Ekman spiral.
Stratification in the Arctic leads to a confinement of the shear stress closer to the airsea interface and also produces greater divergent flow vectors between ice and water (McPhee 2012). In the marginal ice zone or in regions where ice is converging or diverging, these motions, relative to the motion of the water column can produce significant changes in the water column momentum budget as well as air-sea fluxes.
Thankfully, the ITPs can provide us with a measure of the real ice drift.
To generate a more quantitative comparison between the results we utilized the same in vertical temperature structure compared to actual data. In the profile from ITP-13 ( Figure 5) the model again over estimates the temperature beneath the mixed layer, although certain features including the NSTM can still be found near 10 meters, yet not as pronounced since it is very close to PSW. Bearing in mind that density in the Arctic is dominated by changes in salinity, we move forward to density profiles from this point on.
In addition, we note that recent studies show that eddies with diameters of 30 km or less (Nguyen et al. 2012;Spall et al. 2008;Zhao et al. 2014 We are able to discern some broad similarities between the model and ITP density profiles. From November through January, both ITP and model density profiles remain relatively constant. Between February and March, ITP-35 appears to drift through a zone of convection, likely caused by ice formation, with sudden increase of density near the surface. The same feature can be observed in both A1 and A2 density.
However, on a smaller scale, there is significantly more variation in the ITP data than what the model represents.
For exploring the reason behind the density signals, we used the simulated fraction of sea ice cover and ice thickness ( Figure 6). The dominating effect appears to result from sea ice fraction when there is almost continuously covered area. The changes from sea ice thickness can be observed in the volume of fresh water in the water column, as seen by outcropping of the 1022.5 isopycnal coinciding with the increase of sea ice thickness. An increase in near surface density can be seen in late January and early February accompanied by an increase in ice thickness and insertion of brine in the water column. The second peak, which is not as pronounced, happens in late February when ice fraction decreases from 100% to 95% and exposes the surface water to cold atmosphere, leading to production of newly formed sea ice. We further examine these signals in MLD section below.

Mixed layer depth
There are many different methods in the literature for calculating mixed layer depth We compare these 2 methods by applying them to the profiles from Figure 5 and the results are shown in Figure 7. In case (a) and (b) M1 produces a mixed layer depth that is 8 to 12 meters deeper, compared to the other method. A visual examination of profiles indicates that the M1 criteria may be too flexible of a criteria. The results from M1 appear to be intermittently "realistic", whereas M2 can be difficult to implement for data sampled at high vertical resolution as a result of greater small-scale variability. In practice, we find M1 is the most straight-forward to implement.
It should be mentioned that it is difficult to consistently compare performance of the M1(δρ) and M2 methods on ITP and model data, because the model data extends to top 1 m of water column, whereas the ITP data stops at 7 m depth (Peralta-Ferriz and Woodgate 2015). Furthermore, it has been shown that summer mixed layer in the Canada basin can be less than 12 meters (Toole et al. 2010). To account for this effect,

18
we apply an additional restriction wherein any profile whose mixed-layer depth is less than 2 m below the shallowest ITP measurement is discarded. This restriction effectively removes any ML depths shallower than 10 meters due to ITP sampler not resolving the upper 8 m of water column. In some cases, a remnant mixed-layer from the previous winter may exist in the water column. In this case, the methods incorrectly identify the remnant ML as actual ML depth.
To compare the methods over a longer time period, we calculated the mixed-layer depth from model data and ITP-35 data along the ITP-35 drift track. We used M1 to determine the ML depth for A1, A2 and for ITP-35 data ( Figure 8). Both model results show a shallower ML compared to the ITP data; the most prominent feature in late January corresponds to a sudden change in density found in ( Figure 6). Beside the above mentioned peak A1 fails to capture any variability in MLD whereas A2 shows that ML deepens by about 10 meters in mid February corresponding to ice opening occurring during the same time span ( Figure 6). The difference between A1 and A2 and their ability to capture MLD change, can be explained by the capability of a higher resolution model to capture small-scale fractures in the ice cover ( Figure 8), and conversely, the inability of the coarser resolution to do so is due to averaging over a larger grid. The wind appears to be the primary driving mechanism for the divergence in ice cover, which in turn exposes the ocean to the cold atmosphere and leads to a loss of buoyancy and an increase in MLD. With higher resolution these openings can be captured, leading to a better agreement with data in marginal ice zones. The changes in MLD are of first-order importance to the calculation of gas budgets such as the radon deficit. In this regard a fine-scale grid resolution has real advantages through its ability to capture both the ice advection and openings in ice cover that lead to MLD change. Coarser resolution would be justified when the point of interest is sufficiently far away from leads and marginal ice zones where the effect of sea-ice dynamics on MLD is important, so the effects of area averaging would be small enough to omit.
One last important note is the effect of the salt plume parameterization (SPP) on MLD. Nguyen et al. (2009) demonstrated the need to remove the artificial excessive vertical mixing in coarse horizontal resolution models. To rule out the dependency of this parametrization to vertical resolution as a source in MLD bias, we performed a suite of 1D tests, with and without the SPP on a variety of vertical resolutions (not shown here) and sea ice melting/freezing scenarios and confirmed that SPP is not dependent on vertical grid spacing. We also investigated MLD in A3 (no SPP) run compared to A2, and confirmed the average MLD is the same between these 2 runs.

Velocities in the water column
We have very little information from direct observations that permit us to track a water parcel especially beneath sea ice. This is one area where model output could be critical as there are not obvious alternatives. To assess the consistency of the model water current field, we compared 2D model water velocity to data gathered from two sources: (1)

20
We compared the velocity components averaged from 5 m to 50 m to account for flow direction that is moving the water parcels in the mixed layer over the duration of ITP-V working days which is from Oct 9, 2009 to Mar 31, 2010 ( Figure 9). The ITP data has been daily averaged to remove higher frequency information which we do not expect the model to capture due to the low frequency (6-hourly) wind forcing. Both A1 and A2 show less than 0.3 correlations with data with no improvement in respect to resolution.
We further add A3 to our comparison for moorings velocities (Figure 9), and compared velocities at 25m, which is the level that is shared between all our models and removes the necessity of any interpolation. The simulation results show RMSE normalized by data of higher than 5 and correlations of less than 0.3 over 3 moorings and almost two years of data. This result indicates ocean currents are not well captured in the model irrespective of horizontal grid resolution. We must therefore look into the atmospheric forcing as a likely source of error on high frequency water velocities near surface. As noted above, the wind inputs into the model from the reanalyses are available at a 6-hourly frequency. Chaudhuri et al. (2014) and Lindsay et al. (2014) have compared various available reanalysis products over the Arctic which we used to force our model, along with multiple other reanalysis products with available shipbased and weather station data and found out that wind products in all of those have low correlation i.e less than 0.2. To investigate we compared JRA55 (Onogi et al. sensitivity of k to the IOBL forcing parameters. In the event that we can trust the majority of the model outputs, such as the case here with high fidelity in the simulated SI concentration and SI velocities in A1, we conclude that, a numerical model, even a coarse resolution one, can make significant improvement to the estimate of k. The question of constraining the radon budget within a Lagrangian water parcel is somewhat more complicated.

Application of forcings on Radon budget
The results in Section 3.4 showed that the difference between model and data water trajectories accumulated too much error to be useful, and indicate that for a regional GCM to be useful for reconstructing the back trajectories of radon-labeled water parcels, we will need improved wind-forcing fields. With current reanalysis products, finding the back trajectory of radon-labeled water parcel is not feasible. When Typically sea ice velocity is ~ 5 times greater than vertically averaged water velocity in mixed layer . In this regard, it may be acceptable to assume that the water parcel is stationary as long as ice advection is accounted for. Hence spatial averaging, should account for ice drift over the point of radon/radium sampling. The same logic also applies to the changes in the MLD and sea ice concentration.

Summary
We errors. In addition to sea ice concentration, we also found good correlation (82%) between model ice speed and ITP drift.
The estimation of mixed layer depth is challenging due to its dependence on unconstrained density anomaly or density gradient thresholds. No MLD algorithm performs well in all situations. In addition, CTD profiles from drifting buoys often do not include the top 7-10 m of the surface ocean where stratification can be important.
Adding to the challenge is the dependence of the ocean density structure on vertical fluxes. In these model-data comparisons we found model MLD to be consistently biased on the shallower side in all model resolutions. We note however this result can partly be due to the missing upper 7m in moored drifters such as ITPs and thus resulting in a 1-sided bias in the observed MLD. The evolution of the mixing events showed that MLD correlates to sea ice fraction: in areas of nearly full ice cover, small openings may result in exposure of water to the cold atmosphere and the resulting freezing events would deepen the mixed layer via brine rejection. The higher the 25 resolution, the higher the capability of the model to capture these openings and the resulting deepening effects. The usage of the salt plume parameterization does not play an important role in determining the MLD.
The A1, A2 and A3 experiments consistently could not capture the water velocity observed in ITPs or Mooring. We speculate this discrepancy may be the result of the quality of the reanalysis wind products that are forcing these models. The wind products have been shown to have poor correlation with observed data at high frequencies.
Considering that the response of near surface water is almost instantaneous to the wind forcing, low correlation in wind velocity would have direct impact on the modeled near surface water velocities and likely yield low correlations between modeled and observed ocean currents. On the other hand, the same wind fields at lower frequencies and on broader spatial scale have higher accuracy, as evidenced by the high correlation between the modeled and observed sea ice velocity.
Taking into accounts all the misfits through detailed model-data comparisons, we were able to quantify the usefulness of a numerical model to improve gas exchange rate and parameterization methods. We showed an example of how the sea ice concentration, velocity and mixed-layer depth can affect gas-exchange rate by up to 200% in marginal sea ice zones and that the model outputs can help constrain this rate. By finding the low correlation in near surface ocean velocities, irrespective of model horizontal resolution, we concluded that finding the back trajectory of radon labeled water parcels is currently not feasible. Furthermore, we speculate the source for the common errors in our models, namely the high frequency and under-constrained atmospheric forcing fields, as well as identify alternative approaches to enable the use 26 of a model to achieve the back trajectory calculation task. The alternative approach includes using the MITgcm Green's functions and adjoint capability to help constrain the model ocean velocity to observations, and performing the simulations in a smaller dedicated domain based on the specific spatial distribution of data for both atmospheric winds and ocean currents in the mixed layer.     In the open ocean a series of assumptions including horizontal homogeneity and steady state conditions makes determination and interpretation of the biogeochemical budgets straight forward (Peng et al., 1979;Smethie et al., 1985).
Even though gas exchange velocity, due its dependence on wind (Wanninkhof, 2014), can vary with high frequency compared to budget renewal time, it has been shown that by applying a weighted average on this forcing, a steady state condition can be assumed with <20% error . The weights in this method are based on ventilation ratio of a constant ML depth (Reuer et al., 2007). Hence, the steady state condition is a valid assumption where the ML depth can be assumed near constant. This assumption holds well on time scales shorter than a seasonal cycle where atmospheric forcing gradually changes the ML depth (Acreman & Jeffery, 2007;Kara, 2003;Ohno et al., 2008). A shortcoming in using a forward unconstrained 1D model to capture ML depth changes and its effect on gas budgets is the highly unconstrained ambient conditions, such as initial vertical profile of salinity (S) and temperature (T) (Dwivedi et al., 2011), that dictate the initial ML depth and initial gas inventory. An inversion framework, on the other hand, can be used to obtain a set of initial conditions that are consistent with relevant observation . Here we utilize a 1D model based on MIT General Circulation Model (MITgcm) code and its non-linear inversion (also known as adjoint) framework (Forget et al., 2015b;Wunsch & Heimbach, 2007 to estimate the initial T & S profiles and study the effects of environmental forcing, such as changes in ML depth and sea ice melt/freeze on gas budgets. A suitable approach to investigate the individual contribution of environmental forcing on the gas budgets is to quantify the inventory of near surface radon. Radon is an inert gas, does not participate in biology, and the half-life of radon (3.82 days) allows us to limit the simulations of the physical ocean conditions to 40 days i.e. 5 half-lives of radon.
Sampling of radon budgets is part of a geochemical method called the "Radon deficit method" Peng et al., 1979). This method has been used in the open ocean and sea ice zones to determine the gas transfer velocity Rutgers Van Der Loeff et al., 2014). The exchange transfer velocity inferred from these profiles suggests that, in contrast to other independent measurements Prytherch et al., 2017) The MITgcm adjoint code is generated using automatic differentiation tools (Heimbach et al., 2005) Observational constraints used to invert for initial T/S conditions in this study are the final (day 40) hydrography profile measurements. Day-40 model-data misfits are systematically reduced through iterative minimization of a least-square misfit function (adjoint or Lagrange Multiplier method). At each iteration, incremental adjustments are made to uncertain model parameters and input fields (together termed control variables), which are often not well constrained (Fenty & Heimbach, 2012;Forget, et al., 2015a;Stammer, 2005), to bring the model T/S into agreement with the observed hydrography to within data and model representation errors (Forget, Campin, et al., 2015;. For this study, the control 54 variables are the atmospheric forcing inputs: 10-m air temperature, specific humidity, downward long and short wave, precipitation, and 10-m winds.

2-2 Radon deficit method background
The radon deficit is used to determine the air-sea gas transfer velocity (k). The magnitude of k in sea ice zones has been the focus of scientific community in recent years and this method has emerged as one of the promising tools for determining k in In the surface ocean, the vertically integrated conservation of mass for radon in terms of activities can be written as: where λ is the decay constant of radon (0.181 d -1 ), D is an arbitrary depth, k is gas exchange velocity ,and FmM is freezing minus melting rate. A Rn and A Ra are radon and radium activities. A FmM and A surface denote the activity in sea ice and water 55 surface, finally Adv and Mix accounts for effects of advection and mixing. Based on current geochemical approaches, to account for time dependency of a gas budget, a weighted average of k is used (Reuer et al., 2007). For radioactive budgets, these weights are based on fraction of ML that has been ventilated and radioactive decay rate. This weighted average method allows for neglecting the LHS of (Eq-1). Away from sea ice zones, FmM is equal to zero. Vertical mixing leads to rapid air-sea exchange and the lateral differences in gas concentration are assumed negligible from advection i.e. 2 last terms in RHS are negligible. Rapid mixing would results in vertically constant activities in ML and allows us to set D equal to ML depth.
Therefore the conservation of mass for radon, near the surface can be simplified to: where A Rn observed activity of radon and A ra is observed activity of radium. Both activity values can be obtained by a single vertical sampling of radon and radium.
We compare the 1D

2-3 Implication on other biogeochemical gases such as O 2
The calculations offered in this paper are explicitly done for budget of radon, in this section we are going to demonstrate how these calculations can also be applicable to other geochemical gases. We chose oxygen as an example and offer a scaling argument between similar terms in conservation of mass for radon and oxygen.
The conservation for oxygen can be written as: where is concentration of oxygen, G is gross production of oxygen, R is respiration, is oxygen saturation concentration. In this equation the effects of bubble mediated transport are neglected (Kaiser et al., 2005).
There is a similarity between main sources of uncertainty e.g. advection, mixing, freeze/melt effects in Eq-3 and Eq-1. The possible difference, regarding the uncertainties, between these two equations is the ratio between time scale associated with sources of uncertainty and renewal rate of gas budget in ML (productionconsumption). For gas exchange and radon for example a possible scaling is where k is gas exchange velocity, U 10 is wind at 10 meter above the interface, Sc is the Sea ice formation/melt can result in brine/freshwater discharge in the water column, which in turn would deepen/shoal the mixed layer. Our model captures these phenomena. Sea ice cover, on the other hand, acts as a barrier for gas exchange even if formed outside of our numerical domain and advected to the station's location. We input satellite data (Comiso, 2000) as a mask for gas exchange, and this mask does not affect the heat budget of the water column; hence has no effect on ML depth.

3-1 Optimization
The term t f is the time in days when mixed layer was sampled for radon and radium, and we anticipate that the hydrographic trajectory of a forward 1D model integration from t = t f -40 to t = t f would capture the effect of changes in ML depth on the gas budget, and thus the history of 222 Rn over this period. Continuous changes in ML depth have direct impact on the budget of gas available for gas exchange. Thus an accurate representation of the evolution of the ML depth is critical to our ability to quantify a gas budget. In addition to sea ice conditions over the 40 day period, the ML depth evolution depends on the initial density profile at t = tf -40 and the atmospheric conditions.
In the initial forward 1D model integration, we utilized the outputs from our regional model  to establish the initial sea ice cover, as well as temperature and salinity in the water column. However, this produced hydrography at t = t f that was inconsistent with observed T-S profiles from shipboard CTD casts. In a separate test, we have also initialized sea ice cover and hydrography from satellite ice cover and observed hydrography at stations. However, within the 40-day integration period, the evolution of ML depth still diverged from the observed T-S at t = tf. Based on our previous analysis , we attribute most of this divergence to poor short-term accuracy in the reanalysis atmospheric conditions.
The observed divergence between 1D model forward runs between t = tf-40 and t = tf, revealed that more effort was necessary to adequately capture the history of 222 Rn in the ML. Taking advantage of the adjoint capability in the MITgcm and in the 1-D configuration , we optimized for the initial hydrographic conditions using the observed T-S at t = t f and observed satellite ice concentration as constraints.
At each station, for each observed T-S profile at t = t f , we assigned a corresponding uncertainty. The uncertainty for T (σ t ) and S (σ s ) are vertically constant profiles equal to 0.1 and 0.5 respectively. These uncertainty values are estimated from maximum variability observed between three CTD casts on the sampling stations.
The objective of the optimization is then to estimate an initial T/S profiles (t = t f-40 ) such that hydrographic conditions at t = t f are within the prescribed uncertainty. In other words, the adjoint optimization aims to minimize the misfit function J: In the unconstrained forward runs, the misfit J T and J S as defined above were 53.80 and 22.84, respectively, when averaged over 53 stations. The optimization reduced the average J T and J S to 13.3 and 6.3 through adjustment of uncertain surface atmospheric conditions and initial T and S at t = t f-40 (Forget, et al., 2015a;. For full convergence, J should be reduced to 1, i.e., T-S misfits fall within 62 the uncertainty of the observed profiles. However, due to other unknown factors such uncertainty in our estimated  T ,  S and in the atmospheric forcings, J does not reach full convergence here. However, the improvement in the model T and S at tf is significant (75% and 72% respectively) and demonstrates the effectiveness of the optimization procedure (Figure 2).
Station CB-2a is shown (Figure 2) as an example of this optimization process, where the unconstrained run's misfit J S was significantly higher at 98. After 50 iterations the optimized solution's J S is reduced to 8 (Figure 2-b). The optimization process, by varying the near surface heat budget, can indirectly effect the evolution of sea ice cover (Figure 2-d). Although there are no direct cost (J) related to ice cover, we did not observe more than 30% divergence of ice cover from satellite data. This value is close to 20% uncertainty reported for satellite data (Ivanova et al., 2015).
Again, the only metric in water column available to us is the ML depth of ~6 meters from CTD at time of sampling t = t f . While the unconstrained model shows stratification up to the surface, the optimized ML depth is ~7 meters, which is consistent with the observation. This improvement is crucial in studying near surface radon budgets. The corresponding optimized ML depth evolution is shown in Figure   2-c. Utilizing the optimized ML depth evolution, we proceed to estimate the near surface radon budgets and investigate the effect of the environmental forcings on these budgets.

3-2 Radon budgets and environmental forcings
With atmospheric forcing and initial conditions optimized, we can move to assessment of environmental forcing effects on radon's budget. This task is divided into two steps. First, we statistically compare the simulated budget with available data; next we quantify the contribution of each forcing on the simulated budgets or radon.

3-2-1 Model results and data: gas exchange parametrization in sea ice zones
In this section we compare the modeled radon budgets with the profiles of radon deficit from the 53 sampling stations. As an example, we showed the results of a Even considering the range of Ra across all stations (10 to 18) we cannot explain a phenomenon that is observed in 4 of our summer stations, in which the surface water experiences an increase in concentration of Ra during the melt. One explanation for this phenomenon may be the release of shelf sediments that were entrained into sea ice during the time of ice formation (Notz and Worster, 2009).Further sampling of ice cores may lead to answers regarding why this increase happens instead of the expected dilution. and W14 are bounded by data error bars for ~50% of the stations and we found no statistically significant difference between their results. We showed the result from W14 in Figure 3. The data consistently show higher observed deficits with +80% ice covers, with more than 70% of data landing on positive side of the graph, supporting enhancement of gas exchange beyond the linear scaling in this range of ice covers. Figure 3 Ratio of budgets simulated with W14 formulation and radon data vs ice cover, the red line show 100% accuracy.

3-2-2 Errors introduced by simplifying assumptions in the geochemical budget
As we stated in the introduction, the common convention is to neglect effects of freeze/melt, mixing and assuming a steady state when interpreting mixed-layer budgets of e.g. oxygen and radon. To quantify the error introduced by each assumption, we should compare the corresponding terms neglected in main 66 conservation of mass equation (Eq-1). Discretizing Eq-1 in time, we neglect the advection term for now and we denote the vertical integral of the activity (A) to some depth (D) as the budget (B). Finally, if D is allowed to vary in time we arrive at discrete 1D mass conservation of the budget with no further assumptions: where every term on RHS, except the first term, is a function of time. The budget of radium (B Ra ) can be treated as steady state in the water column, because sedimentary sources are not proximal (Kadko and Muench, 2005) ) predicted by weighted average k Loose et al., 2017;Reuer et al., 2007;Rutgers Van Der Loeff et al., 2014)  . Finally for steady state condition, we apply all these assumptions simultaneously. With these assumption mathematically defined, we move to quantifying the effects of advection.
In some oceanographic circumstances (e.g. in an oligotrophic gyre or far from boundary currents), the assumption of horizontal homogeneity may be more applicable. However in marginal sea ice zones, the lateral gradient of mixed layer depth (Timmermans et al., 2017) and sharp changes in the rate of air-sea flux for heat/momentum/gas due to sea ice cover ) produce significant lateral variations in the surface ocean. The 1D version of the MITgcm does not permit us to capture the effects of advection, but we have developed a simple relationship to capture the effects of lateral heterogeneity. To take into account the horizontal advection of radon we can follow a water parcel that enters our domain from a hypothetical neighboring cell, effectively adopting a Lagrangian framework for that specific parcel.
Following a hypothetical parcel as it enters the numerical grid cell, we make two simplifying assumptions. First, the environmental forcing including wind, k and FmM are the same between the numerical cell and the traveling parcel. Second, we can estimate A as B rn divided by H, where H is the same for the traveling parcel and the numerical cell. In other words, we assumed all deviations from horizontal homogeneity are sea ice related, hence every other forcing can be assumed horizontally constant. With this set of assumptions, Eq10 can be used for traveling parcel as well, although, the specific solution to this equation will be dependent on integration time and initial radon budgets.
When solving Eq11 for the numerical cell, we integrate 40 days forward from secular equilibrium i.e. B = B Ra . In contrast, we have to estimate the initial Rn budget (B ) and time span of integration for the traveling parcel. The integration time can be simply estimated by dividing the length of the numerical cell (4 km) by the magnitude of Ekman transport . Initial budgets (B ) were estimated to be related to ice cover. Ice covers at neighboring cells are obtained by finding the minimum and maximum ice cover from satellite data within a circular region centered at observed station location with a radius of 2 km (i.e horizontal spacing of our numerical cell).
This estimation effectively takes into account the balance between the parcels with high/low radon content coming from high/low ice cover and the time it takes for 69 these parcels to equilibrate compared with the parcels inside the numerical cell. We should note that this method gives an upper limit for effects of advection. By limiting the integration time to the length of the numerical cell divided by Ekman transport, we effectively assumed a straight path for the trajectory of the traveling parcel.
We systemically quantified the error introduced by each of the estimations mentioned above for 53 sampling stations. The errors from No Freeze/Melt, Const ML, Const K, and Steady State are sensitive to ice cover (Figure 4). At the open ocean limit, none of the assumptions result in more than 10% error, the assumption of constant ML depth results in only 5% error. This value is consistent with the 20% error estimated by . These results support the validity of these assumptions in the open ocean limit. The most error from these assumptions, except error from neglecting advection which is only a function of the distance to the ice front, occurs at 10-80% ice cover range.
70 Figure 4 Normlized magnitude of error introduced by neglecting melting/freezing, assuming constant ML, constant k and steady state condtion vs ice cover in % for 53 stations.
We depict the error distribution in the collection of 53 radon deficit stations, using histograms. For example, ignoring the processes of sea ice freeze/melt ( Figure   5-b) leads to less than 10% error on 85% of stations and more than 20% accuracy on only 5% of stations, making this assumption a safe estimate. On the other hand, neglecting advection ( Figure 5-a) has more than 90% accuracy on 80% of stations but can also generate up to 90% error. This kind of error distribution, when left unresolved, would make interpretation of any budget unreliable. The method offered here assists in filtering out the stations that are susceptible to effects of advection. Without the 1D numerical model, there no robust way to estimate the surface value from the integrated budget.

Summary and Conclusion
We utilized an optimized 1D numerical model to explain radon profiles observed during three Arctic cruises (n = 53). The optimization based on inversion model reduced the uncertainties in initial conditions and supported the 1D model in capturing the changes in mixed layer depth.
After resolving the environmental forcing and removing the stations affected by advection, we showed current formulation of gas exchange under predict the radon deficit observed in marginal ice zones.
We offered a simple analytical approach to estimate the errors caused by advection and evaluated the accuracy of current assumptions and estimations of near surface gas budgets in sea ice zones. We showed that assumption of constant mixed layer depth would introduce up to 40% error in estimating near surface gas budgets and importance of resolving this environmental forcing.

Introduction
Constraining the magnitude of the gas exchange velocity (k) at the air-sea interface in the Marginal Ice Zone (MIZ) has implications on estimating the fluxes of biogenic and anthropogenic gases such as carbon dioxide (Bates, 2006;Evans et al., 2015;Liss et al., 2004;Smedsrud et al., 2013;) and methane (Elliott et al., 2011;Uhlig & Loose, 2017;Wåhlström & Meier, 2014). Gas exchange in the MIZ is complex and the relationship between sea ice concentration and carbon 84 sink is poorly understood (Parmentier et al., 2013). Depending on how one describes the sea ice effects on gas exchange, either as natural stirrers acting on Arctic surface water or as barriers to wave generation, its role can vary from enhancer to inhibitor of carbon sink. With the Arctic MIZ extent (Strong & Rigor, 2013) and ice thickness both declining (Lindsay & Schweiger, 2015), understanding the relationship between gas exchange and sea ice cover is crucial, because ice cover alters the kinetics of gas exchange (McGillis et al., 2001).
The MIZ is unique compared to the open ocean (Loose et al., 2009; both with respect to the naturally fetch limited wave field (Smith & Thomson, 2016) and existence of environmental forcings such as sea ice induced turbulence (Loose et al., 2014;McPhee, 2008). Fetch is the distance in which wind acts on the water surface and produces waves. When the energy of the wave is limited by this distance, the resulting wave field is commonly known as fetch limited.
Previous studies on effects of sea ice concentration on gas exchange reveal a field of study that is actively evolving. The first paper to consider the effects of sea ice used a linear dampener on the open ocean gas transfer velocity . Two recent studies support this approximation Prytherch et al., 2017). However, there are other observations that show enhancement (Else et al., 2011;Loose et al., 2017) as well as reduction In the open ocean, results of numerous experiments Ho & Wanninkhof, 2016;Sweeney et al., 2007;Wanninkhof, 1992Wanninkhof, , 2014Wanninkhof et al., 2004) have shown that wind speed can explain more than 80% of variability in gas exchange  There have been theoretical studies suggesting k is related to wave field (Woolf, 2005;Zhao & Toba, 2001;Zhao & Xie, 2010) although such frameworks are not readily applicable in open ocean (Shuiqing & Dongliang, 2016). Here, we take advantage of the previous studies that suggest turbulence regulates gas exchange in the interface (Lamont & Scott, 1970), and specifically the relationship between the Turbulence Kinetic Energy (TKE) dissipation rate and gas exchange (Katul & Liu, 2017;Lorke & Peeters, 2006;Zappa et al., 2007). The TKE dissipation rate has been used to study the effects of various environmental forcings of gas exchange in low wind conditions, such as by tides (Zappa et al., 2003) and rain (Ho et al., 1997;Tokoro et al., 2008;Zappa et al., 2009). It is a challenging task to utilize this method when a combination of wind, waves, and other environmental forcings act on the water column at the same time. We accomplish this task by fitting a TKE based gas Finally, in section 3.4; we apply the parametric model alongside a conventional gas exchange formulation to a regional general circulation model and in section 4 we summarize our findings.

Methods
The model introduced in this manuscript is based on the relationship between the TKE dissipation rate (ε) in water and the gas exchange velocity k (Lamont & Scott, 1970), = ( ) / . ,  where α is a constant, ν is the kinematic viscosity, and Sc is the Schmidt number. Zappa et al. (2007) suggests that α is equal to 0.419 and we adopt this value in this study as well.
When subjected to wind and waves, ε in the water column can vary vertically.
Previous studies suggest that the vertical profile of ε can be divided into 3 regimes.
About 10 times the significant wave height (H s ) below the interface, the law of the wall applies, hence proportionality between ε and z -1 is observed, where z is the distance to the interface. Closer to the interface, in the region between 10 to 2 times H s , the effects of wave breaking shift the profile to z -2 (Terray et al., 1996). Recently, it has been shown that above the z -2 layer, in less than two significant wave heights below the interface, another layer of z -1 proportionality exists (Gemmrich, 2010;Sutherland & Melville, 2015). The latter study shows the relationship between ε and z can be written as, above the depth of about 0.3H s where c p is the wave phase velocity at the wave frequency spectral peak, u *w is the friction velocity in water, U 10 is wind speed at 10 m, and κ is the von Karman constant. The air friction velocity (u *a ) is related to U 10 such that (u *a / U 10 ) 2 is the air sea drag coefficient. The ratio c p /u *a is often called "wave age".
The main challenge in utilizing Equation-2 is determining the depth at which ε is evaluated, hereafter called z * . We achieve this empirically, by matching the value of k determined from a combination of Equation-1 and Equation-2, to the existing empirical quadratic formulations of (Wanninkhof, 1992;Sweeney et al., 2007;Ho and Wanninkhof, 2016). Most of these quadratic formulations are based on open ocean data, such as GasEx experiment (McGillis et al., 2001) in North Atlantic, in a fully developed sea with wave ages of 29 to 32 (Shuiqing & Dongliang, 2016). We set the value of z *, such that the resulting k from Eq-1 and Eq-2 is equal to the quadratic formulation of  (Southern Ocean GasEx) at wave age of 32 which has been observed in same region (Sahlée et al., 2012). Here, we use the COARE 3.5 drag coefficient (Edson et al., 2013) to relate the 10 meter wind speed and the wind friction velocity. The drag coefficient is assumed independent of wave age, since its dependence on wave age is not well understood or constrained (Edson et al., 2013). In Section 3.2 we will show that our main results are not affected significantly even if different sea state dependent drag formulations are introduced instead.
This exercise yields z * values that decrease with wind speed, from 3 m at wind speed 2 m/s to 10 -3 m at wind speed 15 m/s. This rapid decrease of z * is expected because Eq-1 suggests that should be proportional to ( ) to be consistent with the quadratic formulations of , but Eq-2 suggests is proportional to only * if z * is independent of wind speed and the wave age is fixed. This suggests that z * is not related to the actual depth where the gas exchange process is controlled, as long as we assume that k and epsilon are averaged over some long period of time. This is not surprising considering the fact that near surface turbulence is often dominated by intense but infrequent wave breaking events, and that time averaged and time averaged are affected by such events in very different ways. We should therefore interpret z * as a tuning parameter to relate the two well established relationships Eq-1 and Eq-2 and the observed quadratic formulation of .
In order to estimate the k using Equations 1 and 2 in fetch limited conditions, such as in the MIZ, we need to estimate the wave age c p /u *a . The relation between the wave age and the fetch is estimated based on the empirical formulation   where g is gravitational acceleration, x is the fetch. This empirical formula outputs the wave period (T) at the spectral peak. The phase speed c p at the spectral peak is then determined using the deep water dispersion relation. This formulation produces wave 90 ages up to 38 for the open ocean and does not produce older seas (wave age > 40). The relation between fetch, wind speed and wave age is depicted in Figure -a.
We assume that fetch in sea ice zones is related to fraction of sea ice cover by the formulation offered by Smith and Thomson (2016) = 162 where f ice is fraction of ice cover between 0 and 1 and x is effective fetch. No upper limit for this formulation is required, because as f goes to 0, x goes to infinity and the wave period goes to T c .
In addition to the turbulence generated by wind and waves, the near surface turbulence in the MIZ can be further enhanced (diminished) by convection (stratification) and friction around ice floes. These effects are parameterized in the same manner as in Loose et al., (2014). This method contains the effect of buoyant convection/stratification by sea ice freeze/melt (Killawee et al., 1998;Loose et al., 2009) and shear induced by the ice movement (McPhee, 2008).
To assess the spatial and temporal distribution of gas exchange velocity k in realistic settings, the parametric model introduced in this paper, alongside the quadratic formulation offered by Wanninkhof (2014), is applied to a regional Arctic model. The Arctic configuration used in this study, which is based on the MITgcm, is described in .

Open ocean
We compare the WAGT output with conventional quadratic parametrizations (Figure -b). The model is tuned to match the parametrization of  (H06) at the wave age of 32 and it indeed follows that parametrization consistently through all wind regimes. The parametrization of (Wanninkhof, 2014) (W14) is also close to WAGT model at the wave age of 32. The parametrization offered by (Sweeney et al., 2007) (S07) closely follows the WAGT model at the wave age of 34 and  is close to the model at the wave age of 29. All of these wave ages are in the range of fully developed seas.
Typically, wave field information, such as the wave age, is either not recorded or poorly constrained during gas exchange experiments. An exception is the study of Zappa et al. (2007), in which the lower and upper bounds of wave age are reported as of a storm that has been reported during their sampling. Unfortunately, the wave age at each sampling location was not reported, and was not made available at the time this manuscript was prepared. We therefore compare all the data gathered during that experiment (Figure -c)

Fetch limited experiments in lakes
We next evaluate our model in fetch limited condition with data gathered in Siblyback Lake, Dozmary pool (Kwan & Taylor, 1993;Upstill-Goddard et al., 1991), Rockland Lake, Mono Lake, Crowley Lake (Wanninkhof et al., 1987), and Pyramid Lake (Wanninkhof, 1992). Dozmary Pool, Rockland Lake, and Siblyback Lake have less than 2 km maximum length (from shore to shore where t is length of averaging block in seconds. U t and U 3600 are wind speed averaged at t seconds and 3600 seconds. With all the wind data transformed into 1 hour averaged blocks, we compare the observational results with our model outputs. With different fetch limits (lake sizes), the modeled gas exchange velocity starts as a single quadratic shape at lower wind speeds (because the wave field is fully developed and is independent of fetch), and then switches to a nearly linear form when the wave field becomes fetch limited at higher wind speeds. The wind speed where this transition takes place increases as the fetch increases. As the wind speed increases further, there is a second transition point where the wave age becomes 9. At wave age 9, the turbulence enhancement coefficient, that is, 21(c p /U 10 ) 3.5 in Equation-2, reaches unity. It is unlikely that the (averaged) dissipation rate becomes less than the wall layer turbulence value below wave age 9. However, it is unclear whether the gas 95 transfer velocity becomes also independent of sea states below wave age 9. We therefore consider two scenarios; letting the enhancement coefficient fall below one or introducing the lower bound of 1 for enhancement coefficient (that is, setting the wave age minimum of 9 as lower limit). In the second scenario, below wave age 9, ε and consequently k are no longer sensitive to fetch.
The results of the smaller lakes are compared with WAGT model estimates at fetch 1000 m for Siblyback Lake (Figure 2-a) and 500 m for Dozmary Pool and Rockland Lake (Figure 2-b). The gas exchange velocities are significantly lower than the open ocean parameterization (W14) and are quite consistent with our model estimates. The second scenario (with wave age minimum 9) appears to be more consistent with the data.
Crowley Lake, Mono Lake (Figure 3-a), and Pyramid Lake (Figure 3-b) have more than enough size to demonstrate the effect of fetch without the uncertainty of wave age below 9. Again, the reduction of the gas exchange velocity relative to W14 is consistent with the modeled fetch effect. There are other fetch limited experiments e.g. in estuaries and laboratories.
Wave age estimated in laboratory experiments are less than the threshold of 9 (Bock et al., 1999;P. S. Liss & Merlivat, 1986) and they show a significant sensitivity to drag coefficient (Figure not shown here). Estuaries have enough fetch but are subject to tidal forcing (Borges et al., 2004;. Unresolved sources of turbulence and enhancement of k rule out the possibility of comparison with estuary data. As for the lake results shown here the lines of minimum and maximum modeled gas exchange contain most of the data points, while W14 over-estimates the gas transfer velocity reduction that is caused by fetch limitation.

Sea Ice zones
The equation connecting effective fetch and sea ice is offered by Smith and Thomson (2016) during a study on Arctic wave field. This relationship (Eq-6) relates effective fetch with fraction of sea ice to power of -0.49. The wave age goes to 32.2 ( Figure 5-a) as the sea ice fraction goes to zero, but it quickly drops to values below 20 with medium sea ice cover and wind speeds. As U ice increases, the sea ice forcing starts to enhance the k eff . Even with this enhancement the k eff from WAGT model remains less than conventional methods until the sea ice velocity gets close to the free drift velocity (U ice = 0.02 U 10 ). With free drifting ice and sea ice covers below 60%, the WAGT model predicts an enhancement of gas exchange beyond the linear relation.
We first compare our model predictions and the eddy covariance data available from  and (Prytherch et al., 2017). We compare the k eff to k open ratio for the EC data and WAGT parametric model output ( Figure 5 There is a great variability in these two data sets and the error bars cover both the linear relationship and the WAGT model with free drifting ice. Some of this scatter in the EC data may reflect the influence of processes that we attempt to capture with the WAGT model. The WAGT results with 0.01 U 10 < U ice <0.02 U 10 mostly bounds the EC data, suggesting that sea ice drift might account for much of the scattering, however we would require more detail into the ice and water column properties in order to further test the ability of the WAGT model or the linear relationship to reproduce the EC data. In summary, with these data sets both the WAGT and linear model can reproduce the spatially averaged EC data.
We next compare our model with estimates of k by the radon deficit method.
Interpreting the radon data in MIZ is rather a complex task. This complexity is due to variability in mixed layer depth and sea ice forcings that are acting on radon budget prior to sampling. Here we simply employ the weighted time average Loose et al., 2017)  In summary, the radon method results show somewhat better agreement with the WAGT model than with the linear relationship, since the former may explain the observed large variability of k and its dependence on sea ice drift velocity.

Arctic Ocean
We apply the WAGT model alongside the formulation from (Wanninkhof, 2014) (K w14 ) to Arctic regional domain of MIT general circulation model (MITgcm).
The outputs of numerical model, including sea ice cover, sea ice velocity, sea surface temperature and salinity, are daily averaged from 2006 to 2012 and then used as an input to the WAGT model. The quantity of interest is the ratio of effective K w14 to the effective gas exchange velocity from the WAGT model. First, we investigate the time averaged spatial distribution of the ratio ( Figure   6-a). The most pronounced features include overestimation by 20% north of Greenland. With immobile ice cover, such as fast ice or packed ice, the distribution shows a ratio higher than one such is the case in central Arctic. But when the ice is sparse enough to free drift, the ratio is near one as is the east of Greenland and Arctic peripheries. The distribution shows the ratio of one in ice free regions such as Norwegian Seas.
In order to elucidate the temporal variation of the ratio, the daily flux ratio over entire region north of 70° has been calculated and averaged over 5 years for each month ( Figure 6-b), this ratio is independent of partial pressure due to common terms in nominator and denominator. There is an obvious correlation between sea ice cover and the ratio, as the cover increases, the ratio also increases and vice versa, albeit with a lag. This behavior should be viewed from packed ice perspective since areas covered by this type of ice are dominant. As the packed floes start to melt, they start to be more susceptible to wind forcing and gaining higher ice velocity, increasing the denominator of the ratio, hence the ratio minimum in July.
Another aspect of temporal variability is how the ratio would change in response to the area of marginal sea ice zones (MIZ) i.e. areas with more than 1% and less than 99% ice cover from 2006 to 2012 (Figure 6-c). A positive correlation can be found between the ratio and MIZ area. As we move toward the ice free zones and a decrease in sea ice extent, the ratio will get closer to one. Based on the results of the parametric model introduced in this paper, the conventional models may be overestimating the flux into Arctic by 10% in the present climate.
We have developed a parametric Wave Age and Gas Transfer (WAGT) model based on the earlier studies relating the turbulence dissipation rate to gas exchange.
First, the WAGT model is tuned to match the quadratic formulation in open ocean, when the wave field is fully developed. Next, the fetch (or wave age) effect is incorporated into the WAGT model, and the model is successfully tested against observational data in lakes. Sea ice cover is then introduced in the WAGT model both as a source for turbulence and as an effective fetch limiter. The effect of the sea ice velocity has been investigated and it is shown that, on free drift speed, sea ice generates enough turbulence to overcome the effects of wave attenuation.
Finally, the WAGT model is applied to a regional Arctic ocean-sea ice model.
Comparison with the conventional quadratic formulation, with linear relationship with ice cover, shows that the enhancement by sea ice motion is dwarfed by the suppression of wave driven turbulence and consequently, when averaged over 70 degree north, the flux of gases into the Arctic Ocean is overestimated by 10% by conventional formulations in the present climate.
Although our WAGT model is constructed by adding many components, we have so far been able to compare only the combined end result of the gas transfer velocity with existing observations. In the future, it will be desirable to test each model component separately with field and lab data that capture both gas exchange and the forcing, to clarify how different ice related physical processes affect gas exchange.