Impact of Nonzero Intercept Gas Transfer Velocity Parameterizations on Global and Regional Ocean–Atmosphere CO₂ Fluxes

Abstract

Carbon dioxide (CO₂) fluxes between the ocean and atmosphere (FCO₂) are commonly computed from differences between their partial pressures of CO₂ (ΔpCO₂) and the gas transfer velocity (k). Commonly used wind-based parameterizations for k imply a zero intercept, although in situ field data below 4 m s⁻¹ are scarce. Considering a global average wind speed over the ocean of 6.6 m s⁻¹, a nonzero intercept might have a significant impact on global FCO₂. Here, we present a database of 245 in situ measurements of k obtained with the floating chamber technique (Sniffle), 190 of which have wind speeds lower than 4 m s⁻¹. A quadratic parameterization with wind speed and a nonzero intercept resulted in the best fit for k. We further tested FCO₂ calculated with a different parameterization with a complementary pCO₂ observation-based product. Furthermore, we ran a simulation in a well-tested ocean model of intermediate complexity to test the implications of different gas transfer velocity parameterizations for the natural carbon cycle. The global ocean observation-based analysis suggests that ignoring a nonzero intercept results in an ocean-sink increase of 0.73 Gt C yr⁻¹. This corresponds to a 28% higher uptake of CO₂ compared with the flux calculated from a parameterization with a nonzero intercept. The differences in FCO₂ were higher in the case of low wind conditions and large ΔpCO₂ between the ocean and atmosphere. Such conditions occur frequently in the Tropics.

1. Introduction

The Earth’s ocean currently absorbs about 26% of all anthropogenic carbon emissions, with an anthropogenic flux of ca. 2.4 Gt C yr⁻¹ [1]. The ocean is a huge reservoir of CO₂, and a better understanding of how the ocean absorbs CO₂ is critical for predicting climate change. Carbon dioxide fluxes (FCO₂) between the ocean and atmosphere are normally calculated from the difference between oceanic and atmospheric partial pressure of CO₂ (pCO₂) and gas transfer velocity (k):

FCO₂ = kα ΔpCO₂,

(1)

where α is the solubility depending on temperature and salinity [2].

The k is typically calculated from parameterizations depending solely on wind speed, as reviewed by Wanninkhof [3]. However, these parameterizations lack in situ data for the low wind regimes, that is, wind speed <4 m s⁻¹ [4]. Therefore, they assume that for wind speed close to zero, the k is also zero.

The nonzero intercept is important for CO₂ budgets inferred from ocean pCO₂-based products, because there are large oceanic regions with prevailing low wind conditions (in Archer and Jacobson [5] and in Wurl, et al. [6]) (e.g., the Mediterranean Sea in summer and the tropical Pacific). Low wind regimes with limited surface mixing can lead to stratification of the upper layer and, given favorable conditions and assuming the nutrients are not limiting, to increased ocean productivity that would decrease the pCO₂ in the surface ocean and would increase the ΔpCO₂. With commonly used k parameterizations, these regions will have a net FCO₂ close to zero during periods of low sea states because the parameterization will predict k to be equal or close to zero. Consequently, the ocean’s source will be underestimated in low wind regime regions, such as the Tropical Ocean. It is important to quantify the associated bias as ocean pCO₂-based flux products offer one approach to estimating the global uptake of anthropogenic carbon [1,7,8].

Parameterizations with a nonzero intercept have been reported in [9,10,11,12,13,14,15], and the nonzero intercept ranges from 1 to 15 cm h⁻¹. Gas transfer velocity is unlikely to approach zero during low wind conditions as the wind speed is not the only main factor driving gas transfer across the surface; surface currents, buoyancy fluxes, micro breaking waves, and the presence of surface films (i.e., sea surface microlayer) affect gas transfer as well [3].

Zhang and Cai [16] reported potential biases based on an assumed nonzero intercept of 10 cm h⁻¹ reported in published papers [10,17,18,19] and pointed out that if such a nonzero intercept is adapted, the global CO₂ uptake for the year 1995 would increase by 0.2 Gt C yr⁻¹. However, they also pointed out that more in situ data are required to compute nonzero intercepts and to combine model and observational studies in order to draw more robust conclusions about relevant cases.

The purpose of this study is to quantify the uncertainty of global and regional carbon fluxes introduced by uncertain gas transfer velocity parameterizations. We first present empirical evidence for a nonzero intercept in gas transfer velocity parameterizations by presenting a database of 245 in situ measurements of FCO₂. We then use this database to evaluate the influence of the nonzero intercept on flux calculations. Hence, the evaluation of these parameterizations should not be seen as a way to promote our new parameterization as we and others have already determined that wind speed is not the only force that determines gas transfer velocities [9,20,21,22] and have identified other parameters, such as turbulent kinetic energy and surfactants, that affect the transfer. We believe the analysis offered by this paper is important because wind-based parameterizations are the most commonly used parameterizations, and because wind speed is an easier parameter to measure and/or is available from different sources. For these reasons, we evaluate the impact of the flux with and without an intercept on the gas transfer velocity parameterization and quantify the differences, that is, by assuming the flux is close to zero when wind speeds are close to zero. In a second step, we quantify the impact of nonzero intercept parameterizations on estimates of the oceanic uptake of carbon using a global pCO₂-based product. To calculate differences in global and regional FCO₂, five different gas transfer velocity parameterizations were used with the global maps of pCO₂ spanning 1982–2018; two commonly used parameterizations had a zero intercept, and three parameterizations included a nonzero intercept. In a third step, we tested the implications of a nonzero intercept in the framework of a global three-dimensional ocean carbon cycle model. We discuss the short-term and long-term evolution of the FCO₂ and the evolution of the natural carbon inventories in the ocean and atmosphere after introducing a nonzero intercept. We present both the global evolution of these variables and spatial patterns of change. The model is run for 5000 years, such that the system has time to approach new steady state conditions.

