An analysis of air-sea gas exchange for the entire MOSAiC Arctic drift

Arctic sea ice is thinner now than at any time in the past three decades. Every summer the thin ice pack becomes increasingly fractured with leads and summer melt ponds that transmit sunlight and exchange with the ocean (Perovich, 2008; Katlein et al., 2019). These thin ice types all contribute to an increase in insolation to the upper ocean, shortening the spatial and geographic extent of light scarcity and increasing photon availability for primary production. This new radiation budget has created conditions that support phytoplankton blooms underneath sea ice (Mundy et al., 2009; Arrigo et al., 2012) and may even rival the high productivity of the springtime marginal ice zone. Apparently, the Arctic Ocean is becoming more productive, partly due to decreases in sea ice duration and extent, but also as a result of penetration of sunlight through the ice pack and under-ice blooms (Ardyna and Arrigo, 2020; Ardyna et al., 2020). This revelation has transformed the perspective on productivity in the Arctic and the potential for an oligotrophic ecosystem to one that may become more productive over time (Horvat et al., 2017).

While the extent and timing of under-ice blooms of Arctic production are not well-established even less clear is whether the changes in ice type and coverage will increase air-sea gas exchange and the Arctic as an anthropogenic CO₂ sink. In-water primary production decreases the partial pressure of CO₂ (pCO₂) and enhances the pCO₂ differential between surface ocean and atmosphere, favoring ocean CO₂ uptake. However, the surface ocean kinetics that move gas across the interface are easily disrupted by a range of scales of processes, such as increased water surface tension from surface-active compounds (McKenna and McGillis, 2004), stratification (Timmermans et al., 2012), and the imposition of sea ice in its many forms (Anderson et al., 2004; Loose et al., 2017). Beyond the kinetic drivers of gas exchange, the sea ice zone is thermodynamically unique as a consequence of the decoupling of the carbonate system during seawater freezing (Dieckmann et al., 2008; Richaud et al., 2023). Collectively, the multitude of ice types and surface ocean conditions in the sea ice zone mean that gas exchange predictions require detailed knowledge of the sea state and the ice state at a range of scales. As the tools for surface ocean observation become available, the mechanistic and empirical models of air-sea exchange must leverage these tools and be informed by process studies to capture the evolving state of Arctic biogeochemistry.

One outcome for the reduction in sea ice thickness and cover is that the Arctic Ocean will evolve toward a typical open-ocean condition where productivity and gas exchange are closely coupled, and new production is accompanied by an uptake of carbon dioxide from the atmosphere (Stanley et al., 2010) but is also closely compensated by seasonal heterotrophy (Duarte et al., 2013). Another possibility would be that the Arctic Ocean carbon system becomes more asymmetric, where thinner ice allows more radiative transmission to the water column, driving under-ice primary production and lowering the pCO₂, while ice cover and meltwater stratification mitigate air-sea gas exchange. Such conditions of slow air-sea exchange would trap organic matter and respired carbon during the dark months of heterotrophy, so that it can sink or be fixed by primary production, suggesting a more closed-loop Arctic biological pump (Yager et al., 1995). The field and experimental work carried out during the year-long Multidisciplinary drifting Observatory for Arctic Climate (MOSAiC) campaign has provided the opportunity to develop water column budgets of climate-relevant gases over the course of an entire year, but these budgets require a reliable constraint on air-sea fluxes as well as water column processes, like entrainment and water mass exchange.

In this study, we analyzed the rate of air-sea gas exchange, based upon the time-series measurements provided by MOSAiC. Gas exchange was estimated using the Wave Age Gas Transfer (WAGT) model that derives bulk estimates of turbulent kinetic energy (TKE) dissipation and the resulting prediction for air-sea gas transfer velocity under real central Arctic Ocean conditions over the course of an entire year. In Section 1.1 we provide a brief theoretical background on gas exchange, before describing the implementation of the WAGT model during MOSAiC (Section 2.3) and its comparison with in-situ profiles of TKE dissipation (Section 2.2). Section 3 details the WAGT model performance against dissipation profiles, considering the different ice types and forcing conditions. Last, Section 4 discusses the role of ice type and coverage before extending the estimates of transfer velocity over the entire MOSAiC drift year and comparing with a wind speed parameterization of gas exchange.

The flux (F) of a gas across the air-water interface is determined by the gas concentration difference across the interface (C_w − C_a) and by a rate term, that is, the gas transfer velocity (k):

