Numerical investigation of the Arctic ice–ocean boundary layer and implications for air–sea gas fluxes

In ice-covered regions it is challenging to determine constituent budgets – for heat and momentum, but also for biologically and climatically active gases like carbon dioxide and methane. The harsh environment and relative data scarcity make it difficult to characterize even the physical properties of the ocean surface. Here, we sought to evaluate if numerical model output helps us to better estimate the physical forcing that drives the air–sea gas exchange rate (k) in sea ice zones. We used the budget of radioactive 222Rn in the mixed layer to illustrate the effect that sea ice forcing has on gas budgets and air–sea gas exchange. Appropriate constraint of the 222Rn budget requires estimates of sea ice velocity, concentration, mixed-layer depth, and water velocities, as well as their evolution in time and space along the Lagrangian drift track of a mixed-layer water parcel. We used 36, 9 and 2 km horizontal resolution of regional Massachusetts Institute of Technology general circulation model (MITgcm) configuration with fine vertical spacing to evaluate the capability of the model to reproduce these parameters. We then compared the model results to existing field data including satellite, moorings and ice-tethered profilers. We found that mode sea ice coverage agrees with satellite-derived observation 88 to 98 % of the time when averaged over the Beaufort Gyre, and model sea ice speeds have 82 % correlation with observations. The model demonstrated the capacity to capture the broad trends in the mixed layer, although with a significant bias. Model water velocities showed only 29 % 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. Overall, we found that even the coarse-resolution model can make a modest contribution to gas exchange parameterization, by resolving the time variation of parameters that drive the 222Rn budget, including rate of mixed-layer change and sea ice forcings.


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 (Bender et al., 2011): (1) parametrization of the turbulent kinetic energy (TKE) at the base of the viscous sublayer, (2) tracing purposefully injected gases (Ho et al., 2006;Nightingale et al., 2000), (3) micro-meteorological methods (Zemmelink et al., 2006(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.
When the ocean surface is not restricted by fetch, TKE is mostly dominated by wind speed and waves (Wanninkhof, 1992;Zemmelink et al., 2006;Wanninkhof and McGillis, Published by Copernicus Publications on behalf of the European Geosciences Union.1999; Nightingale et al., 2000;Sweeney et al., 2007;Takahashi et al., 2009).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., , 2014;;Legge 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 (MLD), and water current speed and direction in the iceocean 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 icecovered 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 the 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 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 and/or 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 the mixed layer, λ is decay constant of 222 Rn (0.181 d −1 ) and h is the MLD.
The MLD, h, is calculated from the measurements performed at the hydrographic stations during 222 Rn sampling process.Gas transfer velocities from Eq. (1) reflect the memory of 222 Rn for a period of 2 to 4 weeks (Bender et al., 2011), which is 4 to 8 times the half-life of 222 Rn (3.8 days).This memory integrates the physical oceanography properties of the IOBL, including sea ice cover, MLD 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 Eq. ( 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).
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 the mixed layer, causing the activity ratio to come closer to unity.These two influences on Eq. (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 MLD 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 the sea ice zone (Fig. 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 MLD 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;Loose et al., 2011), 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 parcels to track the evolution of the mixed layer and its relative velocity 4 weeks prior to sampling.
Although satellite data, ice-tethered drifters (Krishfield et al., 2008) and moorings (Krishfield et al., 2014;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 radon-labeled 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 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 Massachusetts Institute of Technology general circulation model (MITgcm) to adapt a high near-surface vertical resolution to existing configurations.
The remainder of the article is organized as follows: in Sect. 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 Sect. 4 we discuss the results and sources of error and their impact on estimated gas exchange and lastly, Sect. 5 provides the summary of our results.

ECCO model configurations
Three ECCO configurations are used, at horizontal grid spacings of 36, 9 and 2 km, respectively.The models are based on the MITgcm code and employ the z coordinate system described in Adcroft and Campin (2004).Our approach is to first assess the model outputs from the coarse-resolution model using model-data misfits, then to investigate if there is quantitative reduction in model-data misfits with higher horizontal resolutions.Surface forcings are from the 25-year Japanese Reanalysis Project (JRA-25) (Onogi et al., 2007) for 36 and 9 km runs and the European Centre for Medium-Range Weather Forecasts (ECMWF) analysis for the 2 km run.Initial conditions are from World Ocean Atlas 2005 (Antonov et al., 2006;Locarnini et al., 2006) and initial sea ice conditions are from Zhang and Rothrock (2003) for 36 and 9 km, from which the models are allowed to spin up from 1992.The 2 km global run is initialized from a 4 km spin-up version of the ECCO adjoint-based state estimate for January 2011 and covers the period February 2011 to October 2012.The vertical mixing uses K profile parameterization (KPP) developed by Large et al. (1994) and 36 and 9 km runs utilize salt-plume parameterization (SPP) of Nguyen et al. (2009).The horizontal boundary condition for the 36 and 9 km configurations comes from existing global ECCO 2 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 m.In contrast, the 2 km model (called A3) has 25 layers in the top 100 m of water column, starting from 1 m and increasing to 15 m 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 25 km horizontal resolution (Comiso, 2000) is interpolated to a horizontal grid system to facilitate model-data comparison.In addition, sea ice drift gathered by 28 ice-tethered profilers (ITPs) (Krishfield et al., 2008) which have more than 2 months of data in the Beaufort Sea between 2006 and 2013 have been used to do the ice velocity comparison.
We compared near-surface water velocity data from an icetethered profiler with velocity instruments (ITP-V) (Williams et al., 2010) to A1 and A2 and an upward-looking acoustic doppler current profiler installed on a McLane moored profiler (MMP) (McPhee et al., 2009;Cole et al., 2014) to A1, A2 and A3 in order to compute the accuracy and feasibility of calculating back trajectory of parcels located in the mixed layer.We limit our comparison of ITP-V, which runs from October 2009 to March 2010, to A1 and A2 since those models run from 2006 to 2013 and A3 runs from 2011 to 2013.
Using salinity and temperature profiles from ITPs (Krishfield et al., 2008) we calculated MLD and compared it to 2 m vertical resolution model output (A1, A2).Most of the observed data exist in the Beaufort Gyre, hence we mostly focus our comparison to that geographic perimeter.Figure 2 www.ocean-sci.net/13/61/2017/Ocean Sci., 13, 61-75, 2017 depicts the bathymetry and location of most important observations we used to make the comparisons with the model.

Sea ice concentration
For sea ice concentration analysis we introduced a grid system covering the Beaufort Gyre and interpolated the data from satellite (Comiso, 2000) and A1 onto the grid.The analysis grid extends from 70 to 80 • N and 130 to 170 • W, covering most of Beaufort Gyre (Fig. 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 chose three points within this sea ice geography to compare the seasonal and interannual behavior of the model with satellite ice cover.The points are located at 80 The ice cover at P1, P2 and P3 (Fig. 3) can be divided into 3 ice phases: (a) fully covered in ice, (b) open water and (c) a transition between (a) and (b).P3, which is the furthest south, has all three phases.In contrast P1 ice cover only dips below 60 % for two brief periods during the 7-year time series depicted in Fig. 3 -once in 2008 and again in 2012.These three points illustrate where and when the model has the greatest challenge reproducing the actual sea ice cover.At the extremities of the ice pack, where the water is predominantly covered by 100 or 0 % ice (P1 and P3), the model captures the seasonal advance and retreat and the percentage of ice cover itself is accurate.However, in the transition regions that are characterized by marginal ice for much of the year (P2), the model has more difficulty reproducing the observed sea ice cover as well as the timing of the advance and retreat.This behavior is consistent with the description that has been explained by Johnson et al. (2007), that models have a higher accuracy predicting sea ice concentration in central Arctic and less accuracy near its periphery and the lower latitudes.
The spatial sensitivity of the model can be observed using root mean square (RMS) error (Hyndman and Koehler, 2006) Eq. ( 2), calculated over the 1992-2013 period (Fig. 3).The area with the highest misfits coincides with area between the 80 and 60 % contour lines (Fig. 3) and is concentrated primarily in the Western Beaufort.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 packs of drifting sea ice; this error is common among numerical models and has been brought to attention during AOMIP (Johnson et al., 2012). (2) 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., January-March) when sea ice is close to its maximum extent.Comparing data and A1, Fig. 3 depicts an increase in RMSE during July, August, September and October and a minor decrease in May and November.The RMSE appears to be greater during the summer months of ice retreat, and slightly less during the autumn months of ice advance.Overall, the periods of transition (melt and freeze) coincide with the greatest RMSE.
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 satellitederived 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
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 (Cole et al., 2014), 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 air-sea 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 method introduced by Timmermans et al. ( 2011), to compare ice velocity components (eastward-northward) of A1 to ITP velocity and compute the correlation coefficient of each experiment with the dailyaveraged actual drift velocity from the ITPs (Fig. 4).
When averaged over all the ITPs operating in Beaufort Gyre during 2006 to 2013, A1 had correlations of 0.8 with actual velocity components and 0.82 correlation with speed magnitude.RMSE calculated for A1 based on Eq. ( 2) shows an error of 0.043 ms −1 and no significant bias.

Vertical salinity and temperature profiles
We chose four hydrographic profiles in the Beaufort Sea to assess the simulated vertical salinity and temperature.The first two sets of profiles are from ITP-1 winter and summer 2006; the third set is from ITP-43 during winter 2010 and the fourth is from ITP-13 during summer 2008 (Fig. 5).For visualization we linearly extrapolated the profiles from the first layer of the model up to the surface, which occurs over the top 1 m of the water column.
During winter time, the model temperature and salinity profiles show a well mixed layer that extends below 15 m, followed by a very large gradient.The mixed-layer temperature is close to the local freezing point in a condition called "ice bath" (Shaw et al., 2009) (Shimada et al., 2001, p. 201) or near-surface temperature maximum (NSTM) (Jackson et al., 2010) which is a seasonal feature generated by shortwave solar heat diffusion (Perovich and Maykut, 1990).These two well-defined phenomena are broadly descriptive of the summer surface layer in the Beaufort Gyre.They are, however, absent from the ITP data at this location, indicating a different ice and heat budget time history.Data and model profiles in Fig. 5b show better agreement in the shape and the absolute value of the T and S profiles.Both model and ITP data have a 20 m deep mixed layer during 2010 winter.The model in this case does not show as much change in vertical temperature structure compared to actual data.In the profile from ITP-13 (Fig. 5) the model again over-estimates the temperature beneath the mixed layer, although certain features including the NSTM can still be found near 10 m, 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;Zhao and Timmermans, 2015;Zhao et al., 2016) play an important role in transporting Pacific water from the shelf break into the Canada Basin.Adequate representation of ocean eddies and investigation into their roles in setting the water column stratification require a model with finer horizontal resolution.Hence, moving for-ward, in addition to A1, we utilize the 9 km model (A2) to investigate the density profiles as well as study the MLD.

Density profiles
We compared the 36 and 9 km model outputs of density to the time series of density profiles from ITP-35 (Fig. 6) from October 2009 to March 2010.A black mask indicates locations where there is no data from ITP-35 -particularly in the upper 7 m of the water column.As ITP-35 transited through the Canada Basin, density profiles contain both temporal and spatial changes.
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 a 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 (Fig. 6).The dominating effect appears to result from a sea ice fraction when there is an 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 the production of newly formed sea ice.We further examine these signals in the MLD section below.

Mixed-layer depth
There are many different methods in the literature for calculating MLD (Brainerd and Gregg, 1995;Wijesekera and Gregg, 1996;Thomson and Fine, 2003;de Boyer Montégut et al., 2004;Lorbacher et al., 2006;Shaw et al., 2009).The methods can be divided into two main types (Dong et al., 2008): the first type of algorithm looks for the depth (z MLD ) at which there has been a density increase of δρ between the ocean surface and z MLD .A typical range of values for δρ are 0.005 to 0.125 kg m −3 (Brainerd and Gregg, 1995;de Boyer Montégut et al., 2004).The second type uses slightly different criteria, where the base of the mixed layer is determined as the depth where the gradient of density (∂ρ / ∂z) equals or exceeds a threshold; typical numbers for ∂ρ / ∂z are 0.005 to 0.05 kg m −4 (Brainerd and Gregg, 1995; Lor-  , 2006, p. 200).A more sophisticated approach to type 1 of these criteria is to utilize a differential between ρ 100 m − ρ surface as the cut of point (instead of using a fixed δρ) to account for the effects of surface ρ changes during winter and summer (Shaw et al., 2009).Here, we have implemented two of these methods M1 and M2, with M1 using δρ equal to 0.2 of ρ 100 m − ρ surface (Shaw et al., 2009) and M2 with a gradient (∂ρ / ∂z) cut off point equal to 0.02 kg m −4 , which matches innate model parameterization of MLD (Nguyen et al., 2009).
We compare these two methods by applying them to the profiles from Fig. 5, and the results are shown in Fig. 7.In case (a) and (b) M1 produces a MLD that is 8 to 12 m deeper, compared to the other method.A visual examination of pro-files 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 smallscale variability.In practice, we find M1 is the most straightforward 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 the 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 the summer mixed layer in the Canada Basin can be less than 12 m (Toole et al., 2010).To account for this effect, we apply an additional restriction wherein any profile whose MLD is less than 2 m below the shallowest ITP measurement is discarded.This restriction effectively removes any MLDs shallower than 10 m due to the 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 mixed layer (ML) as actual MLD.
To compare the methods over a longer time period, we calculated the MLD from model data and ITP-35 data along the ITP-35 drift track.We used M1 to determine the MLD for A1, A2 and for ITP-35 data (Fig. 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 Fig. 6.Beside the above-mentioned peak, A1 fails to capture any variability in MLD whereas A2 shows that the ML deepens by about 10 m in mid February corresponding to ice opening occurring during the same time span (Fig. 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 (Fig. 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 firstorder 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 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 parameterization to vertical resolution as a source in MLD bias, we performed a suite of 1-D 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 de- pendent on vertical grid spacing.We also investigated MLD in A3 (no SPP) run compared to A2, and confirmed that the average MLD is the same between these two 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 2-D model water velocity to data gathered from two sources: (1) from ADCPs mounted on moorings that were deployed starting in 2008 in Beaufort Gyre (Proshutinsky et al., 2009) and (2) the ITP-V sensor equipped with MAVSs (Modular Acoustic Velocity Sensors) (Williams et al., 2010), which was the only operating ITP before 2013 which had an acoustic sensor mounted on it.
We compared the velocity components averaged from 5 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 was from 9 October 2009 to 31 March 2010 (Fig. 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 (Fig. 9), and compared velocities at 25 m, 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 three moorings and almost 2 years of data.This result indicates that 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 the 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 ship-based 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 JRA-55 (Onogi et al., 2007) and NCEP (Kalnay et al., 1996) to shipboard data gathered during 2014 in a time span of 2 months in the Arctic and found that JRA-55 had −0.20 correlation, RMSE of 7.36 and bias of −1.3, and NCEP had a correlation of 0.10, RMSE of 5.73 and bias of −1.40 when compared with high-frequency data on each cruise, reinforcing our suspicion of high-frequency wind as a source of error in water currents.

Gas exchange estimation
Up to this point we have spent extensive effort assessing the skill of the MITgcm to reproduce the key forcing parameters listed in our introduction.This effort is motivated by the potential for using the MITgcm model output as a tool to improve our ability to model gas budgets in the IOBL and to improve our estimates of k in the sea ice zone, both of which depend on sea ice processes in the IOBL.To illustrate the potential impact that IOBL properties can have on the estimate of k, we perform a simple experiment, using estimates of k over the range of variation in model output at three locations in the sea ice zone.The intention is to illustrate the variability in k and in the radon deficit that can arise as a result of sea ice processes.

Constraining gas exchange forcings
Utilizing the results from Sect.3.1, 3.2 and 3.3, we calculated gas exchange velocities at P1, P2 and P3 (Fig. 10), over the course of the model simulation (i.e., n = 2557 days × 3 = 7671) introduced in Sect.3.1.The MITwww.ocean-sci.net/13/61/2017/Ocean Sci., 13, 61-75, 2017 gcm IOBL properties are fed to the estimator of k, considering sea ice processes (Loose et al., 2014).Our selected points have the mean sea ice concentrations of 96.1, 87.62 and 61.69 %, sea ice speeds of 0.05, 0.086 and 0.10 ms −1 , wind speeds of 8.73, 5.87 and 4.11 ms −1 .The result yields a point cloud of values that varies depending primarily on the range of ice velocity, wind speed and sea ice cover.The values of k range between 0.1 and 14.0, with a mean of 2.4 and standard deviation of 1.55.This exercise demonstrates the 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 Sect. 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 parcels is not feasible.When improved wind fields are available, the Green's functions approach (Menemenlis et al., 2005;Nguyen et al., 2011) or adjoint method (Forget et al., 2015, p. 4;Wunsch and Heimbach, 2013) can be used to reduce misfits between modeled and observed MLD velocity and likely make the model a valuable tool for tracking back trajectories, either in a smaller domain or full Arctic regional configuration.A possible source of wind data can be from shipboard measurements, assuming the measurements persist over 10 days in the given sampling station.
However, it may be possible to improve on the existing approach.When the drift trajectory is not known, one solution is to resort to averaging IOBL properties within a radius that is equal to the 30-day drift track (e.g., as done by Rutgers Van Der Loeff et al., 2014).The averages within this circle are treated as the representative IOBL properties.The radius of spatial averaging should be restricted by the average magnitude of the water parcel's velocity multiplied by the time span of interest.When applying a spatial averaging, if the timescale of changes in forcings is smaller than time span of interest, the time dependency of forcings should be ac- counted for.Typically sea ice velocity is ∼ 5 times greater than vertically averaged water velocity in the mixed layer (Cole et al., 2014).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 and radium sampling.The same logic also applies to the changes in the MLD and sea ice concentration.For example, gas exchange calculated (Eq. 1) based on assumption of constant MLD of 27.5 m with limits of 5 to 50 m (Peralta-Ferriz and Woodgate 2015) would have limits of ±80 %, whereas gas exchange calculated based on model MLD would have ±50 % error and accounts for time variability.With the current level of uncertainty in reanalysis products and inherent heterogeneity of marginal sea ice zones, we suggest a mixed weighted combination of model outputs and shipboard data to be the way forward for constraining gas budget in sea ice zones.

Summary
We have used 36, 9 and 2 km versions of the ECCO oceansea ice coupled models based on the MITgcm to investigate whether numerical model outputs can be used to compensate for lack of data in constraining air-sea gas exchange rate in the Arctic.The goal is to understand if model outputs can improve estimation of gas exchange velocity calculation and to evaluate the capability of the model to fill in the missing information in the radon deficiency method.This systematic comparison of upper-ocean processes has revealed the following.
The coarse-resolution model showed a good fidelity in regard to reproducing sea ice concentration.Depending on the location and/or season, the error of simulated ice concentration varied between 0.02 and 0.2.Away from ice fronts or active melting/freezing zones the model tended to have higher accuracy.Even in the marginal ice zone, due to the potentially high error in the satellite-derived ice concentration, the model can still be used to quantify the air-sea gas exchange rate, though with an expected higher uncertainty due to the combination of model and data errors.In addition to sea ice concentration, we also found good correlation (82 %) between model ice speed and ITP drift.
The estimation of MLD 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 that this result can partly be due to the missing upper 7 m in moored drifters such as ITPs, thus resulting in a one-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 resolution, the higher the capability of the model to capture these openings and the resulting deepening effects.The usage of the SPP 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 that 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 occurs almost simultaneously to the wind forcing, low correlation in wind velocity would have direct impact on the modeled nearsurface water velocities and likely yield low correlations between modeled and observed ocean currents.Conversely, 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 MLD can affect the 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 as to the source for the common errors in our models, namely the high frequency and under-constrained atmospheric forcing fields, and identify alternative approaches to enable the use 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.

Data availability
Due to sheer volume of A0, A1 and A2 model outputs (89 GB per 3-D field), the simulation results are archived at NASA Supercomputing Center and NCAR GLADE storage system and can be extracted upon request by contacting the primary author.The observational data can be accessed through http://www.whoi.edu/page.do?pid=20781.The postprocessing scripts utilized to compare the observational data and simulation results can be found in the Supplement.
The Supplement related to this article is available online at doi:10.5194/os-13-61-2017-supplement.

Figure 1 .
Figure 1.A graphic illustration of two possible back trajectories for a single sampling station.

64A.Figure 2 .
Figure 2. Bathymetry and location of ITP-V and mooring for data comparison.

Figure 3 .
Figure 3. (a) Averaged satellite sea ice cover from 2006 to 2013, with the solid black line marking 60 % cover and dashed black line marking 80 %; blue dots show the analysis grid, stars show the location of the three points Cyan P1, Green P2 and Red P3 where time-series data is graphed in (b).(b) Time history of sea ice fraction from top P1, P2 and P3, with satellite data represented by blue dots, compared with A1.(c) Horizontal distribution of RMS error of A1 sea ice concentration averaged over time from 2006 to 2013; black mask covers the grid points on the land.(d) Spatially averaged annual RMS error of A1 sea ice concentration.

Figure 6 .
Figure 6.(a) Observed upper-ocean density vs. 36 km (A1) and 9 km (A2) resolution MITgcm density along the path of ITP drift; the black mask covers areas where no ITP data is available and solid black line shows isopycnal of 1022.5 kg m −3 .(b) Simulated sea ice fraction and thickness on top of the water column.

Figure 8 .
Figure8.Sea ice cover higher than 0.9 with gray circle marking the area of ITP operation for (a) 36 km (A1) and (b) 9 km(A2) horizontal resolution of the model.A2 captured the ice opening and resulting mixed-layer change while this phenomenon has been averaged out by a (c) coarse-resolution model observed and simulated evolution of mixed-layer depth on the path of ITP.

Figure 9 .
Figure 9. (a) Daily-averaged velocity components from 5 to 50 m observed by ITP-V vs. those simulated by A1 and A2.(b) Dailyaveraged velocity components at 25 m observed by mooring D vs. A1,A2 and A3.

Figure 10 .
Figure 10.Gas exchange estimated model outputs of wind and sea ice speed at locations P1: 77.4 • N, 143.6 • W, P2: 74.8 • N, 163.5 • W and P3: 70.59 • N, 159.4 • W from January 2006 to December 2012, Areas enclose the outputs around the mean and two standard deviation.The size of the points demonstrate the magnitude of the gas exchange velocities normalized by sea ice cover.

A
. Bigdeli et al.: Numerical investigation of the Arctic ice-ocean boundary layer Figure 4. Time series of sea ice velocity components and speed of ITP 53 vs. 36 km horizontal resolution of MITgcm (A1).The correlations between eastward, northward and magnitude of velocity between ITP 53 data and A1 are 78, 75 and 80 %, respectively.however the ITP MLD is deeper by nearly 10 m, indicating more ice formation and convective heat loss over this water column, as compared to the model water column.In summer the model mixed-layer shoals to approximately 5 m depth following two local temperature extrema; the bigger maximum is at ∼ 35 m, generated by intrusion of the Pacific Summer Water (PSW) which is a dominant feature in the Canada Basin.The second smaller maximum happens around 10 m called the summer mixed layer . The ITP profiles are similar; www.ocean-sci.net/13/61/2017/Ocean Sci., 13, 61-75, 2017