Contribution of Sea‐State Dependent Bubbles to Air‐Sea Carbon Dioxide Fluxes

Breaking surface ocean waves produce bubbles that are important for air‐sea gas exchanges, particularly during high winds. In this study we estimate air‐sea CO2 fluxes globally using a new approach that considers the surface wave contribution to gas fluxes. We estimate that 40% of the net air‐sea CO2 flux is via bubbles, with annual, seasonal, and regional variability. When compared to traditional gas‐flux parameterization methods that consider the wind speed alone, we find high‐frequency (daily to weekly) differences in the predicted gas flux using the sea‐state dependent method at spatial scales related to atmospheric weather (10 to 100 km). Seasonal net differences in the air‐sea CO2 flux due to the sea‐state dependence can exceed 20%, with the largest values associated with North Atlantic and North Pacific winter storms. These results confirm that bubbles are important for global gas‐flux dynamics and that sea‐state dependent parameterizations may improve performance of global coupled models.


Geophysical Research Letters
10.1029/2020GL087267 breaking waves (Deike et al., 2016) together with observations and modeling of the wave and wave breaking statistics (see Deike et al., 2017) and validated against field measurements of the gas transfer velocity (Bell et al., 2017;Brumer et al., 2017). In this study we utilize this novel sea-state dependent gas transfer velocity parameterization (see section 2) together with estimates of the global wind and wave fields at three-hourly time increments and ∼50 km (0.5 • ) spatial resolution (see section 3) to investigate the role of bubbles in air-sea CO 2 flux and reveal how variability in the sea-state (due to storm waves, swell waves, etc.) contributes to variability in the flux. Through the application of these methods (sections 4 and 5) we find that accounting for sea-state when computing CO 2 flux changes its value by up to 20% with large-scale coherent spatial patterns that are likely of important consequence for global CO 2 flux modeling. These results demonstrate that considering the sea-state in CO 2 flux parameterizations has important implications for estimating air-sea CO 2 fluxes and for reducing uncertainties in quantifying CO 2 budgets.

Parameterizing the Air-Sea CO 2 Flux
In this section we review the wind only and sea-state dependent (DM18) approaches that are used to parameterize k w in this study.

Wind-Only Approach
Traditional parameterizations for the gas transfer velocity, k w , reflect correlation with wind speed, such as in Wanninkhof (2014, hereafter W14) and Ho et al. (2006): where Sc is the Schmidt number (Sc = ∕D, which relates the fluid viscosity, [m 2 s −1 ], to the gas diffusivity, D [m 2 s −1 ]), U 10 is 10-m wind speed (m s −1 ), and C W14 = 0.251 (cm hr −1 ) (m 2 s −2 ) −1 , yielding k w in conventional units of cm hr −1 . The CO 2 transfer velocity is often given relative to Sc = 660, which is the value of Sc for CO 2 in seawater at 20 • C (k 660 w = k w (Sc∕660) 1∕2 ). The power dependence on U 10 (here quadratic) is an empirical result obtained when considering mean wind speed dependence in laboratory and field studies. The transfer velocity coefficient (here C W14 ) is often constrained using bulk field observations, such as global bomb C14 observations and tracer budget experiments (see the review by Wanninkhof et al., 2009). When the coefficient and power are calibrated using in situ data, the transfer velocity contains the average bubble contribution to gas transfer present during the observation periods. Calibration of the coefficient is common for different wind product applications due to both uncertainty in the large-scale constraints and differences in time-averaging of the mean U 10 .