Close to the air-water interface there is a thin layer about 100 µm thick, called the viscous sublayer, where turbulence dies out from surface tension and all processes revert to the slow molecular rate (Jessup et al., 2009). The argument is that the thickness of this layer can be reduced or breached entirely by turbulent eddies, by microbreaking waves, and by bubbles that cross the interface. This reduction or breaching is how the gas transfer kinetics become enhanced and explains why the rate is so variable (Jähne et al., 1987). The most common way to estimate k is based on wind speed (Wanninkhof, 2014), but this captures only 70% of the variability (Woolf, 2005), even when fetch is not limited. Further, the empirical wind speed relationships have been developed for the open ocean and therefore do not reflect an environment where land or ice limit the area of open water and dampen waves (Loose and Schlosser, 2011; Bigdeli et al., 2018). In mixed ice cover, surface turbulence can be affected by wind, currents, ice-water shear, convection, and stratification, local processes that can lead to an observable decrease in the gas transfer velocity compared with a wind speed model (Prytherch and Yelland, 2021). However, a complete model and a robust set of field measurements that highlight these differences is still work in progress.

As turbulence is the ultimate source of enhanced gas transfer, another approach is to use a model that measures the underlying turbulence impinging upon the viscous surface layer of the ocean, because this turbulence is the ultimate source of enhanced gas transfer. The term that best captures all sources of upper ocean turbulence processes and records them along a common scale is ϵ, the viscous rate of dissipation of TKE. The gas transfer velocity scales as (Lamont and Scott, 1970):

where υ is the water viscosity, and Sc is the dimensionless Schmidt number which is the ratio of water viscosity to molecular gas diffusivity. The relationship between ϵ and k is empirical, with much of the evidence that has developed a scaling relationship between these terms coming from laboratory studies like Zappa et al. (2007) and Loose et al. (2016). Because of ice cover it is possible to have high aqueous TKE dissipation (high ϵ) and no air-sea gas flux if there is no open water. However, the gas transfer velocity itself is an explicit representation of the kinetic motion of gas molecules in the viscous surface layer at the water surface, so it is also possible to have high values of k and yet no flux. The production and flux of gases like oxygen, carbon dioxide, dimethyl sulfide, nitrous oxide, and methane in the Arctic are intimately connected to biology at the ice-water interface and biology in sea ice. These gases can accumulate within and beneath the ice or in the brine (e.g., Fransson et al., 2013; Nomura et al., 2018), and high gas concentrations can also occur at the ice surface, in the brine skim and frost flowers (Fransson et al., 2015) exposed to the atmosphere. However, in this work, we focus on ventilation of the water column, which is determined by upper ocean turbulence. Next, we examine how that turbulence drives gas exchange.

The measurements used in this study were collected during the MOSAiC drift, which took place from October 2019 to October 2020 aboard the German icebreaker R/V Polarstern. The vessel was frozen into the Arctic sea ice in the Laptev Sea and drifted across the central Arctic toward Fram Strait in July 2020, eventually exiting into the marginal ice zone, before repositioning further north in pack ice for the final months of the drift. A complete description of the coupled-system concept for MOSAiC and its accompanying experimental program can be found in Nicolaus et al. (2022), Rabe et al. (2022), and Shupe et al. (2022). Crew and scientist exchanges happened approximately every 2.5 months, so that the drift was broken into 5 legs. Measurements were carried out aboard the Polarstern, from the ice floe in what was named the Central Observatory, as well as from a Distributed Network of buoys and other drifting sensors. This study uses time-series measurements as input for a bulk turbulence and gas exchange parameterization. Measurements of wind speed (U_wind) and air temperature (T_air) from the Polarstern meteorological tower, ice drift speed (U_ice) from the GPS position of Polarstern (Schmithüsen, 2021a; 2021b; 2021c; 2021d; 2021e), and water temperature (T_water) and salinity (S_water) collected from the Polarstern sea chest at 11 m depth (Haas et al., 2021; Kanzow et al., 2021; Rex et al., 2021a; 2021b; 2021c) were subsampled to a common hourly time base and merged with daily estimates of sea ice concentration. The meteorological wind speed from the Polarstern is measured at 39 m, so a logarithmic windspeed profile was used to adjust windspeed magnitudes to a 10 m reference. The resulting time series was then used to produce hourly estimates of bulk surface ocean TKE dissipation. Daily sea ice concentration was determined using the reprocessed AMSR2 Bremen sea ice data product with 3.125 km on edge resolution (Spreen et al., 2008). Daily ice concentration was found by identifying the grid cell closest to the position of the Polarstern, and then taking the average ice cover from the surrounding 36 grid cells to yield sea ice cover on a spatial scale of roughly 18.75 km on edge. Details of the bulk turbulence model are described in Section 2.3. Dissipation profiles measured up to the air-sea interface (Section 2.2) were used as an in-situ benchmark against which the WAGT model was evaluated.

