Dimensionless Parameterizations of Air-Sea CO2 Gas Transfer Velocity on Surface Waves

1 Introduction

Cumulative anthropogenic Carbon Dioxide (CO₂) emissions have been considered responsible for global warming and climate extreme events. The ocean is dynamically exchanging CO₂ with atmosphere and has taken up about 26% of total human activity-induced CO₂ emissions over the last decade (Friedlingstein et al. 2022). The air-sea CO₂ flux (F) in a bulk formula is generally expressed as a function of gas transfer velocity (K_CO2), aqueous solubility (s) and partial pressure difference across the air-water interface:

where pCO_2w and pCO_2a are CO₂ partial pressure in surface water and air, respectively. To reduce uncertainties in the estimation of air-sea CO₂ fluxes for global carbon budget, high-quality data sets of CO₂ partial pressure and accurate parameterizations for quantifying gas transfer velocity are required. Major advances have been made in collecting (Bakker et al. 2016; Takahashi et al. 2019) and reconstructing (e.g., Landschützer et al. 2020) observed ocean surface CO₂ concentration. The data sets with improved accuracy and temporal and spatial coverage provided notable support for the assessment of decadal trends and variabilities of historical air-sea CO₂ fluxes (Landschützer et al. 2016; Friedlingstein et al. 2022).

CO₂ transfer velocity K_CO2 depends on environmental forcings including wind, surface waves, bubbles, surfactants, etc. The parameterizations of K_CO2 based on wind or wave parameters have been derived in previous theoretical, laboratory and field studies. Gas transfer velocity is inversely correlated with a total resistance which is the sum of gas and liquid resistances. Because CO₂ is sparingly soluble in water, transfer resistance lies mainly in the water side (Bolin 1960; Liss 1973). Molecular diffusion (for tranquil water body) and turbulent mixing in water boundary layer are dominant mechanisms for transfer efficiency. Considering that the kinetic energy input to water is originated from wind stress, typical parameterizations of K_CO2 are developed in terms of wind speed (U^m) with power m ranging from 1 to 3. For example, by employing laboratory and field (lake) observations, Liss and Merlivat (1986) proposed a piecewise linear function regarding smooth water surface, wavy surface and wave breaking surface under varied wind regimes. Thereafter, quadratic dependences were developed to estimate K_CO2 under intermediate wind speed (Wanninkhof 2014; Nightingale et al. 2000; Ho et al. 2006). For high wind speed with wave breaking and bubble entrainment, cubic relationship was adopted (Wanninkhof and McGillis 1999). With the development of eddy-covariance method for ocean observations, the gas transfer velocity with wind speed has been actively investigated (McGillis et al. 2001; Miller et al. 2010; Landwehr et al. 2018; Dong et al. 2021; Yang et al. 2022). It is also worth noting that in laboratory, K_CO2 was also measured under extreme wind conditions with 10m wind speed (U₁₀) up to 70 m/s by Iwano et al. (2013) and 85 m/s by Krall et al. (2019). They established a new regime of U₁₀ beyond 33 m/s in which K_CO2 was largely enhanced. The magnitude of power of wind speed and constants used in previous wind-based formulae were tuned for a mean dependence of K_CO2 on wind speed. These empirical relationships can yield reasonable K_CO2 estimations within their scope of validity but have uncertainties induced by complex physical processes in water (e.g. wave breaking, bubbles). Gaps exist across the formulae especially under intermediate to high winds (Wanninkhof 2014; Zhao et al. 2003), which reveals the inadequacy of wind-only parameterizations.

For CO₂ gas transfer, the turbulent intensity in water is of vital importance. With the wind stress acting on sea surface, a part of wind energy is transferred directly to the drift currents through viscous stress while the other part supplies the growth of waves. By applying a rigid boundary and assuming that turbulent dissipation is balanced by the shear production, a “wall layer” with constant stress and logarithmic velocity profile exists at water surface. Turbulent kinetic energy (TKE) in this layer can be simply scaled with the cube of the water-side friction velocity and the TKE dissipation rate is inversely proportional to the depth (Craig and Banner 1994). However, near-surface turbulence can be significantly enhanced by waves and wave breaking. It has been proposed (Qiao et al. 2004; Babanin 2006) and verified (Babanin and Haus 2009; Dai et al. 2010; Babanin and Chalikov 2012) that waves, even without breaking, can produce turbulence through wave orbital motions when water-side Reynolds number exceeds a critical value. The turbulence induced by non-breaking waves can be effectively strong and affect other ocean processes such as the upper-ocean mixing (Toffoli et al. 2012). In spectral waves, energy is transferred among the wave scales (i.e. wavenumbers) through nonlinear wave-wave interactions and subsequently dissipates via wave breaking, the processes of which are in balance with the wind input in equilibrium range of wave spectrum (Phillips 1985). Recognizing the important role of waves and wave breaking for TKE, the vertical profile of TKE dissipation rate and its wave-dependent scaling have been studied (Terray et al. 1996; Babanin et al. 2005; Gemmrich 2010; Sutherland and Melville 2015; Thomson et al. 2016; Lee et al. 2017).

