A Sea State Dependent Gas Transfer Velocity for CO2 Unifying Theory, Model, and Field Data

Wave breaking induced bubbles contribute a significant part of air‐sea gas fluxes. Recent modeling of the sea state dependent CO2 flux found that bubbles contribute up to ∼40% of the total CO2 air‐sea fluxes (Reichl & Deike, 2020, https://doi.org/10.1029/2020gl087267). In this study, we implement the sea state dependent bubble gas transfer formulation of Deike and Melville (2018, https://doi.org/10.1029/2018gl078758) into a spectral wave model (WAVEWATCH III) incorporating the spectral modeling of the wave breaking distribution from Romero (2019, https://doi.org/10.1029/2019gl083408). We evaluate the accuracy of the sea state dependent gas transfer parameterization against available measurements of CO2 gas transfer velocity from 9 data sets (11 research cruises, see Yang et al. (2022, https://doi.org/10.3389/fmars.2022.826421)). The sea state dependent parameterization for CO2 gas transfer velocity is consistent with observations, while the traditional wind‐only parameterization used in most global models slightly underestimates the observations of gas transfer velocity. We produce a climatology of the sea state dependent gas transfer velocity using reanalysis wind and wave data spanning 1980–2017. The climatology shows that the enhanced gas transfer velocity occurs frequently in regions with developed sea states (with strong wave breaking and high significant wave height). The present study provides a general sea state dependent parameterization for gas transfer, which can be implemented in global coupled models.


Introduction
The air-sea flux of gases, such as carbon dioxide (CO 2 ), oxygen (O 2 ), nitrogen (N 2 ) and dimethyl sulfide (DMS), is a key part of the climate system and a major factor in ocean biogeochemical processes (Wanninkhof, 2014).Significant efforts have been made in developing air-sea gas exchange parameterizations to constrain air-sea gas fluxes in a variety of wind conditions, especially for CO 2 (D'Asaro & McNeil, 2007;Deike & Melville, 2018;Ho et al., 2006;McNeil & D'Asaro, 2007;Wanninkhof, 1992Wanninkhof, , 2011Wanninkhof, , 2014;;Woolf et al., 2019).Air sea fluxes of CO 2 are an important component of the global carbon cycle with ocean uptake accounting for approximately 30% of the carbon dioxide released into the atmosphere (Friedlingstein et al., 2022;Resplandy et al., 2018;Woolf et al., 2019).
The air-sea gas flux is obtained from the gas transfer velocity (k w , mks unit:m/s) and air-sea gradient of gas concentration, which can be expressed in terms of the gas solubility (K o , unit:mol/m 3 /atm) and partial pressure (  Δ =    −    , with unit atm, where     is the partial pressure of gas in the water and     in the air near the water surface): Deike and Melville (2018, hereafter DM18) proposed a physics-based framework to model the bubble mediated transfer velocity for intermediate to high solubility (i.e., CO 2 and dimethylsufide (DMS)), using the chemical model from Keeling (1993) together with spectral statistics of breaking waves from Phillips (1985) to describe air bubble entrainment.The gas transfer formulation was validated for CO 2 and DMS against field measurements with wind speeds up to 25 m/s and reconciles recent data sets from two North Atlantic field campaigns (Bell et al., 2017;Brumer, Zappa, Blomquist, et al., 2017).Reichl and Deike (2020, hereafter RD20) implemented a simplified sea state dependent formulation from DM18 into a global wave model and found that the bubble mediated gas transfer contributes about 40% of the global CO 2 flux.Separately, Romero (2019) developed a wave breaking statistic model linking the wave spectrum and the breaking distribution and implemented it into the spectral wave model WAVEWATCH III (WW3, WW3DG ( 2019)).
The combined work from Deike and Melville (2018); Reichl and Deike (2020); Romero (2019) are utilized in the present study.We investigate the sea state dependent gas transfer velocity parameterizations developed by Deike and Melville (2018), including both forms based on the full wavenumber-direction wave spectrum (hereafter the spectral model) and based on only the significant wave height (hereafter the bulk model).The spectral model further requires estimating the wave breaking probability distribution function, which leverages the wave breaking model developed by Romero (2019).We also discuss the classic wind only gas transfer velocity parameterization (Wanninkhof, 2014), such as typically employed in present coupled Earth System Models (e.g., Dunne et al., 2020).All three formulations are systematically compared to an extensive data set from 11 different research campaigns for air-sea CO 2 exchange, as summarized by Yang et al. (2022).Leveraging the comprehensive data set from Yang et al. (2022), we optimize the sea state dependent flux model coefficients to obtain good agreement between the sea state dependent parameterizations and the observations over the full set of environmental conditions.We discuss when differences emerge between the two sea state dependent formulations, particularly in the presence of strong swell at low to moderate wind speed, and note remaining uncertainties at high wind speed (above 20-25 m/s) related to the limited observational data sets.Furthermore, we apply the resulting optimized spectral model for gas transfer velocity for CO 2 from this work to multi-decadal global wave simulations to investigate climatological patterns in the gas transfer velocity, extending the similar investigation using the bulk model in Reichl and Deike (2020).We investigate the climatological role of bubbles in gas transfer velocity for CO 2 and the difference between the three (spectral, bulk, and wind-only) parameterizations in terms of global and seasonal patterns.The results indicate that bubbles play a key role in CO 2 gas transfer velocity in the North Atlantic, North Pacific and Southern Oceans, which is consistent with the analysis of Reichl and Deike (2020).In general, the bulk sea state dependent gas transfer velocity formula prediction is consistent with the spectral model results when U 10 ≤ 20 m/s.Our climatology analysis reveals that the difference between the mean gas transfer velocity from both bulk and spectral formulas is small when U 10 < 10 m/s.The difference between the spectral and bulk mean gas transfer velocity increases to 15 cm/ hr as U 10 increases to 30 m/s, which is a difference of ∼12.5% of the estimated gas transfer velocity.The maximum seasonal mean difference of sea state induced variability from bulk and spectral formula is around 3 cm/hr (∼2.5%).The differences suggest that the bulk parameterization approximates the spectral gas transfer velocity fairy well in most conditions for general use.However, this difference has important implications for complex ZHOU ET AL.

10.1029/2023EA003237
3 of 18 sea states, such as can occur in strong swell in low wind conditions (<10 m/s) to obtain a more skillful gas flux prediction.

Gas Transfer Velocity Model
In this section, we review and describe the three parameterizations for gas transfer velocity used in this study.First we introduce the wind-speed only parameterization from Wanninkhof (2014); Ho et al. (2011) (Section 2.1), second we discuss the bubble mediated sea state dependent model based on the wave spectrum from Deike and Melville (2018) (Section 2.2.1), and finally we present the sea-state dependent parameterization based on the bulk wave properties (Deike & Melville, 2018) (Section 2.2.2).

Wind-Only Approach
The wind only parameterization for gas transfer velocity commonly used in large-scale ocean and climate models (e.g., Dunne et al., 2020) is only dependent on the neutral, 10 m wind speed, such as Wanninkhof (2014, hereafter wan14): here k w,wan14 is typically given with unit cm/hr, such that C wan14 = 0.25 (cm/hr)/(m/s) 2 ,   =   is the Schmidt number, D and ν are the CO 2 molecular diffusivity and water kinematic viscosity (m 2 /s) that depend on salinity and temperature and are computed from empirical formulations, see Wanninkhof (1992).

Bubble Mediated Gas Transfer Velocity
To consider explicitly the role of bubbles, we separate the gas transfer velocity into two components: the bubble and non-bubble gas transfer velocity (Deike, 2022;Deike & Melville, 2018;Keeling, 1993;Woolf & Thorpe, 1991): (3) The non-bubble component is due to molecular diffusion enhanced by turbulence and is written as: where A NB is a non-dimensional coefficient obtained from empirical data on gas transfer and u * is the wind friction velocity.The scaling in Sc number and wind speed can be justified by eddy renewal theory (Garbe et al., 2014).Uncertainty of A NB is reported around 30% from previous literature (Bell et al., 2017;Fairall et al., 2003;Goddijn-Murphy et al., 2016;Hare et al., 2004;Smith et al., 2018;Yang et al., 2014).

Spectral Model
The bubble-mediated gas transfer velocity accounts for the volume of air entrained by wave breaking V a obtained from the statistics of breaking waves in terms of the length of breaking crest distribution Λ(c) introduced by Phillips (1985).The efficiency of the bubbles to transfer gas is then obtained from the chemical model from Keeling (1993) as described in Deike and Melville (2018).The final formula for the bubble mediated transfer reads: where α = K o RT is the Ostwald dimensionless solubility coefficient, K o (mol/m 3 /atm) is the gas solubility in sea water, R is the ideal gas constant (J/K/mol), and T the sea surface skin temperature with unit K.  B is a constant representing the balance between mechanical dissipation and the work done by buoyancy (Deike & Melville, 2018;Deike et al., 2016), with values between  B = 0.1 ∼ 0.2 that is discussed in Section 4.2.2.The wave slope is expressed as , where B T = 1.1 × 10 −3 is a constant representing the breaking threshold following Romero et al. (2012) and B(k) = ϕ(k)k 3 is the saturation spectrum.Here, ϕ(k) is the directionally integrated wave spectrum and k is the wave number.χ is a nondimensional constant of order 1 to ensure a closed energy dissipation budget (Romero, 2019).The breaking wave distribution Λ(c) is computed from the wave spectrum ϕ(k) following Romero (2019) and Shin et al. (2022).The bubble gas exchange efficiency factor ZHOU ET AL. 10.1029/2023EA003237 4 of 18 E(r) is defined as a function of the depth of bubble injection z 0 and an equilibration depth H eq (r) together with the gas diffusivity and solubility.The equilibrium depth H eq is a function of the bubble rise velocity w b (r) and the exchange velocity for a single bubble κ b (r), which depends on bubble radius r, solubility and gas diffusivity of gas in sea water.Here we use the bubble rise velocity w b (r) in quiescent water following Keeling (1993) and Deike and Melville (2018).The equations are , with ∶ (6a) where z 0 is the bubble injection depth, z 0 = H s /2 in this study (Deike & Melville, 2018),  = 10.82∕ χ , and    =  3  2 .The bubble size distribution n(r) of entrained bubbles, defined as the number of bubbles per unit radius per unit volume, follows experimental, computational and theoretical work on individual breaking waves and bubble break-up (Deane & Stokes, 2002;Deike, 2022;Deike et al., 2016;Mostert et al., 2022;Rivière et al., 2022), We note that the volume flux of air entrained by breaking waves V a is the rate of entrained air volume per unit area of ocean surface (unit: m/s), which is a function of wave breaking crest length distribution (Λ(c)) Deike and Melville (2018); Deike (2022), that is calculated in the spectral wave model as we describe in Section 3: which constrains the total volume of air entrained, so that V a = ∫n(r)4/3πr 3 dr.
The details of all parameters are described in Deike and Melville (2018).This spectral sea-state dependent gas transfer formulation requires both the wave spectrum ϕ(k) and the distribution of length of breaking crest Λ(c) both function of wind speed and wave history.

Bulk Model
The wave spectrum and breaking crest distribution are not quantities that are always available during field measurements as they require specific instrumentations (Kleiss & Melville, 2010;Sutherland & Melville, 2013).A simplified sea state dependent gas flux parameterization can therefore provide an approach that is dependent on more commonly observed bulk sea state properties to evaluate air-sea gas flux, such as would be measured on saildrones Nickford et al. (2022).A simplified sea state dependent parameterization of k b is proposed by Deike and Melville (2018) based on the full spectral model described above which constrain the dependency in the wind friction velocity and the significant wave height: This bulk relationship based on significant wave height (H s ) can reproduce the observational data from HiWinGS campaign (Brumer, Zappa, Brooks, et al., 2017) and Knorr 2011 campaign (Bell et al., 2017).Deike and Melville (2018) suggested the best fit coefficients A b = 1.0 ± 0.2 × 10 −5 s 2 /m 2 .The exponent in terms of solubility and Sc number are here updated compared to Deike and Melville (2018); Reichl and Deike (2020) to correspond to the one suggested by the full bubble model from Keeling (1993).In this work we optimized the parameters in this approach based on the available 11 research cruises (see Figure 1).
ZHOU ET AL. 10.1029/2023EA003237 5 of 18 We also note that we will compute the wind frictional velocity u * from 10-m wind speed (U 10 ) using the COARE3.5 drag coefficient (Edson et al., 2013).The drag coefficient itself could also be parameterized as a function of the sea state (making the non-breaking transfer velocity also a sea-state dependent term), but we choose to not consider this effect here to simplify our analysis and due to remaining uncertainty in the sea state dependence of drag coefficients (e.g., Sauvage et al., 2023;Zhou et al., 2022), in particular in high wind speed conditions, or mixed seas conditions.

Summarizing the Three Gas Transfer Velocity Parameterizations in This Study
Results for the gas transfer velocity will be presented in terms of the gas transfer velocity normalized to Schmidt number of 660 which accounts for variations in diffusivity due to temperature, ) 1∕2 (hereafter, wind-only gas transfer velocity).We summarize the primary symbols used in this study in Appendix A.