The dissipation rate was measured during the MOSAiC drift in the water column using a free-rising microstructure profiler as described by Fer et al. (2022b), following methods used for descending profilers by Fer (2009) and Fer et al. (2012). The profiler was a modified, 500-m-rated, Rockland Scientific International vertical microstructure profiler (VMP-250). Turbulent-shear profiles from the up casts can be used to measure turbulence to a noise floor of less than 5 × 10⁻¹⁰ W kg⁻¹, up to the air-sea interface. Importantly, the creation of an “uprising” mode for the VMP allows the profiler to sample at high resolution up to the ice-water or air-water interface in the critical region where turbulence impinges on the viscous sublayer. The WAGT model creates a single “bulk” value for dissipation and friction velocity; in order to compare those to the vertically resolved VMP data, choosing a depth over which to average the VMP profiles is necessary. This choice was made by iteratively comparing the bulk estimate of ϵ from the WAGT model with the vertical average of ϵ from the VMP profiles. We found that the best fit occurred by averaging dissipation data in the upper 4 m to represent the dissipation in the boundary layer at the surface. Note that the boundary layer thickness was estimated to be 4 ± 2 m using the VMP profiles (Fer et al., 2022b).

The WAGT model estimates ϵ from the simplified bulk budget of TKE assuming that the kinetic energy budget of production and dissipation are always in balance (Gaspar et al., 1990). The primary TKE inputs to the ice-ocean boundary layer are from shear at the ice-water and air-water interface and from convection driven by destabilizing buoyancy fluxes that can arise during surface cooling or sea ice formation. Stabilizing buoyancy fluxes from surface warming or sea ice melting act as a sink in the TKE budget and suppress turbulence. The air-water and ice-water shear are weighted by the fraction of open water (f_ow) and the fraction of ice cover (1 − f_ow; Steele et al., 1989). Last, the wave field and microbreaking (Zappa et al., 2004) are important to TKE production and especially to gas exchange, because wave microbreaking is thought to be a key mechanism that disrupts the viscous surface layer, which is a molecular bottleneck for gas and turbulence. As ice cover increases, fetch for wind-driven shear is reduced and this effect attenuates the wave field (Kohout et al., 2020). The WAGT model was first described in detail by Loose et al. (2014), but the present version, which includes a parameterization of fetch limitation based upon wave age, can be found in Bigdeli et al. (2018). Here, we highlight several details about the WAGT model that are pertinent to its comparison with in-situ estimates of TKE dissipation.

Shear stress at the ice-ocean boundary layer is often the leading order source of momentum and turbulence, especially in an ocean that is almost completely ice-covered. The magnitude of the shear depends on the velocity of the ice and the planetary rotation that leads to Ekman-like flow (McPhee, 2008). This shear profile can be determined using Rossby similarity theory to match the velocity over length scales near the ice-water boundary, captured by Ro_* = u_*0/(fz₀) with the outer layer within dimensionless lengths of u_*0/(fL₀) (McPhee, 1994; Shaw et al., 2008). The surface friction Rossby number, Ro_*, is the ratio of the planetary scale to the under-ice roughness length scale, z₀. The term f is the Coriolis parameter and L₀ is the Obukhov length. Here, we use Rossby similarity to estimate the friction velocity (u_*0) at the ice-water interface, in terms of the ice drift speed V₀, the stratification, and an estimate of z₀:

where A and B are fit parameters determined by the similarity scaling, and Ro_* and μ_* are the inner and outer layer length scales described above. The Obukhov length gives a measure of shear versus convection and can also reflect a stratified or unstable water column (Lombardo and Gregg, 1989). Briefly, the approach is to create an initial estimate of friction velocity u_*0 by matching the velocity scales of the inner and outer planetary boundary layer (Shaw et al., 2008). In practice, the Law of the Wall relation yields an initial estimate of u_*0 that is normalized by the Coriolis parameter times Obukhov length, and by the Coriolis parameter times the roughness length. These upper and lower bounds on the planetary boundary layer are used to develop a universal drag relationship with fit coefficients, A and B, which are used to estimate u_*0 from ice drift velocity. The WAGT model uses the surface ocean temperature and salinity, relative to an estimate of ice temperature and salinity, to determine enthalpy and salt budgets at the ice-water interface and estimate turbulent heat and salt fluxes, as well as changes in water column buoyancy. Altogether the inputs to the WAGT model are U_ice, U_wind, T_ice, T_water, S_water, and f_ow. These terms were estimated from Polarstern underway data, and f_ow was estimated from the AMSR2 3.125 km sea ice concentration as described above.