Sea-State Dependent Approach
In Deike and Melville (2018), the effects of sea-state on k w are included by expressing the gas transfer velocity as the sum of non-bubble k wNB and bubble k wB components (following Keeling, 1993;Woolf, 1997Woolf, , 2005: The bubble-mediated gas transfer velocity k wB in (3) is found by DM18 to vary with significant wave height, H s (m), and wind friction, u * (m s −1 ): where A B is a dimensional fitting coefficient (DM18 find A B = 1 ± 0.2 × 10 −5 s 2 m −2 ), R is the ideal gas constant, T 0 is the sea surface temperature (here skin temperature measured from satellites), and g is gravitational acceleration. Wave breaking tendency (responsible for bubble formation) depends on the high-frequency tail of the wave spectrum and can be formulated with additional characteristics of the wave field, such as the period of the most energetic waves (T p ). While DM18 considered such parameters in determining k wB , the leading order sea-state dependent effects were captured from H s alone in their data (related to an inherent relationship of T p to H s and u * ). The non-bubble gas transfer velocity k wNB in DM18 is found from where A NB is an empirical, nondimensional coefficient given as 1.55 ×10 −4 , corresponding to the the COAREG coefficient A = 1.5 (Fairall et al., 2011). Variations of about 30% for this coefficient can be found in the literature (Bell et al., 2017;Blomquist et al., 2017;Fairall et al., 2011;Smith et al., 2018;Yang et al., 2014).
In this study we apply the DM18 approach to investigate sea-state dependence of air-sea CO 2 flux, but other approaches to consider wave breaking and sea-state dependence exist Woolf, 2005;Zhang, 2012;Zhao et al., 2003). For example, Brumer et al. (2017) find similar mean parameter dependency in k w but expressed using wave Reynolds number (Re H = u * H s ). The wave Reynolds number approach is constrained dimensionally to give equal weight to u * and H s , while the larger exponent for u * than H s in k wB of DM18 arises from scaling of the breaking probability distribution function. While this could indicate enhanced sea-state induced variability in wave Reynolds number based parameterizations, we note that DM18 found generally consistent results for CO 2 between the two approaches. More complex approaches to parameterize bubble-mediated air-sea gas flux for less soluble gases like oxygen and nitrogen also separate the potential terms in parenthesis of equation (1) into bubble and non-bubble components (e.g., Leighton et al., 2018;Zhang, 2012;Liang et al., 2013). Such approaches reduce dependence on seawater solubility and permit super-saturated gas concentration in the ocean. This super-saturation effect is less important for more soluble gases, including carbon dioxide (Liang et al., 2013;Woolf, 1993) and therefore is not investigated here.

Data Description
The bulk formula for air-sea carbon dioxide flux (1) can also be written: where the carbon dioxide concentration difference is expressed as the product of solubility K 0 and air-sea partial pressure difference (Δp CO 2 ) (see Wanninkhof et al., 2009). To investigate global fields of k w and air-sea CO 2 flux we combine several observational and model derived data sets needed to apply equations (3), (4), (5), and (6). We obtain global fields needed to model wind friction (u * ), wave height (H s ), solubility (K 0 ), and carbon concentration difference (Δ p CO 2 ) spanning from 1982 to 2015. These fields are described in this section.

Wind Product
Atmospheric wind velocity is used to estimate wind friction (u * ) and as an input to simulate the wave height (section 3.2). We employ the Japanese Meteorological Society Reanalysis product (JRA55-do), which includes 50 years of 10-m wind vectors (U 10 ) at approximately half-degree resolution every 3 hr (Kobayashi et al., 2015;Tsujino et al., 2018). JRA55-do uses hindcasting, data assimilation, and various statistical constraints to produce a realistic atmospheric state estimate (see Taboada et al., 2019) and was selected by the World Climate Research Program's sixth Coupled Model Intercomparison Project for forcing ocean-ice model experiments (CMIP6, see Griffies et al., 2016). We conducted similar experiments from the CORE atmospheric reanalysis (Large & Yeager, 2009) and found qualitatively similar results.
The wind friction velocity u * is estimated from U 10 using a bulk formula parameterization (Edson et al., 2013). The atmospheric stability is also important in determining u * if there is a sufficiently strong buoyancy gradient. However, it is primarily high wind events that contribute to k wB , a regime where the atmospheric stability effect is typically small; therefore, we neglect the stability effect here.  (Chawla et al., 2013). Modern wave model physics parameterizations are often optimized using observed wave height from buoy and satellite observations and therefore provide similar results for wave height in typical conditions. However, spectral properties such as peak period can be more sensitive to choices in such parameterizations, especially under high wind conditions (e.g., Liu et al., 2017). Furthermore, observational challenges dictate that model performance over large swaths of the ocean remains under-analyzed, though methods to observe high wind and remote locations are emerging (e.g., Thomson et al., 2018).