2. Materials and Methods

2.1. In Situ Flux Measurements with the Sniffle Floating Chamber

The field measurements were performed during several cruises (

Table 1. Cruise names, research vessels, dates, areas of observation, number of quality-controlled gas transfer observations, and studies that have discussed the cruises in more detail.

Cruise ID	Research Vessels	Start Date	End Date	Year	Area	# obs	Ref.
M117	Meteor	1 August	12 August	2015	Baltic Sea	19
FK161010	Falkor	12 October	6 November	2016	Timor Sea	60	[9]
HE491	Heincke	10 July	25 July	2017	Norwegian fjords & coastal North Atlantic	47	[9,20]
EMB184	Elisabeth Mann Borgese	1 June	10 June	2018	Baltic Sea	56
		31 March	2 August	2016	Jade Bay *	34	[23]
		3 April	4 September	2017	Jade Bay *	29	[25]

* Single daily samplings from different small boats, not a cruise.

). An autonomous drifting buoy, Sniffle, was deployed to measure the partial pressure of CO₂ (pCO₂) in the air, at a water depth of 1.2 m, and inside a floating chamber. The Sniffle is described in detail in Ribas-Ribas, et al. [23] and a diagram is presented in Figure 1.

In short, we used the floating chamber technique and measured the turbulence inside and outside the chamber simultaneously to see if turbulence affected the measurements. Furthermore, we monitored temperature, humidity, and pressure inside the floating chamber in order to detect any problems in the measurements. This procedure allowed us to overcome the problems associated with floating chamber techniques identified by previous research. Using this approach, we were able to measure the short temporal and spatial resolution needed to understand the rapid gas transfer velocity processes, in contrast to other techniques, such as eddy covariance or dual tracers. pCO₂ was determined using an infrared gas analyzer (OceanPack LI-COR LI-840X; SubCtech GmbH, Germany). The salinity, temperature, and wind speed were obtained from the research catamaran S³ [24] deployed simultaneously with the Sniffle. The CO₂ fluxes were calculated as follows:

$${\mathrm{F}}_{{\mathrm{CO}}_2}=\frac{{\mathrm{dpCO}}_2}{\mathrm{dt}}\frac{\mathrm{V}}{\mathrm{S}\mathrm{T}\mathrm{R}},$$

(2)

where $\frac{dpC{O}_2}{dt}$ is the slope of the pCO₂ change in the floating chamber, V is the volume of the floating chamber, S represents the surface area of the floating chamber, T represents the water temperature at a depth of 1 m from S³, and R is the gas constant. Measurements were excluded when the regression for the slope of the pCO₂ change in the floating chamber was R² < 0.90. The equation of the gas transfer velocity ${\mathrm{k}}_{\mathrm{w}}$ is

$${\mathrm{k}}_{\mathrm{w}}=\frac{{\mathrm{F}}_{{\mathrm{CO}}_2}}{\mathsf{\alpha }\left({\mathrm{pCO}}_{2\mathrm{water}}-{\mathrm{pCO}}_{2\mathrm{air}}\right)}$$

(3)

The solubility coefficient α depends on the temperature and salinity of the seawater and was calculated according to Weiss [2]. Finally, ${\mathrm{k}}_{\mathrm{w}}$ was standardized to k with the following formula:

$$\mathrm{k}={\mathrm{k}}_{\mathrm{w}}{(\frac{660}{\mathrm{Sc}})}^{0.5},$$

(4)

$Sc$ is the temperature-dependent Schmidt number Wanninkhof [3]. More details about the gas transfer velocities and CO₂ data from some of the cruises can be found in Mustaffa, et al. [9], Banko-Kubis, et al. [20], Ribas-Ribas, et al. [23], Stolle, et al. [25].

As a further quality control analysis for the floating chamber technique, we only used data with wind speeds lower than 7 m s⁻¹ (to avoid breaking-wave interference with the chamber), and we deleted outliers, defined as more than three standard deviations from the mean in each wind bin (2.4% of our data). A least squares quadratic regression with wind speed was used to determine the nonzero intercept parameterization that best represented the flux measurements. This new parameterization was compared with the zero intercept parameterization of Wanninkhof [3]. We used the ancillary data measured during the measurement campaigns detailed in Table 1. Accurate uncertainty estimations are difficult. Therefore, we assumed a maximum uncertainty of 10% for all parameters (worst scenario) (i.e., volume and surface of floating chamber, temperature and salinity, and pCO₂) affecting the calculation of k, as done before in [26] and [23]. We then propagated these errors on our calculations according to the method used by Taylor [27], resulting in 10.8% variation for CO₂ fluxes [23].

2.2. Global Surface Ocean pCO₂-Based Flux Product

We downloaded seawater and atmospheric pCO₂ from http://www.bgc-jena.mpg.de/CarboScope/?ID=oc to calculate ΔpCO₂. The pCO₂ observation-based product was obtained by fitting a simple data-driven diagnostic model of ocean mixed-layer biogeochemistry [28] to surface-ocean CO₂ partial pressure data from the Surface Ocean CO₂ Atlas (SOCAT) [29]. Auxiliary datasets for sea surface salinity were taken from the World Ocean Atlas 2013 v2 https://www.nodc.noaa.gov/OC5/woa13/ (climatological mean field without temporal variability). We obtained data on wind speed and sea surface temperature with six-hourly temporal resolution from the ERA-Interim data product [30], which is produced by the European Centre for Medium-Range Weather Forecasts (https://www.ecmwf.int/en/forecasts/datasets/archive-datasets/reanalysis-datasets/era-interim). We then interpolated ΔpCO₂ and salinity on the ERA-Interim grid. We carried out the FCO₂ calculation in MATLAB (MathWorks) from wind speed, sea surface temperature, salinity, and ΔpCO₂ at each ERA-Interim grid point using the “CO₂flux” function available from https://github.com/mvdh7/co2flux. We then normalized each grid point from Rödenbeck, et al. [31] dataset by ocean surface area. Finally, we integrated FCO₂ both globally and across each region of interest. We tested for five different parameterizations, as shown in