The atmospheric pCO₂ is available directly from measurements made during the drift (Angot et al., 2022). Seawater pCO₂ was determined by solving the carbonate equation system using total inorganic carbon (DIC) and total alkalinity (TA) measurements that were also made throughout the drift (Ulfsbo et al., 2023). The carbonate system solution was based on the equilibrium coefficients from Sulpis et al. (2020), which are the default settings in the Open Source carbonate system library pyco2sys. To generate a single average pCO₂ differential over the course of the drift, we first averaged over values in the top 50 m at every profile time and then averaged all the profiles in Ulfsbo et al. (2023) over time.

The three primary diagnostic terms used by the WAGT model are depicted in Figure 1 for three individual month-long periods during the MOSAiC drift. Sea ice cover was 95% or greater for all the drift, except for a period in April and from approximately July 15 to July 30, when the Polarstern entered the marginal ice zone; during this period, average ice cover decreased to 86%. During the majority of the drift, the sea ice drift velocity agreed qualitatively well with the relationship U_ice = 0.02 × U_wind (Cole et al., 2014). This relationship holds except in the period when ice cover decreased below 90%, which coincided with a drift velocity that exceeded the 2% relationship: from July 15 to July 30 the average U_ice was 5% of wind speed or 0.24 m s⁻¹. This exception also coincided with the passage of Polarstern through Fram Strait, a narrow exit point from the central Arctic Ocean that is known for high ice drift speeds (Rabe et al., 2022).

There were 235 individual profiles of TKE dissipation in the top 50 m during the MOSAiC drift. Many of these occurred during early spring and summer, when floe sizes were diminished and ice cover ranged from 86% to 100%. To compare estimates between VMP profiles and the bulk surface methods used by the WAGT model, we adopted the same evaluation criteria as used by Fer et al. (2022b). Specifically, an iterative method was used to compute u_*0 from the dissipation profiles against the theoretical relationship $u*03/κ|z|$⁠. Samples with no data in the top 2 m of the water column were excluded, which yielded a total of 167 estimates of u_*0 from dissipation profiles. The values of u_*0 from these VMP profiles ranged from 0.1 cm s⁻¹ to 1.5 cm s⁻¹.

Removing the oscillating inertial and tidal motions from the background “mean” ice drift velocity is often desirable or necessary (McPhee, 1988). The underlying velocity is primarily driven by and thus closely correlated to wind speed. Rossby similarity typically describes the free-drift response, hence the typical use of filtered, background drift velocity, whereas separating the inertial and tidal motion from net drift can provide useful insight into turbulent mixing processes (McPhee, 1988). Initially, we followed the same approach, using the complex-demodulated background ice drift velocity to reproduce boundary layer conditions under the ice and in the Rossby similarity estimates. The complex-demodulated velocity exhibits the best qualitative consistency with the wind speed, reflecting the clear connection between wind speed and the Ekman velocity exhibited by sea ice. However, the VMP profiles exhibited sub-daily variations in the surface values of ϵ_VMP that ranged over nearly two orders of magnitude, variations that were not apparent in either wind speed (Figure 2) or air temperature (not shown). This absence of sub-daily variations led to a reexamination of the pre-processing of inputs to the WAGT model, which determined that the sub-daily oscillations in ice velocity brought on by tide and other periodic motions are an important contributor to turbulence and appeared to couple closely with the in-situ dissipation profiles.

A comparison of u_*0 from Rossby similarity with u_*0 from the VMP is depicted in Figure 3. The instrumentation needed to determine the roughness length (z₀) in situ was not deployed during MOSAiC. Instead, we determined an average value that can be applied to all estimates throughout the drift. To determine the best estimate of z₀ we solved the Rossby similarity solution for a range of 100 values of z₀ from 0.1 mm to 10 cm. At each value of z₀, the values of u_*0 were regressed against the VMP estimates of u_*0 and the value with a slope closest to 1 was chosen. This value coincided with z₀ = 2 mm, which is toward the low end of the range and diagnostic of a relatively smooth under-ice surface farther away from large roughness features such as ice keels (Shaw et al., 2008). The mean of the residual differences between u_*0 from the VMP and u_*0 from Rossby similarity is 0.0014 m s⁻¹, which is equivalent to a 29% misfit between the WAGT and VMP estimates. We considered this agreement to be sufficient to utilize Rossby similarity to reconstruct friction velocity for the duration of the drift.