Wave Height
For this simulation the WW3 spectral resolution is 25 frequencies by 24 directions with spatial resolution equivalent to the JRA55-do forcing grid (half-degree) to fully exploit the relatively high reanalysis resolution. Inputs to the wave model include JRA55-do wind vectors and monthly varying sea-ice from a coupled ocean sea-ice simulation forced with the JRA55-do fields (see Adcroft et al., 2019). Ocean currents induce additional sea-state heterogeneity due to wave refraction and other wave-current interactions but are not included in this study. The wave height statistics from the JRA55-do simulation agree well with Chawla et al. (2013) when considering small differences in the wind (not shown). Since a focus of this study is to reduce uncertainty in gas transfer velocity estimation due to sea-state, we provide the significant wave height from these simulations at doi.org/10.5281/zenodo.3626120.

Carbon Dioxide Partial Pressure ( p CO 2 ) Climatology
The air-sea partial pressure difference of carbon dioxide is not observed globally, but many discrete time series exist (see Rödenbeck et al., 2015, for discussion on various products and interpolation methods). The available oceanic pCO 2 data are synthesized in the Surface Ocean CO 2 Atlas (SOCAT) database Sabine et al., 2012). Landschützer et al. (2014) used a neural-network method to construct global estimates of this quantity at one-degree resolution. Their method is used here to provide monthly Δ p CO 2 spanning from 1982 to 2015 (https://www.nodc.noaa.gov/ocads/oceans/ SPCO_2_1982_2015_ETH_SOM_FFN.html).

Global Field Analysis
Snapshots via punctual model outputs and climatology computed from this study are given in Figure 1. The snapshot fields demonstrate significant spatial and temporal variability associated with atmospheric weather events (50-100 km of horizontal scale) that dominates instantaneous wind, wave, and transfer velocity patterns (Figures 1a-1c). The mean climatology of these fields demonstrates how these events correspond to mean global wind patterns that are reflected in the wave and transfer coefficient fields (Figures 1d-1f). General correlation between wind, wave, and gas transfer patterns is seen, particularly in climatological averages. The wind and wave field patterns are least correlated in the presence of high wind storm events (where sea-state varies along the storm path) and nonlocal swells (where wave heights are not driven locally). These events indicate where the sea-state dependence of the gas transfer velocity k w is expected to be most important for determining the CO 2 flux. These differences are largest in the punctual fields, but key large-scale differences exist in the climatology due to fetch effects associated with prevailing wind direction and land configuration. For example, the eastern North Atlantic and North Pacific show enhanced wave heights associated with the dominant westerly winds.
The bubble contribution to the gas transfer velocity is proportional to k wB ∝ u 5∕3 * gH s 2∕3 while the non-bubble contribution is proportional to k wNB ∝ u * . The bubble contribution to total gas transfer is thus stronger in high wind and wave conditions, physically linked to increasing wave energy and air entrainment. The DM18 bubble contribution to gas transfer velocity, k wB , exceeds the non-bubble contribution, k wNB , for wind speeds greater than 17 m s −1 (see DM18 and section 4). The percentage of time where the wind exceeds 17 m s −1 therefore indicates where bubble-mediated gas transfer is the dominant component of the transfer velocity.  (3)). January (g-i) and July (j-l) climatological mean for the U 10 > 17 m s −1 occurrence (g, j), used as a proxy for where bubbles account for 50% of total gas transfer velocity; the ratio of bubble contribution to total gas transfer velocity k 660 wB ∕k 660 w (h, k), and the net air-sea CO 2 flux (from equation 6) (i, l).
The January and July mean of this threshold occurrence exhibit strong seasonality, with values exceeding 20% in high-latitude winters (Figures 1g and 1j). These regions play a significant role in atmosphere-ocean carbon dioxide exchange, suggesting important bubble induced contributions (e.g., Landschützer et al., 2014;Takahashi et al., 1997). From equation (3), we diagnose the fraction of the gas transfer velocity due to bubbles (k wB ∕k w , Figures 1h and 1k), indicating similar regional and seasonal patterns to those observed on the 17 m s −1 wind speed threshold occurrence.
However, the January and July mean climatology of the net air-sea CO 2 flux (Figures 1i and 1l) suggest that the sea-state dependence of the gas transfer velocity alone cannot completely describe the sea-state effects on CO 2 flux patterns. The sign and magnitude of the carbon concentration differences (Δp CO 2 ) are not correlated to the wave patterns, which suggests that certain regions with high wind and wave activity will contribute more to global net CO 2 fluxes (see further discussion in section 5.2). Thus, the wave contributions to the transfer velocity and the wave contributions to the net CO 2 flux must be considered together.