Table 2. List of gas transfer velocities parameterization (with equation, units, and references) used in Sniffle, the surface ocean observation–based method (SOCAT), and Bern3D evaluations.

Parameterization Name	Reference	$\mathbf{Equation}*{\left(\frac{\mathbf{S}\mathbf{c}}{\mathbf{660}}\right)}^{\mathbf{-}\mathbf{0.5}}\mathbf{[}\mathbf{cm}\mathbf{s}{}^{\mathbf{-}\mathbf{1}}]$	Used in:
Wanninkhof (W)	[3]	kW = 0.251 U102	Sniffle & SOCAT
This study		kM = 5.7 + 0.23 U102	Sniffle & SOCAT
Wanninkhof + Intercept (W+I)	[15]	kW+I = 10.7 + 0.30 U102	SOCAT
Krakauer (K)	[32]	kK = 2.275 * U10	SOCAT & Bern3D
Krakauer + Intercept (K+I)		kK+I = 11 + 2.275 * U10	SOCAT & Bern3D

.

Numerous parameterizations have been developed and reported, and we selected these five parametrizations for the two reasons: (i) The W is an update of a quadratic parameterization [33], commonly used since 1993 to estimate FCO₂. K has a shorter history and has not been used as commonly as other parameterizations, but it is implemented in the Bern3D model used in our study. K is an extensive data-assimilation effort based on ¹⁴C isotopes. We then added the nonzero intercept to these two commonly used parametrizations (W+I and K+I) to evaluate the direct difference on adding the intercept. Finally, we evaluate the parametrization obtained from this study. (ii) Computing time is expensive, which limited us to investigating only five parameterizations.

2.3. Bern3D Ocean Model

The Bern3D Ocean Model of Intermediate Complexity is a three-dimensional frictional geostrophic ocean model [34] with an isopycnal diffusion scheme and Gent–McWilliams parameterization for eddy-induced transport [35]. It is coupled to a sea ice component, a single-layer energy moisture balance model of the atmosphere [36], a prognostic marine biogeochemistry module [37,38], and a sediment module [37]. The horizontal resolution is 41 by 40 grid cells, and there are 32 logarithmically scaled vertical layers. All model components share the same horizontal resolution. A monthly wind stress climatology [39] is prescribed at the surface. Model evaluation is presented in [40,41,42].

Ocean-atmosphere gas exchange and carbonate chemistry were modeled based on the Ocean-Carbon Cycle Model Intercomparison Project (OCMIP)-2 protocols [43,44], with updates for the calculation of the Schmidt number [3] and carbonate chemistry [45]. The gas transfer was computed according to Equation (5). In contrast with the OCMIP-2 protocols, the gas transfer velocity k in Bern3D scales increases linearly with wind speed, following [32], and the global mean ocean-atmosphere transfer rate is reduced by 19% compared with OCMIP-2 to match observation-based estimates of natural and bomb-produced radiocarbon [46]. The k (in m s⁻¹ as implemented in the model) is therefore

$$\mathrm{k}=0.81\left(1-{\mathrm{F}}_{ice}\right)\ast 7.798\ast {10}^{-6}\ast wind{\left(\frac{\mathrm{Sc}}{660}\right)}^{-0.5},$$

(5)

where F_ice is the modeled sea ice concentration, wind is the prescribed climatology of monthly mean absolute wind speeds, and Sc is the Schmidt number following [3]. The normalization factor ensures the conservation of the globally integrated rate of annual mean gas transfer with respect to the original OCMIP formulation.

In a sensitivity experiment, a nonzero intercept is added to the gas transfer velocity function of CO₂, such that

$$k_{K+I}=k_K+3.06\ast {10}^{-5}[\mathrm{m}{s}^{-1}],$$

(6)

or, to unify with units

$$k_{K+I}=k_K+11[\mathrm{cm}{h}^{-1}],$$

(7)

Simulations start from preindustrial boundary conditions (A.D. 1765) with pCO₂ prescribed at 278 ppm; model drift is negligible. The nonzero intercept is added to the k at every grid cell starting from the beginning of a run, after which CO₂ is allowed to evolve freely. The simulation is run for 5000 years and changes in CO₂ feedback on simulated climate.

3. Results

3.1. Flux Measurements with Sniffle

3.1.1. Parameterization

Figure 2 shows a mean value of the gas transfer velocity calculated from Sniffle of each binned wind speed normalized at 10 m and the wind speed distribution during six expeditions (see Table 1). Gas transfer velocities generally increased with wind speed. The data show the usual scatter from field measurements (see, for example, Johnson [4]). The excluded outliers did not arise from measurement technique problems (as we did a first-step quality control analysis to detect lifting of chamber and error in measurements inside the chamber). The outliers were probably real data (but mathematically outliers) due to biogeochemistry or physics factors, such as chemical enhancements, water or air side convection, and turbulence. During the measurement campaigns, we encountered wind speeds in a range from 0.4 to 10.5 m s⁻¹, although we only used k values in which wind speed was lower than 7 m s⁻¹ to avoid data with breaking waves in the chamber. At low wind speeds (U₁₀ ≤ 4 m s⁻¹), the k values varied between 0.5 and 14.6 cm h⁻¹. At higher wind speeds (U₁₀ > 4 but ≤7 m s⁻¹), the k values varied between 0.9 and 32.8 cm h⁻¹. The black line is the model that best fits the field data, and it is k₆₆₀ = 5.71 (± 1.28) + 0.23 (±0.06) * U₁₀². The blue line is the Wanninkhof [3] parameterization with a zero intercept in the y-axis (zero flux when wind is zero).