The 167 VMP profiles with acceptable boundary layer conditions, as chosen by Fer et al. (2022b; see Section 3.1), yielded a mean surface ocean dissipation with upper and lower 95% confidence intervals of $ϵVMP=5.81±[4.10,7.51]×10−7Wkg−1$⁠. The mean dissipation with confidence intervals from the WAGT was $ϵWAGT=5.78±[4.09,7.47]×10−7Wkg−1$⁠, or within 1% of ϵ_VMP while also showing a similar overall dynamic range. Fer et al. (2022b) subdivided the VMP profiles by a classification describing the conditions that the VMP encountered at the surface. During the drift, the VMP reached the surface beneath pack ice, beneath open leads, and beneath leads that were recently covered by thin ice, leading to three categories of profiles: pack, lead, and thin. That classification was used here to further explore where the bulk turbulence formulation by the WAGT model agrees and where it disagrees.

The measured average VMP surface TKE dissipation is nearly identical in leads as compared with under pack ice: $ϵVMP=6.82×10−7$ and 6.81 × 10⁻⁷ W kg⁻¹, respectively. In comparison, dissipation under thin ice was significantly lower, with an average of 2.12 × 10⁻⁷ W kg⁻¹, likely as a result of reduced momentum and buoyancy exchanges with the atmosphere (see Section 4.1 for further discussion). The WAGT model, by contrast, estimates that the average dissipation in leads was 35% larger than under pack ice: $ϵWAGT=7.38×10−7$ versus 5.37 × 10⁻⁷ W kg⁻¹. Like the VMP, the average WAGT model dissipation from instances coincident with the thin ice profiles was lowest: $ϵWAGT=3.64×10−7W kg−1$⁠. While the average dissipation values from WAGT and VMP agree with each other within approximately 25%, how well the WAGT model can predict individual estimates of dissipation, measured on a sub-daily frequency, is worth considering. This degree of predictability can be conveyed using the mean of the normalized residuals (⁠$r^$⁠), where the individual residuals are the absolute difference between the WAGT and VMP dissipation, divided by the mean of VMP dissipation: $r1^=|ϵWAGT−ϵVMP|ϵ^VMP×100,$ expressed in percent.

Using the values of $r^$ under different ice types, the best agreement emerges under pack ice, where $r^$ is 79% (Figure 4). Dissipation from the WAGT model exhibited greater than 100% variation on average for conditions of both open water in leads and thin ice in leads. The poorest agreement emerged under thin ice conditions, where $r^$ was 196% for thin ice (Figure 4). Further, the estimates of ϵ_WAGT in thin ice skew above the 1:1 line, suggesting that Rossby similarity and the WAGT model overestimate TKE dissipation in recently refrozen leads. Below we explore the conditions during which these overestimates may arise.

Surface TKE dissipation under thin ice in refrozen leads is expected to be significantly reduced, compared with the same conditions if that lead contained open water. Thin ice and pancake ice also have a strong mitigating effect on surface waves and sea state (Kohout et al., 2014; Kohout et al., 2020). A similar condition has been studied by Smith and Thomson (2019) while observing momentum transfer in the marginal ice zone, where they estimated that dissipation was 20% or less when compared to the momentum input by the wind. Under these conditions the WAGT model likely overestimates dissipation, primarily because lead refreezing is a small-scale process that is difficult to diagnose without direct on-ice observation and is not classified in the remote sensing data products that are available to the WAGT model. Referring to Figure 4, the WAGT model showed the greatest departure from the VMP estimates under conditions of thin ice (with this departure measured by comparing the average residuals as proportions of the mean of ϵ_VMP). Under thin ice, the WAGT model overestimated dissipation in just over half of the measurements (58%) and also showed the greatest deviation from the 1:1 correspondence with the in-situ VMP dissipation estimates (see dashed mustard line in Figure 4).