In addition to the turbulence, wave breaking also generates bubble plumes which offer extra interface with increased partial pressure in bubbles and are important especially for sparingly soluble gases (Bell et al. 2017). The behavior of injected bubbles is subject to the turbulence in the vicinity. Melville et al. (1992) found that up to 50% of the dissipated energy from wave breaking could be expended on air entrainment against buoyancy forces. The large bubbles are broken up into smaller ones by turbulent fluctuations. This process stops at the Hinze scale (Hinze 1955; Garrett et al. 2000) at which surface tension of bubbles prevents further breakup. Through dimensional analysis, Garrett et al. (2000) proposed that the bubble size spectrum was related with local TKE dissipation rate, air supply rate and bubble radius. Deike et al. (2016) further parameterized the total air entrainment with breaking wave parameters for practical application.

The correlation between gas transfer velocity and wave-related dynamic processes (i.e. turbulence, bubbles) has been investigated by previous studies. Based on surface renewal theory, Lamont and Scott (1970) found that gas transfer velocity was proportional to TKE dissipation rate with a power of 0.25. This relationship has been supported by both theoretical and field studies (e.g., Lorke and Peeters 2006; Zappa et al. 2007). Jähne et al. (1987) suggested that mean square slope of the waves was a proper parameter to characterize gas transfer velocity. In laboratory, Zappa et al. (2001; 2004) observed that the fractional area coverage of microwave breakers under low to intermediate wind speed was another suitable parameter. At ocean surface, whitecaps are commonly observed in the process of wave breaking. Zhao et al. (2003) and Woolf (2005) employed whitecap coverage as a proxy to scale gas transfer velocity. Furthermore, whitecap coverage in these studies was scaled with a wind-sea Reynolds number rather than wind speed alone. The Reynolds number, denoted as R_H or R_B, is non-dimensional and can be written as

where u_* is wind friction velocity, H_s is significant wave height, ν_a is air kinematic viscosity, and ω_p is peak angular frequency of wave spectrum. R_H and R_B are interpreted as a measure of turbulence generated by wind waves and wave breaking. Their formulae imply that both wind force and sea state are equally important. Brumer et al. (2017) adapted the Reynolds number by changing ν_a to ν_w (water kinematic viscosity) and reconciled gas transfer data sets collected from field observations. Toba et al. (2006) discussed the similarity between air-sea momentum flux and gas flux, and proposed a non-dimensional treatment (divided by wind speed) for gas transfer velocity which was further correlated with the wind-sea Reynolds number. Li et al. (2021) observed the CO₂ gas transfer in laboratory under scenarios of mechanically-generated modulational wave trains without wind forcing, wind-generated waves and the combination of modulational waves and superimposed wind. They formulated non-dimensional parameterizations in which wave parameters were dominant.

Bubble-mediated gas transfer has been modeled and observed in previous studies (Woolf and Thorpe 1991; Keeling 1993; Woolf et al. 2007; Bell et al. 2017; Zavarsky et al. 2018) and is commonly scaled with the wind speed (Stanley et al. 2009; Liang et al. 2013) or whitecap coverage (Woolf 2005). Deike and Melville (2018) proposed a spectral framework for bubble-mediated gas transfer, in which wind friction and significant wave height were used. Based on this formula, Reichl and Deike (2020) estimated that global CO₂ gas fluxes via bubbles accounted for 40% of the total air-sea fluxes.

In the present study, we propose wave-dominated formulae and conduct validations using combined laboratory and field data. A simple bubble-mediated gas transfer scaling is obtained through field measurements and integrated into our formula. In addition, the efficiency of previous parameterizations is evaluated.

2 Data and Methods