Further evaluation of the relationship between field gas transfer velocity and other environmental factors (wind speed, surfactants, turbulence) has been evaluated in greater detail in other studies [9,20].

3.1.2. Comparison of Zero Intercept and Nonzero Intercept Parameterizations

We calculated FCO₂ using the parameterization proposed by Wanninkhof [3] (F_W) and the one from this study (F_M) detailed in Table 2 and obtained from the measurement campaigns (Table 1). Figure 3 shows the relationship between the FCO₂ calculated from Wanninkhof (F_W) and this study (F_M). F_W values ranged from −21 to 28 µmol m⁻² min⁻¹, whereas F_M values ranged from −29 to 43 µmol m⁻² min⁻¹. For our observations (i.e., U₁₀ < 7 m s⁻¹), positive FCO₂ were always higher with the new parameterization, taking the nonzero intercept into account (i.e., FCO₂ were above the 1:1 line), which means that source regions became stronger sources. On the other hand, negative FCO₂ were always lower with the new parameterization (i.e., FCO₂ are below the 1:1 line), which meant that sinks become stronger sinks. The conventional parameterization calculated an FCO₂ equal to 0 (ranging from −1.7 to 2.2 µmol m⁻² min⁻¹) for low wind conditions (U₁₀ < 2 m s⁻¹) (Figure 4b), and the modified parameterization with the nonzero intercept had a range from −12 to 20 µmol m⁻² min⁻¹ for low wind conditions (Figure 4e).

We then present the difference between F_W and F_M and calculate their ratio, as it indicates how different F_M is from F_W. When F_M/F_W is equal to one, the two FCO₂ are identical. If the F_M/F_W ratio is greater than 1, the gas transfer parameterization with a nonzero intercept leads to higher FCO₂ than the parameterization without an intercept.

The ratio between F_M and F_W ranged from 1.4 to 124 (Figure 5). Ratios lower than 2 corresponded to wind speeds higher than 4.6 m s⁻¹. All FCO₂ computed with the nonzero intercept parameterization (F_M) were higher than those computed with the conventional parameterization (F_W). The highest ratios (from 6.6 to 124) corresponded to wind speeds lower than or equal to 2 m s⁻¹. F_M/F_W ratios higher than 25 corresponded to wind speeds lower than 1.0 m s⁻¹, where F_W ranges from −0.1 to 0.2 µmol m⁻² min⁻¹ (Figure 4b) even if ΔpCO₂ ranges from –300 to 211 µatm (Figure 4a and Figure 5a). F_M exhibited a linear relationship with ΔpCO₂ (Figure 4d) and a scatter for all wind speed ranges (Figure 4e).

We compared the in situ fluxes measured directly with Sniffle with the fluxes we calculated with different parametrizations. These comparisons are shown in Figure 4c, f and

Table 3. Results of the comparison between direct measurements and calculated fluxes. r is Pearson’s correlation coefficient, RMSE is root mean square error, and R-squared is the coefficient of determination. The equation that better fits the relationship with in situ fluxes with slope P1 and intercept P2 is shown with 95% confidence bounds in parentheses.

	Comparison Statistics			FX = P1 * F Insitu + P2
	r	RMSE	R-Squared	P1	P2
FM	0.9122	6.9142	0.832	0.881 (0.831, 0.931)	2.928 (2.165, 3.691)
FW	0.8262	10.2064	0.682	0.398 (0.364, 0.432)	0.657 (0.133, 1.180)

. Overall, the slope was close to 1 for F_M (Table 3), where 1 represents perfect agreement. For F_W, the slope was 0.459 (Figure 4c and Table 3). R² and r were also higher for F_M than for F_W. F_M used a k calculated from in situ fluxes so the two values were not completely independent. However, the differences were higher than the 10.8% uncertainty estimated for our measurements that this bias could cause.

The analysis of the Sniffle data demonstrated that the use of a nonzero intercept parameterization might improve the calculated flux. In the following sections, we tested the implications in two different frameworks: the surface ocean observation-based method and the Bern3D model. Note that the three independent approaches used in this study use very different time spans (from minutes to thousands of years). We compared the different gas transfer velocity parameterizations within one approach, not between approaches, which meant that the time span differences were not problematic.

3.2. Application of Nonzero Intercept Parameterizations within Different Frameworks

3.2.1. Use and Description of Parameterizations

We evaluated surface ocean observations with five different parameterizations and Bern3D observations with two parameterizations (see Table 2 for details). We tested the original parameterization described by Wanninkhof [3] (W) and the parameterization from this study (M) with slightly lower slopes than W based on in situ Sniffle data. We then tested (i) the original parameterization described by Wanninkhof [3] (W), (ii) the parameterization from this study (M), (iii) the parameterization described by Wanninkhof [3] with the addition of a 10.7 cm h⁻¹ intercept (W+I) as reported in Ribas-Ribas, et al. [15] (iv) the original parameterization described by Krakauer, et al. [32] (K), and (v) the parameterization described by Krakauer, et al. [32] with the addition of an 11 cm h⁻¹ intercept (K+I) for the surface ocean observation-based method. Finally, we modified the Bern3D Ocean Model [34] with K and K+I. Zhang and Cai [16] used a nonzero intercept of 10 cm h⁻¹ from various studies. This is in agreement with the nonzero intercept of 10.7 cm h⁻¹ reported by Ribas-Ribas, et al. [15]_. This is the reason why we added this to W and K. For the Sniffle data analysis in Section 3.1, for clarity and brevity, we only compared the nonzero intercept parameterization from this study with the conventional parameterization from Wanninkhof. We used the five parameterizations associated with surface ocean observations presented in Table 2.