WAVEWATCH III Model Setup
In this study, the WAVEWATCH III (WW3) spectral wave model v5.16 (WW3DG, 2019) is used to simulate the wave spectrum in a global configuration.The WW3 spectral resolution is 42 frequencies by 36 directions with spatial resolution equivalent to 0.5°.The resolved wavenumber range is chosen to span from 0.0016 rad/m  2018)), which has approximately half-degree resolution every 3 hr.Inputs to the wave model also include monthly varying sea-ice from a coupled ocean sea-ice simulation forced with the JRA55-do fields (see Adcroft et al. (2019).
In this work, we also need estimates of the wave breaking statistics.Romero (2019) developed a framework for spectral wave models to compute Λ(c) and V a in Section 2.2.1 from a single framework by linking the wave dissipation source term with wave breaking statistics.Note that (Romero, 2019)  , where M W is a wind modulation function, B br is a breaking threshold, and M L is a longer waves modulation to breaking.All of above terms can be calculated from the directional wave spectrum F (k, θ), as described in Romero (2019).In Romero (2019); Shin et al. (2022), the lower limit of integration of Λ(c) in Equations 8 and 5 is 2 m/s.In this study, we assume Λ(c) is a constant when c ≤ 2 m/s, the lower limit of the integral through c is 1.57m/s corresponding to the resolved highest wavenumber in the WW3 model set-up employed in this work.

Observational Data
Yang et al. ( 2022) summarizes hourly gas transfer velocity of CO 2 by using eddy covariance technique from 11 research cruises (9 data sets) (Figure 1).The full data set includes field campaigns in the North Atlantic Ocean, Southern Ocean, tropical Indian Ocean and Arctic Ocean in different seasons and years as shown in Figure 1a (Bell et al., 2017;Blomquist et al., 2017;Butterworth & Miller, 2016;Landwehr et al., 2018;Miller et al., 2009;Yang et al., 2021Yang et al., , 2022;;Zavarsky et al., 2018).Their study estimates a bulk u * with the COARE3.5 model (Edson et al., 2013) from the shipboard wind speed measurement.We evaluate the JRA55-do wind forcing based on the wind observations in all above cruises (Figure 1b).The root mean square error (RMSE) is 0.82 m/s between JRA55-do U 10 and observational U 10 , which suggests the wind forcing is sufficiently accurate for the wave simulations to recreate the wave fields from these ship campaigns.Figure 1c shows the comparison between the frictional velocity u * calculated from JRA55-do wind forcing and the shipboard wind speed measurements in Yang et al. (2022).Since there is no observations of the significant wave height H s during the 11 cruises, Yang et al. (2022) used the H s from ECMWF product for the 11 research cruises (Janssen & Bidlot, 2021).Figure 1d shows the WW3 simulated significant wave height (H s ) against the H s from ECMWF ecWAM CY47R1 based hindcast using ERA5 forcing for the 11 research cruises (Yang et al., 2022).The RMSE between our simulation and the ECMWF data set is 0.18 m.The agreement between the two distinct models indicates a reasonable wave model set-up.
The observational data set also provides sea surface temperature measurements, which is used to estimate gas diffusivity and water viscosity, and the associated Schmidt number Sc and the gas solubility α.The CO 2 solubility and Schmidt number as a function of temperature and salinity is calculated following Wanninkhof (2014) and Shin et al. (2022).The viscosity is computed following Nayar et al. (2016).For the gas transfer velocity simulation for the 11 cruises, the sea surface salinity is assumed as a constant 35 psu.

Optimizing Uncertain Coefficients in the Bulk and Spectral Approach From Observations
We optimize remaining uncertain parameters within the spectral and bulk models for parameterizing the air-sea gas transfer velocity (described in Section 2) by performing a statistical analysis against the observational data set.We consider three statistical metrics: the root mean square difference, the linear fit coefficient and the R 2 .
For the spectral approach, we explore multiple combinations of A NB over a range of  B from = 0.1 to 0.15.We compare   660  computed using multiple values of  B and A NB to field measurements   660  , testing multiple suitable combination of  B and A NB (24 combinations, not shown) that give the most comparable spectral CO 2 transfer velocity to the observational CO 2 transfer velocity.We consider minimizing the root mean square difference, optimizing the slope of the linear fit, and maximizing R 2 together between observations and the parameterization to chose the best fit for parameters A NB and  B for the spectral method.We illustrate this discussion by comparing four candidate combinations in the upper row of Table 1.We use A NB = 1.38 × 10 −4 for non-breaking component in this study, and we apply  B = 0.12 in the spectral sea state model (Equation 8).
Next the optimization for A b in the bulk sea state CO 2 transfer velocity model proceeds by testing A b ranging from 1.0 × 10 −5 to 2.0 × 10 −5 s 2 /m 2 .Table 1 shows the four examples of the statistical analysis of the bulk sea state CO 2 transfer velocity model-data comparison (lower row).Setting A b = 1.2 × 10 −5 s 2 /m 2 in the bulk sea state model (Equation 9) for the bubble mediated gas transfer velocity yields the reduced error statistics considering the error statistics and the linear fit.The uncertainty from above coefficient is discussed further in Section 5.2.
Figure 2 shows the spectral sea state dependent gas transfer velocity , where the slope coefficient can reach to 0.75.This coefficient is closer to unity than the coefficient 0.69 for the wind-only parameterization in panel c in Figure 2, and is similar with the coefficient 0.74 for bulk sea state dependent parameterization in panel b in Figure 2. We note that the agreement between the gas transfer velocity parameterization from Wanninkhof (2014) and the data can also be improved by calibrating the prefactor to the data set as done by Yang et al. (2022).Here we are interested in comparing the formulation commonly used in ocean and climate models (Wanninkhof, 2014) and sea-state parameterizations that can explain variability at the same wind speed.We note that choosing a higher  B would lead to a linear fit closer to 1, but would have a higher root mean square difference.In general, the gas transfer velocity from the three parameterizations agrees well with available observations from 11 cruises with the error bars almost always intersecting the 1 : 1 line (black line in all 3 panels in Figure 2).

Uncertainties on the Mean Related to Model Parameters
Now that we have established the optimized coefficients to be used, we evaluate the uncertainties around the mean for these parameters.This discussion is sensitive to the way the data averaging is done.In this study, we do averaging over all data sets including 11 cruises in each wind bin.Yang et al. (2022) averaged for each cruise in each wind bin and then average the averaged value from each cruise.
The panel a in Figure 3 shows the mean value of   660  with  B = 0.12 in each 2.5 m/s wind bin as a function of wind speed.The red vertical errorbar is the uncertainty of   660  when  B ranges from 0.1 to 0.2.The blue vertical errorbar is the uncertainty of   660  when injection depth z 0 ranges between 0.5H s and H s .The black symbol in Figure 3 is the mean observational gas transfer velocity   660  in each 2.5 m/s wind bin, the black vertical errorbar is the standard deviation of   660  in each wind bin.This result from the panel a in Figure 3 shows that the spectral sea state dependent parameterization yields higher gas transfer velocity than the wind-only dependent    with  B = 0.1 (lowest limit of red errorbar in panel a in Figure 3) is still above the curve from (Wanninkhof, 2014).It is also more consistent with the observational mean   660  .For the bulk parameterization, the uncertainty is mostly within the evaluation of the coefficients A b and A NB .The panel b in Figure 3 shows the mean gas transfer velocity for CO 2 from the bulk sea state dependent formula   660  with A NB = 1.38 × 10 −4 and A b = 1.2 × 10 −5 s 2 /m 2 in each 2.5 m/s wind bin as a function of wind speed.The magenta vertical errorbar indicates the uncertainty of   660  for A b = 1.2 ± 0.2 × 10 −5 s 2 /m 2 , the bottom limit shows the mean value of   660  with minimum A b = 0.8 × 10 −5 s 2 /m 2 .The upper limit of green vertical errorbar is the mean value of   660  with maximum A NB = 1.55 × 10 −4 in each wind bin.Even the minimum mean   660  is also larger than wind-only parameterization   660 14 .Both spectral and bulk sea state dependent gas transfer parameterizations yield larger values than the traditional wind-only parameterization at high wind speed even with considering uncertainties from A NB ,  B , z 0 and A b .Overall, the two sea state dependent approaches better represent the observational data sets obtained by eddy covariance measurements summarized by Yang et al. (2022).

Comparison for Individual Cruise Data
We now investigate how well all three parameterizations perform in different wind regimes.Figure 4 shows the spectral sea state gas transfer velocity (red dots), bulk sea state gas transfer velocity (green dots), and observational gas transfer velocity gray dots) as a function of wind speed.The 1 m/s wind bin averaged spectral sea state gas transfer velocity (red square), bulk sea state gas transfer velocity (green square) and observational gas transfer velocity (gray square) are also shown.Vertical error bar indicates the 95% confidence level in each wind bin.The black curve indicates the wind-only gas transfer velocity from Wanninkhof (2014).The data including all 9 data sets from 11 cruises are shown in panel a in Figure 4.The individual data sets are shown separately in the different panels (panel b to k). Figure 4 shows that the modeled data from all three parameterizations reproduces the trend and range of field measurements well.The traditional wind-only parameterization  2014) underestimates the gas transfer velocity for CO 2 when U 10 > 10 m/s.While the sea state dependent formulation reproduces some variability around the mean, the observational data set from the HiWinGS experiment (panel (h) in Figure 4) clearly has larger variability than any of the ones obtained from modeling.
Figure 4 suggests the main difference between the sea state dependent gas transfer velocity and observational gas transfer velocity is from two cruises: Knorr11 and HiWinGS.The root mean square error of Knorr11 and HiWinGS are 16.35 cm/hr and 23.61 cm/hr (not shown).Both cruises happened in the winter storm season in North Atlantic (Yang et al., 2022).For the high wind condition with U 10 > 20 m/s, when the observational   660  > 100 cm/hr, uncertainty becomes large due to limited observations.