Global Constraints and Flux
The coefficient in the sea-state formulation (A B in (4)) is determined using field data where the gas transfer velocity is computed from fluxes via eddy covariance methods (Bell et al., 2017;Brumer et al., 2017). The eddy covariance flux methods yield higher transfer velocity at a given wind speed than methods calibrated to bomb C14 or tracer measurements that integrate the flux over longer periods (12 hr to several days) and larger spatial areas, though the reasons for the differences remain poorly understood (Ho et al., 2006(Ho et al., , 2011Wanninkhof et al., 2009;Wanninkhof, 2014). Furthermore, in DM18, all observed data (including wind and wave observations) are averaged into relatively high-frequency bins (20-30 min) compared to the temporal variance captured by the atmospheric forcing (3 hr) used in this study. In order to focus on sea-state added variability, we calibrate the bubble coefficient A B to yield global mean k w values in agreement with the one from Wanninkhof (2014), used in recent global carbon budget estimates (e.g., Landschützer et al., 2014;Le Quére, 2018;Rödenbeck et al., 2015).

Calibrating the Sea-State Formulation Using Global k w Constraints
The average wind speed dependence of the gas transfer velocity from the original sea-state formulation (DM18) is shown in Figure 2a (diamonds) and exceeds the wind-only parameterization (W14) at all wind speeds (thick red line). As discussed in Wanninkhof et al. (2009), even for wind-only gas transfer velocity, the coefficient in equation (2) is routinely adjusted for the wind product to obtain consistent global mean gas transfer velocity. A canonical value for the global mean ⟨k 660 w ⟩ ≈ 16 cm hr −1 (see Ho et al., 2006;Landschützer et al., 2014;Wanninkhof et al., 2009;Wanninkhof, 2014) is often used, but the coefficient remains highly uncertain (Woolf et al., 2019). The DM18 climatological global mean without calibration is ⟨k 660 wDM18,original ⟩ = 19.19 cm hr −1 , which is significantly higher than both the canonical value and the global mean using the wind-only W14 formulation in our study ⟨k 660 wW14 ⟩ = 14.90 cm hr −1 (see Figure 2b, red line). We therefore calibrate the sea-state dependent DM18 formulation (coefficients in equation (4)) such that the global mean gas transfer velocity ⟨k 660 wDM18 ⟩ is consistent with this W14 value. This is done by multiplying a factor of 77.5% times both the bubble and non-bubble components of the gas transfer velocity (Figure 2a, squares) leading to ⟨k 660 wDM18 ⟩ = 14.88 cm hr −1 (Figure 2b, black line). Results presented hereafter include this 77.5% calibration factor. Note that a similar mean global gas transfer velocity could be obtained by making different calibration choices.

Globally Integrated Transfer Velocity and Net CO 2 Fluxes
The calibration described in the previous section dictates that similar global mean gas transfer velocity is found between the sea-state (DM18) and wind-only (W14) formulations by construction (solid lines in Figure 2b). Figure 2c shows the net globally integrated carbon dioxide flux F CO 2 estimated in this study from 1982 to 2015 is also made similar by the calibration. The average value of the net CO 2 flux over the whole period using the calibrated sea-state method is −1.24 ± 0.48 pg C year −1 , in the range of previous studies (e.g., Landschützer et al., 2014).
The average contribution by bubbles to the net CO 2 flux is diagnosed from the sea-state dependent approach as −0.50±0.16 pg C year −1 . The bubble contribution to the gas transfer velocity (k wB ∕k w ) is roughly 30% over the entire time series (Figure 2d, dashed line), in agreement with the estimate of Woolf (1997). The contribution of bubbles to global CO 2 flux from F CO 2 , however, is on average about 40%, with significant monthly to yearly variability (solid line). The relative contribution of bubbles to the CO 2 flux is higher than their contribution to the gas transfer velocity largely because of the sign bias in the flux. This bias arises because in the equatorial regions, CO 2 is degassed to the atmosphere with little bubble contribution. Therefore, when considering the net flux, the net contribution of the degassing region is to increase the percentage of bubble contribution. The interannual variability in the CO 2 flux is not seen in the gas transfer velocity and therefore must instead be primarily related to long-term trends and internal variability in Δ p CO 2 (Rödenbeck et al., 2015).

Influence of Sea-State Dependence on Parameterized Air-Sea CO 2 Flux
The different spatial patterns of the wind and wave fields (see Figure 1) result in differences in the gas transfer velocity when comparing the sea-state dependent and wind-only parameterizations. We now investigate in detail the differences between the gas transfer velocity predicted by the sea-state dependent formulation and the wind-only formulation and the implications for air-sea CO 2 flux.

Seasonal Climatology
At a given time, spatial patterns of large transfer velocities are complex and follow mesoscale atmospheric storm events (10-100 km) characterized by strong winds and waves. The differences of the gas transfer velocity from the sea-state dependent formulation (DM18) and the wind-only formulation (W14) are shown for the climatoligical January and July mean (k DM18 w −k W14 w , Figures 3a and 3d). The patterns in the difference between the two formulations show some prevailing differences due to basin-scale variability in sea-state, such as enhanced values in the eastern North Atlantic. However, the climatological differences also reflect any differences in the shape of the mean gas transfer velocity as a function of wind speed.
The shape difference is explained because the wind-only gas transfer velocity is strictly a quadratic function of wind speed. For the sea-state dependent formulation, the non-bubble gas transfer dominates at low wind and is linear with wind speed, but the bubble gas transfer goes as u 5∕3 * Hs 2∕3 and dominates at high winds. Thus, even though the global mean in k w is similar between the two formulations, the mean value of each at a given wind speed is different and dominates the signal seen in Figures 3a and 3d. To isolate the sea-state effect from the above bias, we introducek DM18 w as the average value of the gas transfer velocity at a given wind speed (e.g., the black circles in Figure 2a) and solve for the sea-state induced difference as Δ w k w = k DM18 w −k DM18 w . The January and July climatological mean sea-state differences, Δ w k w , are shown in Figures 3b and 3e. The patterns show basin-scale patterns of enhanced and suppressed gas transfer velocity that are solely due to prevailing wind and sea-state patterns. The large regions shaded in red (e.g., North Atlantic and North Pacific in January and Southern Ocean year-round) highlight regions where bubble contributions lead to enhanced gas transfer. This sea-state dependent effect reaches up to 10% of the local gas transfer velocity.
The root-mean-square difference between the sea-state dependent formulation (k DM18 w ) and the wind mean (k DM18 w ) indicates the degree of added variability due to the sea-state and is given by (Figures 3c and 3f). The global patterns for √ ∑ (Δ w k w ) 2 are similar to those for Δ w k w , but with values up to three times higher (Figures 3b and 3e). This indicates that induced accumulated Note that these figures all use the 77.5% scaling factor in the DM18 approach. Colored boxes in panels (c) and (f) denote the regions for averaging (North Pacific, North Atlantic, and Southern Ocean) in Figure 4. gas transfer variability due to differences in sea-state is up to 30% in regions of high wave activity in higher latitudes.

Regional Climatology
We further characterize the role of the sea-state in the transfer velocity by diagnosing their seasonality globally and in three regions of high wind-wave activity. The three regions are the North Atlantic (20-80 • N), North Pacific (20-80 • N), and Southern Ocean (20-80 • S), which are marked in Figures 3c and 3f. The monthly averaged gas transfer velocities for the sea-state dependent method (DM18 w/ 77.5% scale factor) and the corresponding net CO 2 flux are given in Figures 4a and 4b. These quantities show significant seasonal variability that differs in magnitude and phase depending on the region. The gas transfer velocities and net fluxes are enhanced in the winter regions due to higher mean winds and increased frequencies of storms. We next investigate the sea-state induced variability to both the gas transfer velocity (defined as spatial averages of √ ∑ (Δ w k w ) 2 in Figure 4c) and net CO 2 flux (defined as spatial averages of √ ∑ (Δ w F CO 2 ) 2 in Figure 4d). The added variability also follows a seasonal cycle in these regions, showing that sea-state induced variability increases in winter due to storm activity.
We observe that the sea-state driven variability of the gas transfer velocity is about 8% globally with small seasonal variation (Figure 4e). Similar numbers are found for the Southern Ocean, with only a small seasonal cycle. However, both the North Pacific and North Atlantic exhibit enhanced variability of the gas transfer velocity to about 10% in winter time, with reduced variability of about 5% in summer. This seasonality in the northern basins again suggests that winter storms play an important role for bubble-driven gas transfer velocity.
We compute similar metrics of seasonality both globally and regionally for the CO 2 flux. We find striking differences in CO 2 flux from gas transfer velocity due to the weighting by the different CO 2 uptake and emission patterns dictated by the sign of the partial pressure Δ p CO 2 . We also note significantly more interannual variability for CO 2 flux than gas transfer velocity (see Figure 4b), which can result in different variability of ) 2 . Panels (e) and (f) show the root-mean-square difference expressed as percentage (middle panels divided by top panels) for (e) gas transfer velocity and (f) net CO 2 flux (net CO 2 flux percentages are masked when the mean is less than 0.01 pg C year −1 ). Markers are colored as blue for northern hemisphere winter (JFM), red for southern hemisphere winter (JAS), and green for other months. The dark black error bars show the 25th-75th percentile range, and the light gray error bars indicate the 5th-95th percentile range.
the sea-state effect (panel d) compared to the gas transfer velocity (panel c). The added variability by considering sea-state when computing CO 2 fluxes is increased to about 10% in the Southern Ocean and to between 10% and 20% in the winter seasons in the North Atlantic and North Pacific. In the summer months, higher fractional variability is observed in the North Atlantic (Figure 4f) that relates to the small regionally averaged flux (Figure 4b). To avoid overstating these fractional variabilities we mask out months where the net flux is approximately zero (August and September). Note that these estimates are limited to an extent because the data set only has monthly means of Δ p CO 2 , and because the gas transfer velocity and the partial pressure are computed independently from one another while in reality they form a coupled system.