Synchronous measurements of gas transfer and environmental forcings from four field projects are collected in this study. To avoid the inconsistency caused by different observational methods, air-sea CO₂ fluxes across the projects were obtained through direct eddy-covariance measurements and the surface elevation for wave field was measured by Riegl laser altimeter. The air-sea CO₂ partial pressure was recorded by similar underway equilibrator system (e.g., (Pierrot et al. 2009). The local wind speed was measured by sonic anemometer. In addition, the ship motion, air and water temperature, air pressure and relative humidity (RH) were synchronously recorded during those campaigns. Except for the adopted hourly results from previous studies, variables are calculated in 10-minute segments with 50% overlap, which results in 11 pieces in one hour. Then the segments are averaged hourly to reduce the instability and deviation. The hourly results are further averaged into equidensity bins containing 15 data points.

The Southern Ocean Gas Experiment (SOGASEX) was conducted in Southern Ocean near South America from February to April in 2008 (Figure. 1 (a)). The region was selected for large air-sea CO₂ partial pressure difference to ensure a high signal-to-noise ratio of the eddy-covariance flux. CO₂ flux was measured by both closed and open path nondispersive infrared gas analyzers (NDIR). The 10m wind speed was up to about 20 m/s with mean value of 9.7 m/s (Edson et al. 2011). The mean RH was 87% and more than 25% of records were higher than 96% which is the upper limit of reliable RH measurements used. The final CO₂ transfer velocity computed following equation (1) is provided by Brumer et al. (2017).

Figure 1

The High Wind Gas Exchange Study (HIWINGS) took place at Labrador Sea from October to November in 2013 (Figure 1 (b)). The region at that period was subject to high wind speeds, frequent storms and a well-known sink for atmospheric CO₂ with large partial pressure gradient. The wind speed was up to about 25 m/s. The wave field was measured by a Riegl laser altimeter and a Datawell Waverider buoy (model DWR-G4). Bubbles were also measured by a bubble camera (for big bubbles) and a bubble resonator located on spar buoy (for small-sized bubbles). Bubble measurement was operated for 40-minute duration at 3-hour intervals over a maximum of 48 hours in a single deployment. The CO₂ transfer velocity and bubble injection rate are provided by Blomquist et al. (2017).

CO₂ gas transfer was also measured in the project CAPRICORN during the Voyage IN2016_V02 from March to April 2016 in Southern Ocean near Australia (Figure 1 (c)). This project aims to explore Southern Ocean cloud systems, aerosol properties, surface energy budget, upper ocean biological aerosol production, and atmospheric composition. CO₂ flux was measured along the campaign through open path NDIR (LI-COR, LI-7500). CO₂ partial pressure was continuously measured by using a General Oceanic/Neill system which was equipped with a NDIR (LI-COR, LI-7000), shower head equilibrator and Nafion Dryer (Moreau et al. 2017). The wind speed was up to about 25 m/s. The wave field was also monitored through a Riegl laser altimeter. The observed CO₂ flux was processed following the method used by Blomquist et al. (2017).

DYNAMO program (Figure 1 (d)) was conducted during year 2011–2012 in tropical Indian Ocean where Madden Julian Oscillation (MJO) as a result of ocean atmosphere interaction can affect the earth’s climate. The ocean in this region was a weak source for atmospheric CO₂ with highest wind speed up to about 15 m/s. CO₂ flux was measured by both closed (LI-7200) and open (LI-7500) path NDIR. In this work, we employed the results from one LI-7500. The assessment of CO₂ flux follows Blomquist et al. (2014). CO₂ partial pressure was monitored through a system comprising an equilibrator, a CO₂ analyzer (LI-COR 840) and a Nafion air dryer.

Laboratory measurements of CO₂ transfer under different wind and wave conditions from Li et al. (2021) are also used to evaluate parameterizations. The experiment was designed to explore the role of waves in air-water CO₂ gas exchange. Three types of waves were forced in a flume including modulational wave trains generated by mechanical wave maker, pure wind waves and the superposition of mechanical waves and wind waves. Surface elevation was recorded by wave gauges. Wave breaking was also captured by cameras. CO₂ was measured by Apollo system (model AS-P2) which incorporated a shower head equilibrator, a multiposition valve, and a CO₂ analyzer (CRDS, model G2301 by Picarro). The wave sizes in the laboratory are evidently smaller than that observed in ocean. An appropriate non-dimensional scheme should be able to reconcile all data sets.

3 Revised Parameterizations for CO₂ Transfer Velocity

The evaluations of equation (4) and (9) are presented in Figure 2 (a) and (b). The error bars denote standard deviation from the average of binned data. In panel (a), equation (4) is able to collapse all data well although scatters still exist. The fitting coefficients α and β are modified as 9.57 · 10^–11 and 0.876, respectively. The coefficient of determination (r²) and root-mean-square error (rmse) are also computed. Laboratory data (blue circles) are particularly in a good agreement with HIWINGS observations (red circles), for which r² is high as 0.96 (not presented). The scatter is mainly introduced by SOGASEX (green squares) and DYNAMO data (pink triangles) of which the deviations denoted by error bars are also relatively significant. Equation (9) in panel (b) apparently fails to reconcile all data sets. Gas transfer velocity from field observations are higher than the laboratory measurements under similar forcings computed by equation (9). The initial motivation for formulating equation (9) in laboratory is to highlight the importance of wave breaking. Gas transfer rate is regulated by the frequency of breaking events (b_T). However, R_HM is interpreted as an indicative turbulence parameter without accounting for the bubble’s contribution (Li et al. 2021). Bubble-mediated gas transfer could be more evident in ocean than in laboratory because the magnitude of wave breaking in ocean can be larger. Thus, we propose to introduce bubble’s effect into equation (9). Because the injected bubble volume and size distribution are closely correlated with environmental turbulence, it was suggested that bubble volume could be scaled with wave breaking dissipated energy (Fairall et al. 2011; Long et al. 2011). Based on the bubble injection rate V_b (unit m/s) from HIWINGS, we can simply parameterize it with R_HM. From dimensional considerations, V_b is firstly scaled with wave orbital velocity U_wm (m/s) and the dimensionless injection rate is denoted as Ṽ_b. Ṽ_b is found proportional to $R_{HM}^3$ in Figure 2 (c). Although the correlation is evident (78%, not shown), r² is 0.47 with limited data points. Then, the bubble’s effect in terms of $R_{HM}^3$ is incorporated into equation (9) and the new formula is written as

where R_HM has a total power of 4. Figure 2 (d) shows a substantial improvement for equation (10) in reconciling all data. Computed r² is 0.7, higher than that of equation (4) in panel (a). The adjusted fitting coefficient α and β are computed as 2.82 · 10^–11 and 0.260, respectively.

Figure 2

With the collected data, we can also evaluate other dimensional parameterizations for CO₂ transfer velocity. In Figure 3 (a), evaluations of several popular wind-only formulae (Liss and Merlivat 1986; Wanninkhof 2014; Nightingale et al. 2000; McGillis et al. 2004; Edson et al. 2011) are presented. Although these formulae generally comply with field observations, they can not reduce the disparities between field and laboratory results. Indeed, even for measurements of two groups of experiments (blue and green circles) in laboratory, different wave states result in clear gaps under similar wind speeds. The results of DYNAMO (pink triangles) also show significant deviations from the predictions of wind-based formulae. The combined evidences could demonstrate the insufficiency of wind-only schemes. The hybrid parameter R_H in equation (2) is evaluated in Figure 3 (b). We also demonstrate the transfer velocity predicted by R_H-based formulae from Brumer et al. (2017) and Woolf (2005) in panel (b). Woolf (2005) divided the full expression of gas transfer velocity into two parts which comprised non-breaking and wave breaking contributions. R_H was used for characterizing wave breaking contribution while the non-breaking formula was a linear function of wind friction velocity. We compute the total transfer velocity (green dashed line) here and the non-breaking contribution becomes insignificant under intermediate to high wind speeds (not presented). R_H in panel (b), to some extent, can collapse the data sets and the fitting result (black solid line) is close to the predictions of R_H-based formulae (blue and green dashed lines). However, these formulae capture relatively low variabilities. The computed r² are 25% and 27% for equations from Brumer et al. (2017) and Woolf (2005), respectively.

Figure 3

4 Discussion and Conclusion

The formulae (4) and (5) were initially developed in laboratory to account for the direct impacts of wave-related mechanisms on CO₂ gas transfer. They are fully dimensionless equations which are suitable for further application to open ocean conditions. With combined laboratory and field measurements, we adjusted and validated the formula (4). Its dominant term R_HM is built upon wave parameters and is a measure of wave induced turbulence (Babanin 2006). The effects of wind, on the other hand, is scaled as an enhancement factor (Ũ) to account for the marginal contributions from the wind. We also established equation (10) to emphasize effects of wave breaking through probability b_T and bubble’s contribution scaled with $R_{HM}^3$. The utilization of laboratory measurements provides the complement to diversified ocean wave conditions in the validation of our formulae. The formulae are also effective when they are evaluated with field observations. In Figure 4, we show an example of comparison between equation (4) and dimensional K₆₆₀ correlated with 10 m wind speed by using HIWINGS observations. The fitting coefficient r² for equation (4) in panel (a) is higher than that for using wind speed in panel (b) (0.92 vs. 0.87). The obtained r² for equation (10) is 0.89 which also indicates a better fitting than that in panel (b). In additon, similar conclusions can be achieved by using four field observations.

Figure 4

In this paper, b_T in equation (10) is estimated through wave measurements. But when it is not measured, parameterizations can be employed. b_T is determined by wave, wind and bottom-proximity properties. Banner et al. (2000) and Babanin et al. (2001) proposed dependences of b_T on significant wave steepness (ɛ) of the spectral peak in both deep and finite-depth water. Moreover, their studies show that the waves start to break (b_T > 0) when is higher than 0.055. The performance of equation (10) is limited under the condition of b_T = 0. Thus, equation (4) could be a supplement to equation (10) using the criteria of wave steepness,

The other motivation for proposing equation (11) is its potential application in modeling since the variables can be readily computed in spectral wave models (e.g., WAVEWATCH III (Rogers and Zieger 2014)).

The wave breaking induced air entrainment through bubbles may be more notable in open ocean than in laboratory wave tank because of the difference in wave scales and breaking strength. If the bubble’s effect becomes evident for CO₂ transfer in ocean, it is reasonable that equation (10) captures more variabilities compared with equation (4). Due to the paucity of bubble measurements in this study, the power dependence of bubble injection rate on R_HM could be further improved with more data. Additionally, bubble-mediated gas transfer is actually complicated, which is relevant with wave breaking strength (Manasseh et al. 2006), bubble injection depth, rise velocity and transfer rate for bubbles with different sizes (Woolf et al. 2007). Physical parameterizations for bubble-mediated transfer should be involved in the future study.

The results of DYNAMO project in our study have higher deviations than that of other projects, which could be attributed to low signal-to-noise ratio under mild wind and wave states. The observed partial pressure and fluxes for CO₂ are generally less significant than that of other cruises. Diurnal heat fluxes at daytime causing surface stratification and subsequent buoyancy fluxes at night (McGillis et al. 2004) are also important. These factors may all contribute to the deviation of predicted gas transfer velocity. CAPRICORN project was not specially designed for gas exchange experiment and only had one equipment for CO₂ flux. Thus, limited results are obtained after data processing. Our formulae only reflect the role of dynamic processes related with wind and wave, and can not represent other mechanisms that might significantly influence CO₂ gas fluxes, such as surfactants (Frew 1997) and rainfall (Takagaki and Komori 2007) which can suppress and intensify surface turbulence respectively. All the other influencing factors could lead to divergence of our formulae. In addition, the proposed parameterizations in this study are restricted to CO₂ or other sparingly soluble gases, rather than soluble gases. Air-side dynamics are the main constraint for soluble gas exchange, for which suitable parameterizations should be built upon wind parameters. Although Brumer et al. (2017) used a single parameter R_H (or R_B) and reconciled both CO₂ and DMS measurements, the obtained fitting may exhibit more discrepancies. The role of bubbles for gases with different solubilities needs to be properly represented in gas exchange scheme as well.

To summarize, CO₂ gas transfer models are proposed and validated through combined data from laboratory and field campaigns. Non-dimensional formulae (4), (10) and (11) depend on wave-related terms and can successfully bring together all data sets which span a wide range of wind and wave conditions. Bubble-mediated gas transfer is incorporated in formula (10), the performance of which is thus enhanced.

Data Accessibility statement

Data sets in this study are publicly available online. The laboratory data is available at Li et al. (2021); The field data from HIWINGS, CAPRICORN and DYNAMO campaigns can be found at ftp1.esrl.noaa.gov/psd3/cruises; The SOGASEX data can be accessed at http://www.bco-dmo.org/project/2064

Competing Interests

The authors have no competing interests to declare.