Comparison Between Bulk and Spectral Approach
The difference between the spectral and bulk sea state dependent CO 2 transfer velocity for all 9 data sets from 11 cruises are shown in Figure 5.The difference between the spectral sea state dependent gas transfer velocity   660  and bulk sea state gas transfer velocity   660  for cruises ANDREXII, HiWinGS and JR10087 can reach up to 50 cm/hr (see panel a in Figure 5).The percentage of the above difference of spectral sea state dependent gas transfer velocity can reach to 50%, as shown in panel b in Figure 5.
Figure 6 shows the percentage of difference with the color shading indicating the significant wave height (panels a, c and e) and air entrainment volume (panels b, d and f) for the three cruises.The significant differences happens when the wind speed is lower than 10 m/s.The differences occur when the significant wave height is high, but the wave breaking statistics are low.The high significant wave height with low wind condition indicates the presence of non-local swell, which is not expected to contribute to enhanced wave breaking in nature.The swell induced high significant wave height can therefore introduce overestimation of the gas transfer velocity from the bulk sea state dependent gas transfer velocity model.For the HiWinGS cruise, the difference of   660  from the two sea state dependent models happens when 7 ≤ U 10 ≤ 15 m/s.This is not likely due to the same swell effect as in ANDREXII.The major reason for this difference is that the wind condition during the HiWinGS cruise increases suddenly, this rapid increase in wind can lead to young and steep waves that break frequently, which is captured by the spectral approach but not the bulk approach (Brumer, Zappa, Brooks, et al., 2017).
The key differences between the spectral and bulk approach under high wind conditions suggests the high wind condition needs more observational data to better evaluate and constrain the sea state dependent models.In general our results suggest that the bulk sea state gas transfer velocity works well to represent the gas transfer velocity predicted by the spectral model from the wave breaking statistics.However, these special conditions (strong swell, rapid increasing wind phenomenon, and high wind speeds) will require further evaluation in future studies.We note that the significant wave height is mostly determined by the low frequency part of the wave spectrum (being related to the 0 moment of the wave spectrum), and is not necessarily the most relevant parameter to describe the breaking wave statistics.It is however one of the simplest wave metric and one that should become widely accessible from ongoing and future satellite missions.An integrated parameter more directly describing the breaking waves could be the mean square slope (being related to the second moment of the wave spectrum) see the discussion in Wu et al. (2023).Such a parameter could be computed from the WW3 and another gas transfer parameterization could be proposed, and we expect would have behavior more similar to the spectra formulation.