Many of the ϵ_VMP profiles under thin ice occurred in late August and September when Polarstern had repositioned near the North Pole and the onset of seasonal cooling was taking place. Air temperature decreased from near 0°C on September 1 to a minimum of −12°C on September 10 (Figure 5). Under those conditions, thin ice formation can produce salt fluxes and convective overturning in the water column (McPhee, 1992). During these conditions, the sub-daily values of ϵ_VMP still ranged over several orders of magnitude, and the more energetic of these dissipation profiles may be associated with salt fluxes.

Perhaps the strongest discrepancy appears when comparing the much narrower dynamic range of variation in ϵ_WAGT, compared with ϵ_VMP. For example, the box plots of daily values of TKE dissipation in Figure 5 show an ϵ_VMP interquartile range (IQR) of more than 10× and as much as 50× for daily values, whereas the IQR for ϵ_WAGT never exceeded 5×. This difference implies that the WAGT model misses the sub-daily intermittent variations in surface turbulence, because the scale and intermittency of thin ice formation is not reflected in the water column properties that were recorded at 11 m by the Polarstern thermosalinograph. In fact, the cooling events at the end of August coincided with an overall decrease in salinity and an increase in temperature measured at 11 m. This example of the complex density structure of the upper water column suggests that thin ice had begun to break down stratification, but that the first outcome was to mix the near surface heat and freshwater down to 11 m or below (Smith et al., 2023). Despite the existence of a salt and buoyancy flux parameterization within the WAGT model, the model is only as detailed as the input parameters. The bulk underway measurements do not capture these intermittent overturning events as a loss of buoyancy from the water column, which means that the WAGT model will not register an increase in TKE dissipation. Here, even if dissipation in the WAGT profiles did increase, this increase would not lead to enhanced gas exchange if the lead was covered with thin ice, because that ice would completely disrupt mass transfer across the ocean viscous layer. This blind spot emphasizes that one of the most valuable tools for diagnosing air-sea exchange is a more detailed view of ocean surface conditions, including the wave field and ice types that may be present at the horizontal scale of 100 m or less.

The degraded agreement between ϵ_WAGT and ϵ_VMP under thin ice conditions highlights the challenge of accurately resolving ocean surface conditions, the presence of sea ice, and even the different ice types. The WAGT model was designed to estimate the area of open water from sea ice concentration data products, such as those provided by MODIS and other passive radiometers. Indeed, tremendous advances have been made in the algorithms used to reconstruct ice types, reducing the on-edge pixel resolution from 12 km on edge to 1 km on edge in just a few years (Spreen et al., 2008). At the same time, the radiometer data treat any ice cover that is 20 cm or thinner as open water (Gunnar Spreen, personal communication, 15/02/2022). This distinction is particularly critical for gas exchange, because gas diffusion through sea ice is many orders of magnitude slower than gas crossing the air-sea interface (Lovely et al., 2015), which means that gas exchange is effectively shut down as soon as even a thin layer of ice cover forms on the surface. The same effect for gas exchange can be expected for pancake and even disaggregated frazil ice crystals. In fact, Matsumura and Ohshima (2015) have shown how seawater with frazil ice crystals acts as a two-phase flow with a corresponding increase in fluid viscosity, and Loose et al. (2023) observed how frazil ice crystals appeared to shutdown air bubble injection—another gas exchange mechanism that arises during capillary-gravity wave microbreaking (Zappa et al., 2004). Recently, there has been some progress in observing these ice types from space-based platforms: Bradtke and Herman (2023) demonstrated how polynya open water area as well as frazil ice streaks can be measured using visible imagery from passive microwave radiometers. Krumpen et al. (2021) also demonstrated how leads and open water can be inferred using sea ice deformation fields, derived from Sentinel-1 SAR imagery. There are still significant restrictions on the timing and coverage of these techniques, and they are not yet at the stage where they can generate a data product reliably; however, they point the way toward better ocean surface characterization than what is possible using sea ice concentration alone.

Converting the dissipation values to gas transfer velocity (k) using Equation 2 makes apparent how restricted gas transfer is during the ice conditions that characterized the MOSAiC drift. The average of the 167 dissipation profiles yielded an average k_VMP of 0.051 m d⁻¹ with an IQR of 0.041 m d⁻¹, compared with the k_WAGT estimate of 0.045 and IQR of 0.043 m d⁻¹. The averages differ by 12% but show a very similar range of variability, and relatively consistent agreement with an R² of 0.82 for regression between k_WAGT and k_VMP (Figure 6). For reference, the global average transfer velocity is 3.5 m d⁻¹, which was computed from the global radiocarbon inventory (Sweeney et al., 2007). The average open water area during the time these profiles were gathered was 3.6%, whereas the mean of k_WAGT is only 1.4% of the global average transfer velocity, revealing that k is lower than would be predicted by a linear scaling with ice cover. Perhaps by coincidence, the average U_wind during the MOSAiC drift was 6.7 m s⁻¹, which is very close to the global average wind speed over the ocean of 6.64 m s⁻¹ (Archer and Jacobson, 2005), rendering the global average transfer velocity a suitable comparison in terms of wind forcing.