Figure 6 shows the distribution of gas transfer velocities versus wind speed, together with the wind speed histogram. We categorized the wind speed distribution into two regimes, in which different parameterizations (Table 2) produced similar k values, as follows:

• From 0–11 m s⁻¹, parametrizations with a nonzero intercept were similar, and parametrizations without an intercept were similar. This was the wind region for which the implementation of a nonzero intercept was relevant for global FCO₂.
• Above 11 m s⁻¹, the groups were quadratic parameterizations with higher slope versus linear or quadratic with lower slope. The nonzero intercept had no impact on FCO₂, because bubble-mediated transport may dominate the FCO₂. Therefore, this wind regime is not further discussed, and details on gas transfer at the high wind speed regimes can be found in Krall and Jähne [47], McNeil and D’Asaro [48].

3.2.2. Surface Ocean Observation-Based Method

Table 4. Annual mean ocean–atmosphere CO₂ fluxes (Gt C yr⁻¹) from the different regions (latitude and longitude given) and different parameterizations from the SOCAT-based globally gridded dataset analysis from 1982 to 2018. The final column presents the differences between the flux calculated from Wanninkhof’s (F_W) and this study’s (F_M) parameterizations.

	Areas	FW	FK	FM	FW+I	FK+I	FW-FM
Global		−2.56	−2.27	−1.83	−4.28	−3.52	−0.73
Southern	60° S – 80° S	−0.02	−0.02	−0.02	−0.03	−0.04	0.00
	180° W – 180° E
North Atlantic	65° N – 30° N	−0.18	−0.17	−0.15	−0.32	−0.28	−0.03
	55° W – 15° W
Tropical	14° N – 14° S	0.28	0.41	0.43	0.70	0.78	−0.16 *
	180° W – 180° E
South Atlantic	20° S – 40° S	−0.15	−0.15	−0.13	−0.27	−0.25	−0.01
	40° W – 10° E
South Pacific	20° S – 40° S	−0.34	−0.36	−0.31	−0.64	−0.59	−0.03
	155° W – 80° E

* Positive flux.

shows the annual mean ocean–atmosphere CO₂ fluxes on the global scale and for the different regions. Integrated globally, parameterizations with a nonzero intercept yielded larger carbon fluxes between the ocean and atmosphere compared with their zero intercept parameterizations. With a nonzero intercept, the W yielded an FCO₂ of −4.28 Gt C yr⁻¹; without the intercept, the exchange flux was −2.56 Gt C yr⁻¹ (Table 4). In general, the W+I parameterization gave higher absolute flux values except in the Tropical areas. W+I had larger k values than the other parameterizations throughout the entire wind speed range (Figure 5) except from 0 to 6 m s⁻¹, where K+I is higher. The K gave results similar to M except for the global ocean, where M was 0.4 Gt C yr⁻¹ lower than K (Table 4 and Figure S1). The Atlantic and Pacific (Table 4) act as a sink of CO₂, and the addition of the intercept increased the sink by a factor of 1.8. The Tropical Ocean was different: ΔpCO₂ was positive, so it acted as a source of CO₂ to the atmosphere for the entire evaluation period (Table 4). The nonzero intercept in M indicated that fluxes were 2.2 times larger than W for regions in the Tropical Ocean (Table 4). We further evaluated the Tropical Ocean as a “case study” in the next section.

3.2.3. Tropical Ocean

Figure 7 shows the different parameterizations used in this study (as in Figure 6) with the wind histogram for the Tropical Ocean. The wind distribution differed from that of the global ocean as it had a binomial pick located at wind speeds of 2–3 m s⁻¹ and a second pick similar to the one observed in the global distribution (Figure 6) around 6–7 m s⁻¹. If we multiply the k parameterizations with the wind speed histogram, we obtain the relative importance of the parameterizations normalized to the actual wind speed distribution. The higher frequency of low wind speeds gave the parameterization with a nonzero intercept a higher importance than for the lower 4 m s⁻¹ wind speed regime (see the pick in the blue square). M was close to the W+I and K+I (based on the sensitivity experiments). CO₂ fluxes calculated from M are almost double F_W, mostly due to the pick of wind frequency at lower wind speed. The parameterization M performed well, as it gave reasonable numbers in the global estimation and other parts of the ocean. W parameterization gave the lowest FCO₂—lower than K and M because of the lower importance of high wind speed. The FCO₂ calculated from K+I were the highest in the Tropics (Table 4), as it had larger k values from 0 to 6 m s⁻¹ and low wind speeds occurred most frequently in the Tropics.

3.2.4. Bern3D Ocean Model

We ran a sensitivity simulation in a well-tested ocean model of intermediate complexity to test the implications of different gas transfer velocity parameterizations for the natural carbon cycle in the absence of anthropogenic CO₂ emissions. We did not test the uncertainty of the oceanic sink for anthropogenic carbon. The importance of the gas transfer velocity for estimates of the uptake of anthropogenic carbon in ocean models has been assessed in previous studies [49,50,51], which have found that the sensitivity of the uptake to the gas exchange rate is small.

3.2.5. Global Evolution