Sea State Dependent Gas Transfer Velocity
We now apply the two sea state dependent parameterizations (as optimized in Section 4) and the wind-only parameterization for gas transfer velocity for CO 2 along with historical global wave simulations from 1980 to 2017.The output fields are stored from WW3 in this study at 3 hr intervals.The time varying Sea Surface Temperature (SST) and Sea Surface Salinity (SSS) are also needed to calculate the kinematic viscosity and Schmidt Number in the sea state dependent CO 2 gas transfer velocity parameterizations (see Section 2).These fields are provided at daily resolution from a corresponding Geophysical Fluid Dynamics Laboratory OM4.25 ocean/sea ice model simulation (Adcroft et al., 2019).
Figure 7, panels (a) and (d), show the mean spectral sea state dependent gas transfer velocity   660  in Boreal winter (DJF) and summer (JJA) over the full time period.The global pattern of gas transfer velocity shows high gas transfer velocity in three regions of high wind-wave activity: North Atlantic, North Pacific and Southern Oceans, which is consistent with Reichl and Deike (2020).The mean difference between the spectral and bulk sea state parameterizations,   660  −  660  , and the spectral sea state parameterization and wind-only parameterization,   660  −  660 14 , in Boreal winter (summer) are shown in Figure 7, panels b and c (e and f).The differences suggest that the bulk sea state dependent parameterization predicts a lower gas transfer velocity than the spectral gas transfer velocity in the region with high wind-wave activity by ∼ 3 cm/hr, but a larger value in the tropical and subtropical ocean as shown in Figure 7 (panels b and e).The spectral sea-state dependent approach is ∼8 cm/ hr higher than the wind-only parameterization in the high wind and wave regions, as shown in Figure 7, panels (c) and (f).This result suggests that the sea state dependence of the bubbles induced by wave breaking are important for gas fluxes and their variability, as suggested by Reichl and Deike (2020).Furthermore, it suggests that accounting for the spectral properties in addition to the wave height better accounts for the sea-state variability in these high wind and wave regions and can appreciably contribute to variability in gas transfer velocity.There may therefore be some benefit to using the full spectral approach when specific studies require more accurate estimates of the gas transfer velocity and its variability.) in winter (c) and summer (f).