Another way to compare the in-situ transfer velocities estimated from dissipation is to use a commonly applied wind speed scaling relationship, the quadratic scaling term from Wanninkhof (2014), hereafter designated as W14: $k=0.06Uwind2$s ⁠, where k is in meters per day and U_wind is the 10 m wind speed in meters per second. Using the values of U_wind and the average ice cover corresponding to the days of the VMP profiles, the W14 scaling predicts k = 0.07 m d⁻¹, which is a 40% overestimate, although admittedly this comparison is between very small values of k.

If we extend the comparison between the W14 wind speed parameterization for gas exchange and the WAGT model values over the entire MOSAiC drift, the discrepancies between the models grow more pronounced. For much of the MOSAiC drift the ice cover was 100% and there was no calculable gas exchange by either method, until March 2020. However, during March and April a period of ice pack divergence created lead openings and an average ice cover around the MOSAiC Central Observatory that diminished to 75% and below. During that period, W14 overpredicts k, compared with the corresponding values from the VMP and WAGT model (Figure 7A and D), and the discrepancy becomes more extreme as open water increases. During the period when the Polarstern drifted close to the marginal ice zone in late June and July, the ice drift speeds increased and sea ice cover decreased below 90%. Under those conditions, the intermittently light winds lead to W14 predictions of k that are lower than those predicted from the WAGT model and the VMP, and gas transfer was enhanced compared to W14. When the Polarstern repositioned to the North Pole in August and September, ice cover was over 90% again and transfer velocities diminished with respect to the prediction from wind speed (Figure 7C and F). Examining the population of gas transfer estimates during the MOSAiC drift, we found that transfer velocity was reduced compared to wind speed prediction for roughly 89% of the drift and enhanced compared to wind speed for 8.6% of the drift, with a nearly perfect match for 2% of the drift.

We used the values of transfer velocity (k_WAGT) to compute the CO₂ flux F between ocean (oc) atmosphere (atm) over the duration of the MOSAiC drift:

The WAGT model outputs k normalized to a Schmidt number of 660, but it can be converted to k_CO2 using the Schmidt number equivalency:

The seawater pCO₂ averaged over the MOSAiC drift was 344.0 ppmv and the average Δ[CO₂] (or [CO₂]_oc − [CO₂]_atm as in Equation 4) between the ocean and the atmosphere over the duration of the MOSAiC drift was −4.9 μmol kg⁻¹ (ΔpCO₂ = −70.32 ppmv); this value was determined using surface layer seawater measurements of total dissolved inorganic carbon (DIC) and alkalinity for calculating the oceanic pCO₂ (Ulfsbo et al., 2023) and atmospheric pCO₂ (Angot et al., 2022) measured throughout the drift. The negative sign on Δ[CO₂] indicates CO₂ uptake by the ocean, a condition that persisted for almost the entire year. Using the average CO₂ difference and the time series of k_CO2 over the entire drift yields maximum fluxes of just higher than −2.5 mmol m⁻² d⁻¹ and a mean value of −0.085 mmol m⁻² d⁻¹. During much of the early drift period from October 2019 to April 2020 the ice cover was so complete that the transfer velocity was zero, and therefore the flux as well, and these values impact the annual flux (Figure 8).

When the time series of CO₂ flux was summed discretely over the entire drift, the total annual CO₂ flux was −0.35 grams carbon per square meter (g C m⁻²; Figure 8). To put this number in context, we can compare this annual flux to estimates of net community productivity (NCP), which is a measure of gross photosynthesis minus gross respiration, from similar regions in the Arctic Ocean. Ulfsbo et al. (2014) found the range of seasonal average NCP to be 10–15 g C m⁻² in the Amundsen and Eurasian basins, with lowest values in the ice-covered central Arctic and values as high as 72 g C m⁻² in the marginal ice zone. A 2016 study using vertical mixing and nitrate fluxes estimated 1–3 g C m⁻² (Randelhoff and Guthrie, 2016), and a benchmark study from 2003 that used a phosphate budget found approximately 0.5 g C m⁻² (Anderson et al., 2003). Collectively, these studies reveal that the production of organic carbon for export (NCP is equivalent to export production) for this region far outweighs the magnitude of CO₂ flux, demonstrating how productivity in the ocean may increase with a thinner, more translucent ice cover, but the drawdown will not necessarily be accompanied by an oceanic uptake of atmospheric CO₂, at least not within the sea ice zone.