The globally averaged FCO₂ fluctuated around zero for preindustrial conditions, implying an equilibrated carbon cycle (Figure 8a). In the first few years after adding a nonzero intercept to the k, the FCO₂ changed from a globally positive anomaly to a globally negative anomaly and back to globally positive anomaly (Figure 8a). After ~16 years, a globally positive FCO₂ anomaly established itself, and it takes ~1000 years for this FCO₂ anomaly to re-approach values of zero. Given sufficient equilibration time, the ocean approached a new equilibrium at which a zero net exchange of CO₂ between ocean and atmosphere reestablishes itself. The carbon inventories in the atmosphere and ocean changed according to the perturbed global FCO₂ and feedback phenomena triggered by sediment interactions. After a small decrease in atmospheric pCO₂ in the first few years, atmospheric CO₂ steadily increased by up to 8 ppmv, ~17 Gt C, after 5,000 years (Figure 8b).

A generally positive FCO₂ anomaly over most of the simulation period caused a net transfer of carbon from the ocean to the atmosphere. The global oceanic dissolved inorganic carbon (DIC) inventory slightly increased in the first few years and decreased by ~60 Gt C as the model approached a new equilibrium state (Figure 8c). The oceanic loss of carbon was amplified by sediment interactions, and the ocean lost additional carbon to the sediments (~40 Gt C) driven by the enhanced burial of calcium carbonate. This occurred because with smaller dissolved DIC concentrations, a larger area of sediments was exposed to higher oversaturation with respect to calcium carbonate, such that more calcium carbonate was preserved within sediments. With new carbon inventories in the ocean and atmosphere, a zero net exchange between the ocean and atmosphere reestablished itself. The increase in atmospheric CO₂ caused a slight warming (ΔT_atm = 0.125 K) with negligible impacts on the oceanic carbon cycle and circulation (not shown). The model did not fully equilibrate after 5000 years, as exchange fluxes with sediments had an equilibration timescale on the order of 50,000 years (see, for example, Roth, et al. [41]).

3.2.6. Regional Patterns

We showed regional patterns of the FCO₂ (Figure 9a) depicting the regional preindustrial steady state exchange of FCO₂ between the ocean and atmosphere. The model simulated CO₂ sunk in the gyre regions of all basins and in the North Atlantic and North Pacific. The Tropical Ocean and the Southern Ocean were source regions for atmospheric CO₂. As stated above, globally averaged, the FCO₂ was zero for preindustrial conditions.

The regional model response to increased gas transfer velocities was to amplify the general direction of preindustrial CO₂ exchange between ocean and atmosphere, that is, source regions became stronger source regions (more positive anomalies), and sink regions became stronger sink regions (more negative anomalies) (Figure 9b). These anomalies were strongest in the first few years after the perturbation (Figure 9b). At the end of the simulation, when the FCO₂ approached equilibrium values, there were still pronounced regional anomalies, with source regions generally becoming stronger sources and sink regions becoming stronger sinks (Figure 9c).

These anomalies caused changes in the carbon inventories of the ocean and atmosphere. The ocean globally loses 60 Gt C, with regional increases and decreases depending on the prevailing FCO₂ at the surface, subsequent water transport pathways, and amplifications through sediment interactions.

In the mid-latitude and North Atlantic (i.e., generally a CO₂ sink region), negative FCO₂ anomalies were simulated for the natural carbon cycle. This means that the region acted as a stronger CO₂ sink, which led to increased DIC concentrations in the ocean (Figure 9d). As the North Atlantic is a region of deep water formation, increased DIC anomalies spread in the Atlantic basin along the flows of deep water formation.

In the Southern Ocean, a CO₂ source region, positive FCO₂ anomalies were simulated. This meant that the region acted as a stronger CO₂ source, which led to decreased DIC concentrations (Figure 9d). As the Southern Ocean ventilates the deep South Atlantic and most of the deep Pacific, negative DIC anomalies spread along the flows of Southern Ocean deep water (Figure 9d). This caused higher calcium carbonate preservation and additional carbon loss to sediments.

The tropical oceans, which were simulated to become stronger source regions, led to lower DIC concentrations of the surface ocean. These anomalies did not spread in the deep ocean, as these regions are characterized by upwelling.

Regions of intermediate water formation, around 30° S, which became stronger CO₂ sinks, caused positive DIC anomalies in the upper ocean of the Atlantic and Pacific Oceans.

Despite small changes in the globally averaged exchange FCO₂, regionally large anomalies were simulated. As the model approached a new zero net exchange equilibrium state, more pronounced regional differences in positive and negative FCO₂ were simulated.

4. Discussion

We evaluated the impact on the CO₂ flux using gas transfer velocity parameterization with wind speed with and without an intercept, as wind speed is the easiest and most available parameter. By doing so, we are not proposing the use of a new wind-based parameterization, as we and others have already shown that wind speed is not the best stand-alone parameter for determining gas transfer velocities [9,20,21,22], and other parameters, such as turbulent kinetic energy and surfactants, will affect the transfer. However, we maintain that this analysis is important, as assuming fluxes are close to zero when wind speed is low introduces substantial uncertainty in regional CO₂ fluxes. Some of the disagreement found could be due to different spatial and temporal resolutions of observational data used for the empirical fits which make some of the most common used parametrizations appropriate for regional-to-global estimates of CO₂ fluxes, but these have also been used on smaller scales. Therefore, extrapolation from a local study to regional or global scales could be erroneous due to different wind speed variability and local water chemistry.

4.1. Global Ocean-Atmospher Carbon Fluxes Based on Observations of pCO₂