Sea State Induced Variability in Gas Transfer Velocity
Following Reichl and Deike (2020), we introduced the sea state induced gas transfer velocity variability: sea state dependent gas transfer velocity against U 10 .The difference between the mean spectral and bulk sea state dependent gas transfer velocity is less than 2 cm/hr when U 10 < 10 m/s, but the difference becomes larger with wind speed increasing from 15 m/s to 20 m/s.The mean bulk sea state dependent gas transfer velocity is larger than the mean value from the spectral approach with wind speeds higher than 22 m/s (Figure 8, right panel).The wind-only dependent gas transfer velocity  (  660

𝑤𝑤𝑤𝑤𝑤𝑤𝑤𝑤𝑤14
) is smaller than the mean value of both sea state dependent gas transfer velocity approaches when U 10 > 10 m/s.In our result, the wind only parameterization of gas transfer velocity tends to fall below the observed bin mean   660  in this wind speed range (see Figure 3), so our analysis indicates that the wind-only formula of gas transfer velocity tends to underestimate the gas transfer velocity for U 10 > 10 m/s compared to the sea state dependent methods.Note, Yang et al. (2022) use a different averaging method for determining the bin mean observed   660  (by averaging first per field campaign and then perform an ensemble average of the campaigns) that instead shows better agreement with   660 14 in this wind speed range.
The climatological mean sea state induced variability from the spectral formula in Boreal winter and summer are shown Figure 9 (panels a and c).The global pattern is similar to that demonstrated from the sea state effect by Reichl and Deike (2020).The sea state introduces stronger gas transfer velocity in the eastern North Pacific, eastern North Atlantic region and Arctic open ocean in Boreal winter.The sea state induced variability is negative in most coastal regions and in the tropical oceans, which can be regions of low waves.In Boreal summer, the sea state induced variability is small in most regions, but the gas transfer velocity is a little stronger on average in the Southern Ocean, south of 40°S (Figure 9, panels c and d).The climatological sea state induced variability from the spectral formula is not as large compared to the bulk formula in the high-latitude northern hemisphere (Latitudes >40°N).This result also shows that the spectral sea state induced variability of gas transfer velocity is weaker than the sea state induced variability in the bulk model in the western margins of the basins in winter (Figure 9, panels a and b).The difference between the spectral and bulk sea state induced gas transfer velocity variability is generally smaller in summer comparing with winter season (Figure 9, panels c and d).