Discussion and Conclusion
In this study we apply globally the sea-state dependent parameterization of the gas transfer velocity proposed by Deike and Melville (2018) and investigate how it differs from a traditional wind-only parameterization. The global mean gas transfer velocity from the sea-state dependent formulation is constructed to be similar to the global mean gas transfer velocity of the wind-only formulation that is commonly used when computing global gas exchange budgets. This allows us to quantify the induced variability of the sea-state formulation relative to a wind-only formulation and to estimate the bubble fraction of the gas transfer.
We demonstrate that bubbles contribute about 30% of the gas transfer velocity globally and between 35% and 50% of the globally integrated flux of carbon dioxide, with significant interannual variability due to the variability of the carbon dioxide partial pressure difference. We find that bubbles contribute more than 40% of the net seasonally averaged gas transfer velocity in the Southern Ocean and North Atlantic and Pacific, thereby playing a significant role in air-sea flux of gases such as carbon dioxide in certain regions over seasonal time frames. Spatially there are many tropical regions where the bubble contribution can be significantly less than this (or nearly absent), confirming significant variability to the mechanisms driving air-sea gas flux. The bubble contribution to net CO 2 flux suggests that improved k w parameterizations with variability due to surface waves, for CO 2 but also other gases, may improve coupled model fidelity. We find two primary effects of sea-state dependent parameterization of air-sea carbon dioxide fluxes: 1. Large root-mean-square sea-state variability in both air-sea CO 2 flux and gas transfer velocity that reflect daily and small scale (50-100 km) variability associated with strong winds and transient wave evolution during atmospheric storms (Figures 3c and 3f). This effect indicates that sea-state differences contribute variability up to 30% locally and 10% over large regions in high-latitude winters for high-frequency (weather scale) events that have potential to impact regional fluxes. 2. Coherent basin-scale spatial patterns of sea-state dependence (Figures 3b and 3e) that primarily reflect differences associated with large-scale wind and wave patterns. This indicates that sea-state effects favor gas transfer in regions where waves are more developed due to wind direction and persistence, such as eastern regions of the North Atlantic and North Pacific and the Southern Ocean.
These two primary results indicate that sea-state dependent gas transfer parameterizations can significantly reduce uncertainties in air-sea gas-flux estimation ranging from time and spatial scales associated with atmospheric weather to global and climate scales.
There are limitations of our approach that must be investigated in future work. First, the resolution of wind and wave data used here does not fully resolve weather scale variability and therefore likely underestimates weather scale induced variability in gas transfer. The implications of high-frequency variability for parameterizing transfer velocity and carbon dioxide flux in coupled models is not yet understood. The accumulated flux variability in this study is likely amplified due to computing the air-sea carbon dioxide flux using global reconstruction of the partial pressure difference Δ p CO 2 . It therefore remains unclear how the sea-state dependence of the gas transfer velocity will effect coupled simulations. To investigate such implications, studies resolving the necessary high-frequency temporal resolution could be achieved by specifically designed field campaigns and through sensitivity tests in coupled simulations. The effects on ocean interior CO 2 concentrations in models are particularly worth investigating, since regions with high variability of the CO 2 flux due to sea-state in this study include the Southern Ocean and the North Atlantic, both regions of deep-water production and CO 2 export (e.g., Sabine et al., 2004). Importantly the largest bubble impacts occur in winter, when ocean mixed layer depths are deepest, which could lead to enhanced model ocean interior CO 2 uptake through nonlinear rectification effects related to winter mode water ventilation. Finally, we note that for this experiment we used the bulk parameterization from Deike and Melville (2018) that only accounts for the significant wave height and significantly reduces the computational requirement. This approach may underestimate the variability due to sea-state in the fluxes versus using the high-resolution spectral representation of the wave and wave breaking statistics. This parameterization is based on a more complete formulation that accurately represents high-frequency (full spectral) complexity of the wave and wave breaking fields. These points will be investigated in the future using recent progress in modeling the breaking statistics and the wave spectrum (Romero, 2019).