We compared our data-based estimates for the global ocean atmosphere exchange fluxes with those obtained by the Global Carbon Project (GCP). The GCP data-based products are based on the work of Rödenbeck, et al. [31] and Landschützer, et al. [52]. For the period 1982–2017, the mean oceanic sink was −1.85 Gt C yr⁻¹ according to Rödenbeck, et al. [31] and −1.42 Gt C yr⁻¹ according to Landschützer, et al. [52]. Both of the data-based products used here [33] were scaled to match the global average piston velocity. These data-based products reported in GCP yielded a CO₂ uptake close to the flux calculated with the parameterization proposed by this study (F_M) (i.e., −1.83 Gt C yr⁻¹). Another source of data-based products with which we can make direct comparisons is Takahashi, et al. [53]. Takahashi, et al. [53] reported a total global FCO₂ of −1.4 ± 0.7 Gt C yr⁻¹, which is lower than the data-based products reported in GCP by Rödenbeck, et al. [31] and F_M, probably because the Takahashi, et al. [53] database excluded observations made during the El Niño periods in the equatorial Pacific and those made in coastal zones. Takahashi, et al. [53] used a quadratic parameterization with a scaling factor of 0.26, which is close to W’s scaling factor of 0.251.

Globally, the F_W reveals a mean FCO₂ of −2.56 Gt-C yr⁻¹ for the period 1982–2018. The FCO₂ in this study using M, with a nonzero intercept, was −1.83 Gt-C yr⁻¹. Our results from Sniffle and the SOCAT-based globally gridded dataset show that including a nonzero intercept in wind-based parameterization will make a difference in the computed ocean-atmosphere CO₂ fluxes, as factors other than wind will not allow gas transfer velocities to approach zero at very calm sea conditions, such as convective mixing. However, the intercept cannot simply be added to existing parameterizations. By doing so, W+I and K+I reveal a global sink of ~−4 Gt C yr⁻¹, which is double the mean GCP estimate. For example, the M includes the intercept, but the slope is lower, and the lower slope leads to a good match with the mean GCP data-based estimate. On a global scale, the higher slope in the W compensates for the missing intercept. However, that does not hold for the tropical oceanic regions, which exhibit a frequent occurrence of low wind regimes; that is, the application of W to tropical regions introduces an error (which is discussed further below).

On the global scale, F_M reduced the oceanic uptake by 28% compared to the F_W. Even with the nonzero intercepts, M gave a noticeably lower FCO₂ than the K+I and W+I, due to the lower gas transfer velocities in the 5–10 m s⁻¹ range. However, the M computed higher k values at wind speeds <5 m s⁻¹, resulting in a better agreement of the FCO₂ with the W and K than with the W+I and K+I (Figure 7 and Table 4). Generally, K+I revealed higher FCO₂ than F_W. Even so, W computed the highest k values above ~12 m s⁻¹; this parameterization estimated lower k values in the low wind regimes (<4 m s⁻¹) than do the K and M.

We suggest that the behavior of parameterizations in the low wind speed regime is crucial in the computation of regional FCO₂ due to the relatively high frequency of the occurrence of wind speeds below 4 m s⁻¹ compared with wind speeds above 10 m s⁻¹ [5]. For the global analysis of the data product, the nonzero intercept parameterization led to a lower CO₂ uptake (Table 4).

It is interesting to note that W and K reveal similar results, despite the extensive debate in the literature about whether parameterizations should be based on linear, quadratic, or cubic relationships with wind [22]. In future work, it is recommended to run similar analysis with other parametrizations, such as with parameterizations developed for coastal waters.

4.2. Regional Patterns

FCO₂ estimations from the surface ocean observation-based method range in general from −0.02 to −0.04 Gt C yr⁻¹ (Table 4). This range is in close agreement with the mean annual flux of −0.06 Pg C yr⁻¹ in the ice-free zone of the Southern Ocean [53] but an order of magnitude lower than the FCO₂ estimate from McNeil, et al. [54]. Different parameterizations give different ranges, which means that the FCO₂ estimation will be stronger in sources or sinks depending on the season. Takahashi, et al. [53] also claimed that their flux was small due to the cancellation of the summer uptake of FCO₂ with the winter release of CO₂ caused by deep-water upwelling.

In the Tropical Ocean, wind speeds were more frequently below 5 m s⁻¹ compared with other global oceanic regions, with the consequence being that ignoring nonzero intercepts led to differences of 32% (Figure 7). Takahashi, et al. [53] reported that the equatorial Pacific is the major oceanic source of atmospheric CO₂ (+0.48 Pg C yr⁻¹), which agrees well with F_M for tropical oceanic regions (+0.43 Pg C yr⁻¹) (Table 4). From Sniffle and the regional results of the surface ocean observation-based method, we can conclude that the parameterizations without an intercept underestimate the FCO₂ in low wind speed regimes. For example, the difference between F_W and F_M in the tropics is 0.16 Gt C yr⁻¹, in agreement with the value of 0.2 Pg C yr⁻¹ reported by Zhang and Cai [16] with an estimated nonzero intercept of 10 cm h⁻¹. That means a 57% increase of FCO₂ in the tropical region (Table 4).

4.3. Inference of Nonzero Intercept Parameterizations Based on Observations