Discussion and Conclusion
In this study, we revisit the sea state dependent gas transfer velocity parameterizations developed by Deike and Melville (2018).We incorporate the wave breaking statistics determined from the wave spectrum (Romero, 2019;Shin et al., 2022) with WW3 into the bubble mediated gas transfer model from Deike and Melville (2018) based on the chemical model from Keeling (1993).We compare two approaches for the sea state dependent parameterization, one based on the wave spectrum (spectral approach) and one on the significant wave height (bulk approach), for the gas transfer for CO 2 from Deike and Melville (2018).The two approaches are evaluated against the existing wind only parameterization using a large range of eddy covariance field measurements in various ocean conditions (see Yang et al., 2022).The two bubble mediated gas transfer velocity approaches are optimized using these field measurements, and we demonstrate that both sea state dependent formulation agree well with the observational data sets, performing better than wind only formulation.The two sea state dependent approaches capture some of the variability observed in the field data at given wind speed.Our results confirm that the bubbles contribute significantly to the total gas transfer velocity with increasing wind speed for CO 2 .The optimized coefficients from the field data are provided in Table 1 for the bulk and spectral k w models.
We archive the 3 hourly climatology of sea state dependent gas transfer velocity spanning 1980 to 2017 using both the spectral sea state dependent parameterization and the bulk sea state dependent parameterization.The climatology shows that sea state dependent bubbles play a significant role for the North Atlantic, North Pacific and Southern Ocean for net CO 2 gas transfer velocity.The sea state induced variability can approach 3 cm/hr, around 6% of the total sea state dependent gas transfer velocity (often ∼50 cm/hr) in regions that are north of 40°N and south of 40°S.This variability could potentially influence the global CO 2 distribution, particularly when considering the global ocean circulation and interior water mass formation regions, as discussed in Reichl and Deike (2020).
We analyze how the sea state dependent gas transfer velocity models differ from a traditional wind-only parameterization, as well as the difference between the spectral and bulk sea state dependent gas transfer velocity parameterizations.The spectral formulation accounts for the high frequency breaking dynamics of the wave field, controlled by local and non-local dynamics, and presents the best comparison to the field data.The bulk formulation, based on the significant wave height, performs generally quite well in our analysis.We do identify a high gas transfer velocity bias in the bulk approach at low wind speeds and in the presence of large swell that is not found in the spectral approach, likely due to the bulk formula not distinguishing between non-breaking swell and the wind sea.The gas transfer velocity at high wind speed, typically above 20 m/s, remains uncertain at least partially due to the limited observational data sets.This is a typical phenomena in Southern Ocean, especially in winter.We note that tropical cyclone wind speed conditions (U 10 > 35 m/s) gas transfer velocities have not been measured in the field, while some laboratory experiments have been reported (Krall et al., 2019;Stanley et al., 2022).
These results suggest that using the fully spectral model should be considered in future coupled ocean-wave models.Separately, the bulk formulation presents the advantage of being easier to implement and only requires knowledge of the significant wave height that can be provided by climatology (e.g., as obtained from the present study), and is could be available through global satellite coverage.We note that the sea-state dependent formulations we present here yield similar annual mean gas transfer velocity to the wind only formulations used widely in global models.The main difference resides in the ability to capture added variability; as well as providing an understanding of the dependencies in wind and wave to isolate contributions from the physical-chemical parameters of the gases, in order to provide a universal framework for all gases.
In summary, the spectral and bulk sea state dependent gas transfer models both increase the overall agreement between the modeled and observed gas transfer velocity compared to the wind only method.Both approaches exhibit some of the variability observed in field measurements of the gas transfer velocity, increasing confidence in their ability to parameterize physically meaningful processes.While the sea-state dependent gas transfer velocity calculation requires more steps than the wind only approach, the routines remain very small compared to other aspects of coupled Earth System Models.Therefore, the cost of coupling the sea state dependent gas transfer velocity is almost entirely if the model is coupled to a wave model (which may also be desired for other aspects of the simulation).However, there are alternative approaches to using a full numerical wave model, such as parameterized wave spectra or offline coupling with preexisting wave data sets.
To conclude, we highlight remaining limitations of the present approach and discuss future research avenues.The wind friction velocity is derived from the wind speed at 10 m through the drag coefficient, here using the COARE3.5 (Fairall et al., 2003) formulation without considering sea state effects (such as misalignment between wind and waves).Incorporating sea state effects into the drag coefficient formulation is an active area of research and could influence the non-bubble (and bubble) variability especially at low wind speed (Janssen & Bidlot, 2023;Sauvage et al., 2023).We found three phenomenons: (a) the strong swell, (b) rapidly increasing wind, and (c) high wind conditions (U 10 > ∼20 m/s), that can introduce a difference between the spectral and bulk sea state gas transfer velocity, as discussed in detail in Section 5.2.Observational data at high wind speed (U 10 > 20 m/s) remain limited due to the difficulties in taking measurements in such conditions, introducing relatively large uncertainties when comparing model and observations.Efforts to collect more observational gas transfer velocity data sets under high wind condition would be helpful to constrain the uncertainty found in this and previous studies.Effects of surfactants on gas exchange have been discussed (Bell et al., 2017), but still need to be incorporated within this modeling framework.Finally, the dependencies on gas solubility needs to be further investigated, in particular to incorporate effects of asymmetric gas exchange (Leighton et al., 2018;Liang et al., 2020;Stanley et al., 2009Stanley et al., , 2022)).