The in-situ measurements of turbulent kinetic energy dissipation reveal that the bulk turbulence model for air-sea gas exchange does a reasonable job of capturing the turbulence budget under a range of fetch-limited conditions and is significantly more representative of turbulence conditions when compared to scaling a wind speed parameterization by the ice cover. Considering the relatively few inputs needed to drive the WAGT model (wind speed, ice speed, air temperature, surface water temperature, surface water salinity and ice cover), implementing it is not a computational or observational challenge. The conditions where the WAGT model performed worst in comparison to the in-situ turbulence measurement-based estimates were in leads that had recently refrozen. Under these conditions, WAGT tended to overestimate the turbulence and missed the wide range of variation in turbulence. Because gas exchange is governed by a surface ocean bottleneck, the biggest advance in improving estimates of gas transfer velocity would come from more detailed coverage of thin ice types at a sub-100 m scale. The ability to resolve water surface conditions in a lead would be particularly transformative.

Overall, the mechanisms driving air-sea gas exchange during the MOSAiC drift revealed that gas flux was weak or zero for most of the drift, even during the two periods of ice divergence and low ice cover in April and July. When the transfer velocities were combined with the average difference in CO₂ concentration between the air and the ocean, the annual air-sea CO₂ flux was less than 10% of the annual net community production estimates in this region, suggesting efficient export production that is accompanied by very little atmospheric CO₂ uptake in the ice-covered ocean. However, CO₂ is not thought to be the limiting nutrient on primary production in the ocean, and the low air-sea fluxes are a result of weak vertical turbulence in the water column over the year, suggesting that stratification and weak mixing can also allow organic matter flocculation and sinking that is not slowed by resuspension. This scenario leads to an increasingly closed-loop biological pump, although the factors affecting the strength of the biological pump are broader than just the balance between photosynthesis/respiration and air-sea gas exchange.

The VMP data are available from Fer et al. (2022a); the atmospheric CO₂ data are available from Angot et al. (2022), the in-water CO₂ data are available from Ulfsbo et al. (2023), and the WAGT model can be found at https://github.com/bloose/keff_in_SeaIceZone. The measurements of continuous underway air temperature and wind speed are available by leg (Schmithüsen, 2021a; 2021b; 2021c; 2021d; 2021e), as are the continuous underway measurements of water temperature and salinity (Haas et al., 2021; Kanzow et al., 2021; Rex et al., 2021a; 2021b; 2021c).

Data used in this manuscript were produced as part of the international Multidisciplinary drifting Observatory for the Study of the Arctic Climate (MOSAiC) with the tag MOSAiC20192020 and the Project_ID: AWI_PS122_00. This work was supported by the National Science Foundation Award OPP-1341630 and by the Fulbright US Scholars Program (B. Loose), and by the Research Council of Norway, grant no. 294396 (I. Fer). We thank all those who contributed to MOSAiC and made this endeavor possible (Nixdorf et al., 2021). We would like to thank Brent Else and two anonymous reviewers for their insightful comments which have improved the manuscript.

The authors declare no competing interests.

Concept presented in this article: BL, IF.

Developed MOSAiC carbon system data: AU, MC, ESD, DN, AF, MH, ST-V.

Developed MOSAiC in-situ VMP profiles: IF.

Developed WAGT model: BL.

Wrote the manuscript: BL.

Reviewed and edited the manuscript: BL, IF, AU, MC, ESD, DN, AF, MH, ST-V.

How to cite this article: Loose, B, Fer, I, Ulfsbo, A, Chierici, M, Droste, ES, Nomura, D, Fransson, A, Hoppema, M, Torres-Valdés, S. 2024. An analysis of air-sea gas exchange for the entire MOSAiC Arctic drift. Elementa: Science of the Anthropocene 12(1). DOI: https://doi.org/10.1525/elementa.2023.00128

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Jean-Éric Tremblay, Department of Biology, Université Laval, Quebec, Canada

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: The Multidisciplinary Drifting Observatory for the Study of Arctic Climate (MOSAiC)