The nonzero intercept, reported here between 5 and 11 cm h⁻¹, may vary among different regions with changing biogeochemical forcing. For example, Frankignoulle, et al. [13] reported a nonzero intercept of 15 cm h⁻¹ in coral reef systems. On the other hand, some observations have indicated a smaller nonzero intercept (~1–3 cm h⁻¹) for nonequatorial regions [11,14,55]. Encountering low wind regimes during research cruises is challenging to plan for, and some techniques lack the sensitivity to record the slow and small change in gas transfer velocities under low wind conditions. In addition, the operations of research vessels are expensive, with the general aim of covering a large observational region within the shortest time frame possible. Therefore, future oceanic studies should be conducted in regions that offer a high probability of encountering low wind speeds and large ΔpCO₂ conditions. Furthermore, multidisciplinary research is required to advance the understanding of the mechanisms of processes driving gas exchange in low wind regimes. Crusius and Wanninkhof [56] found that transfer velocity and low wind speeds do not have a unique relationship, and they suggested that other processes, such as convective cooling [57], chemical enhancement [58], variations in fetch, the presence or absence of organic surface films [59], or rain [60], also affect the k. In particular, the buoyancy fluxes in low wind regimes can be important, as pointed out by MacIntyre, et al. [61]. For example, McGillis, et al. [10] suggested that their 8 cm h⁻¹ nonzero intercept was due to the buoyancy fluxes associated with a strong diurnal warming cycle. Buoyancy fluxes will affect surface turbulence, and they can be parameterized based on heat and momentum fluxes [62].

Zhang and Cai [16] estimated the ocean-atmosphere CO₂ flux by adding a constant of 10 cm h⁻¹ from previously published papers [10,17,18,19]. In agreement with our study, they reported that under low wind speed conditions, frequently found at low latitudes, wind speed may not even be the most dominant parameter for k. Zhang and Cai [16] also pointed out that not only the high ΔpCO₂ (due to upwelling) in the low latitudes but also the large oceanic coverage at low latitudes will significantly underestimate the CO₂ flux in the equatorial areas.

The observations with F_M/F_W ratios higher than 25 correspond to coastal stations (offshore Australia, Jade Bay, the Baltic Sea, and Norwegian fjords) (Figure 5). It is not reasonable to assume that the flux is close to zero when large ΔpCO₂ and low wind speeds are observed (Figure 4 and Figure 5) due to the presence other factors driving gas exchange, as discussed above. The mismatch between oceanic and coastal observation, based on our F_M/F_W ratios, could be one of the reasons why there has been extensive debate in the literature about whether coastal regions act as a source or sink of CO₂ [63]. Coastal systems are highly dynamic, and different conditions, such as upwelling, tidal influence, fresh water discharge, and fetch, are likely to be driving forces of gas transfer velocities. In addition, authors have claimed that specific parameterizations for gas transfer velocities should be used for coastal/estuarine regions [64]. For example, Borges, et al. [65] found that water currents contribute significantly to the k in the Scheldt estuary and that the k has a good linear relationship with wind speed and an intercept of 4 cm h⁻¹. Borges, et al. [12] also found a significant linear relationship between k and wind speed for the Randers Fjord and the Scheldt and Thames estuaries. The y-intercept varied from 1.2 to 9.7 cm h⁻¹, and these differences are related to varying contributions of tidal currents to water turbulence at the interface and fetch limitation.

The way the nonzero intercept is calculated or reported has varied between authors. Some authors have used the individual wind speed data [9] and others have binned the wind speed into different regimes [62,66,67]. Another difference is how the y-intercept is calculated. Most commonly, it is done by fitting the gas transfer velocities to the whole range of wind speed [3,9,12,60]. However, this fitting approach could be influenced by data in the higher wind regimes and, therefore, lead to an overestimation of the intercept. In order to reduce this artifact, Butterworth and Miller [11] determined the y-intercept by linear extrapolation of the two lowest wind speed bins.

4.4. Implications in Ocean Models

Although an accurate representation of the gas transfer velocity is crucial to estimating the current uptake of carbon from pCO₂-based observations [7,8,53], it has been shown to be of secondary importance for estimating the anthropogenic uptake of CO₂ in ocean models; the oceanic uptake of anthropogenic CO₂ in ocean models is driven primarily by increasing atmospheric CO₂ and is limited mainly by its transport from the surface to deep water [49,50,51]. ¹⁴CO₂ and ¹³CO₂, on the other hand, which have longer equilibration timescales with the surface ocean compared to CO₂ [68], are affected by different gas transfer velocities in models, and uncertainties in the gas transfer velocity are therefore critical for the analysis of carbon isotopes. Three-dimensional global ocean models have therefore been applied in combination with a range of suitable data-based metrics to infer the gas exchange velocity parameterization. Krakauer, et al. [32], for instance, used oceanic and atmospheric observations of ¹⁴C and ¹³C to constrain the gas exchange parameterization, utilizing the distinct signatures left by the ocean-atmosphere exchange of ¹⁴CO₂ and ¹³CO₂ in their model. Bern3D has adopted the relationship proposed by Krakauer, et al. [32] and further updated the relationship proposed by Müller, et al. [46] to best reproduce a suite of data-based metrics.

The simulations presented in this study assess how the natural carbon cycle in an ocean model adjusts to a perturbation of the gas transfer velocity. In response to a faster transfer velocity, source regions tend to become stronger sources and sink regions tend to become stronger sink regions. Transiently, this leads to a net transfer of carbon from the ocean to the atmosphere, but the system eventually reaches a new equilibrium at which the global sources and sinks balance each other again. Towards the end of the simulations, an equilibrated system at the global scale is reached. Regionally, however, stronger gradients are simulated between the source and sink regions.

5. Conclusions

It is clear that a nonzero intercept is crucial to computing CO₂ fluxes in low wind regimes based on observations of pCO₂, especially in the Tropics. Nevertheless, more systematic and multidisciplinary studies of in situ gas transfer velocities, together with investigations of buoyancy fluxes, surface films, and rain in different regions and conditions, are required to understand the mechanisms and variability of gas transfer across the ocean’s surface. Autonomous platforms and satellite missions have the potential to fill temporal and spatial gaps in future observations. Uncertainties in the gas transfer velocity are of secondary importance for estimating the anthropogenic uptake of carbon in ocean models.