Appendix A: Table of Symbols
All symbols in this study are described here.See Table A1  for all 11 cruises are available on Zenodo (Zhou, 2023b).The WW3 climatology data set for this study is available on Princeton DataSpace (Zhou, 2023c).Code to compute bubble mediated gas transfer velocity can be found at GitHub (Zhou, 2023a).

Figure 1 .
Figure 1.(a) Tracks of all research cruises in Yang et al. (2022) with climatological mean significant wave height from WW3 simulations from 1980 to 2017 as background.(b) JRA55-do 10 m wind speed against in-situ observed 10 m wind speed in all cruises in panel (a).(c) The surface wind frictional velocity calculated from JRA55-do U10 and COARE3.5 against the surface frictional velocity in Yang et al. (2022).(d) The significant wave height (H s ) from WAVEWATCH III simulations against H s from ECMWF hindcast which is used in Yang et al. (2022).The root mean square error (RMSE) between different products are labeled.
is the first explicit implementation of the wave breaking statistics within a spectral wave model, computing the wave breaking distribution Λ(c) from the wave spectrum computed within WW3.The WW3 model from Romero (2019) has been tuned against ST4 solutions of H s and T p as stated in Romero (2019), so that both formulations lead to similar wave fields.The distribution of length of breaking crest can be expressed in terms of the wave phase speed c or wavenumber k, with Λ(k)d(k) = Λ(c)dc, with   = √ ∕ given by the gravity wave linear dispersion relation.Romero (2019) proposed the following functional form  Λ() =   exp(−∕())() () A NB = 1.48 × 10 −4 A NB = 1.33 × 10 −4 A NB = 1.38 × 10 −4 A NB = 1.38 × 10 −4

Figure 3 .
Figure 3. Spectral gas transfer velocity for CO 2 with uncertainty from the parameter  B in Equation 5 (red errorbar) and z 0 (blue errorbar) (panel (a)).Bulk gas transfer velocity for CO 2 with uncertainty from coefficients A b (purple errorbar) and A NB (green errorbar) as a function of wind speed U 10 (panel (b)).Black symbol is observed gas transfer velocity from the 11 cruises in Figure 1  (  660  ).The black curve indicates the wind-only dependent gas transfer velocity   660 14 .

Figure 2 .
Figure 2. Comparison of gas transfer velocity for CO 2 between (a)   660  , (b)   660  , (c)   660 14 and observational   660  from all cruises in Figure 1(a).The bin-averaged data with k w,obs every 20 cm/hr are shown as gray symbols.The errorbar of bin averaged data is 95% confidence level.The root mean square error is labeled in each panel.The linear fit between modeled gas transfer velocity and observational gas transfer velocity are shown as blue line in each panel, and labeled as blue text.The correction coefficient squared (R 2 ) is also labeled as red text in each panel.The color of round circles indicate the wind speed.
(2014), even accounting for uncertainties in the coefficients  B and z 0 , since the lowest value of   660

Figure 4 .
Figure 4. Modeled and measured   660  from all 11 research cruises (9 data sets) (panel (a)), individual data sets are shown in different panels.The black symbols are the observational CO 2 transfer velocity   660  (panel (b) to panel (k)).Red symbols are the spectral sea state dependent gas transfer velocity   660  .Yellow green symbols are the bulk sea state dependent gas transfer velocity   660  .The bin-averaged label is averaged over each 2.5 m/s.Errorbar indicates the 95% confidence level for each wind bin.

Figure 5 .
Figure 5.The difference between spectral and bulk sea state dependent CO 2 transfer velocity (   660  −  660  , (a)) and percentage of difference between spectral and bulk sea state dependent CO 2 transfer velocity (   660  − 660   660  , (b)) against wind speed U 10 for all 9 data sets (11 cruises).Different color indicates different data sets as shown in labels.

Figure 6 .
Figure 6.The percentage of difference for spectral and bulk sea state dependent CO 2 transfer velocity  (  660  − 660

Figure 8 .
Figure 8. (left) Climatological mean gas transfer velocity in 1 m/s wind speed bins from spectral sea state gas transfer velocity model  (  660  ) and bulk sea state gas transfer velocity model  (  660  ).The black curve indicates the wind-only

Figure 9 .
Figure 9. Climatological mean sea state induced variability from the spectral gas transfer velocity  ( (  660  ) ′ ) in Boreal winter (December, January and February (DJF), panel (a) and summer (June, July, August, (JJA), panel (c); The climatological mean sea state induced variability from bulk gas transfer velocity  ( (  660  ) ′ ) in Boreal winter (panel b) and summer (panel d).
The same error metrics for bulk CO 2 transfer velocity model and data comparison with variable A b in units s 2 /m 2 (bottom row).P(1) and P(2) indicate the slope and intercept of the linear fit, respectively.The Correlation Coefficient Squared (R 2 ), Root Mean Square Error (RMSE) and the Coefficient for Linear Fit Between the Spectral Sea State Dependent Gas Transfer Velocity Model With Variable  B and A NB and Observational CO 2 Transfer Velocity (Upper Row) .