The FuGas 2.5 Updated for the E ﬀ ects of Surface Turbulence on the Transfer Velocity of Gases at the Atmosphere–Ocean Interface

: Accurately estimating air–water gas exchanges requires considering other factors besides wind speed. These are particularly useful for coastal ocean applications, where the sea-state varies at ﬁne spatial and temporal resolutions. We upgrade FuGas 2.5 with improved formulations of the gas transfer velocity parametrized based on friction velocity, kinetic energy dissipation, roughness length, air-ﬂow conditions, drift current and wave ﬁeld. We then test the algorithm with ﬁeld survey data collected in the Baltic Sea during spring–summer of 2014 and 2015. Collapsing turbulence was observed when gravity waves were the roughness elements on the sea-surface, travelling at a speed identical to the wind. In such cases, the turbulence driven transfer velocities (from surface renewal and micro-scale wave breaking) could be reduced from ≈ 20 cm · h − 1 to ≤ 5 cm · h − 1 . However, when peak gravity waves were too ﬂat, they were presumably replaced by capillary-gravity waves as roughness elements. Then, a substantial increase in the turbulence and roughness length was observed, despite the low and moderate winds, leading to transfer velocities up to twice as large as those predicted by empirical u 10 -based formulations. and R.C.; Project administration, M.M. and F.L.; V.M.N.C.S.V.; draft, V.M.N.C.S.V.; Writing—Review & editing, V.M.N.C.S.V., M.M.,


Introduction
The dynamics of atmosphere-ocean gas exchanges are fundamental to the Earth's climate, because the ocean acts as both sink and source of greenhouse gases and dimethyl sulfide (DMS) to the atmosphere. It is generally assumed that the open ocean uptakes CO 2 from the atmosphere, despite the observed seasonal, inter-annual and regional variability [1][2][3]. On the other hand, the balances and fluxes of CO 2 , CH 4 , N 2 O and DMS across the coastal oceans' surface are heterogenic, due to processes such as upwelling, plankton productivity, and continental loads of organic matter and nutrients, from both natural and anthropogenic sources. Consequently, the coastal ocean can be a source of CO 2 [1][2][3][4][5][6][7][8][9][10], CH 4 [2,9,11,12], N 2 O [2,9,[13][14][15] and DMS [16] to the atmosphere, at least on a seasonal basis. The uncertainty about the atmosphere-coastal ocean exchange of CO 2 led scientists to question the overall CO 2 budget of the global ocean [3,17].

Air-Flow Regime
The u * can be estimated from several principles of fluid mechanics applied to the surface boundary layer and air-water interface. These require a preliminary definition of the air-flow regime over the sea-surface. Rough air-flow (or turbulent flow) is characterized by the air moving with strong lateral mixing between adjacent layers. The occurrence of eddies and swirls enhances exchanges perpendicular to the flow direction, i.e., vertically. The roughness length (z 0 ) and u * scale with this turbulence (hence z 0 and u * are positively correlated); z 0 being the theoretical height at which u * becomes zero (see Section 3.4 for a more detailed explanation on z 0 ). In smooth air-flow (or laminar flow) the air moves along parallel layers, with the least lateral mixing between adjacent layers. Due to the lack of flow separation, z 0 is expected to increase with decreasing u * [35,[55][56][57][58]96]. The air-flow regime is usually estimated from the roughness Reynolds number (R r = z 0 ·u * /ν a ), a dimensionless quantity measuring the ratio of inertial forces (z 0 ·u * in m 2 ·s −1 ) to viscous forces (υ a also in m 2 ·s −1 ). Exclusively turbulent rough flow has been proposed above R r = 1 [63], R r = 2 [96], R r = 2.3 [52,65,[97][98][99], R r = 2.33 [100], R r = 2.5 [101], or R r = 3 [102]. Below this threshold is a transient region between rough and smooth air flows [52,58]. Smooth air-flow has been proposed to occur below R r = 0.1 [63], R r = 0.11 [103], R r = 0.135 [65] or R r = 0.5 [100]. Its predominance over the sea-surface increases the shorter the fetch [52], and thus transient and smooth air flow conditions may be particularly important at the coastal ocean.

Friction Velocity
The friction velocity (u * ) at the air-water interface is the fluctuating component of the velocity of turbulent wind. Reworking the traditional formulation for wind stress τ = u * 2 ρ, the u * is a function of the shear stress (τ) exerted by air dragging over the sea-surface and of the fluid's density (ρ). Due to momentum conservation, the shear stress must be equal at both sides of the interface, the result being that u *a × ρ a 0.5 = u *w × ρ w 0.5 , where subscripts a and w stand for the air and water sides, respectively.
Therefore, u *w = u *a (ρ a /ρ w ) 0.5 [20,21,24,30,53,93]. Henceforth, u * refers to u *a . In this section are presented formulations for the estimation of u * , valid for any air-flow regime or exclusive to a specific air-flow regime.

Friction Velocities under Rough Air-Flow
The FuGas allows u * to be estimated from the uw and vw wind components measured by the EC method. One frequent option is Equation (2a) [26,52,66,67,95,100,104,105]. The u', v' and w' are the longitudinal, lateral and vertical wind fluctuations, with the overbars corresponding to the bin averaged second order crossed central moments. This formulation is implemented in the COARE 2.5. Because it may lead to biased estimates of individual u * measurements [64], Equation (2b) has been used instead [20,28,39,45,46,51,106,107] and implemented in the COARE 3.0 [64]. Because the estimates of u * from E-C measurements at different heights may show subtle differences [46,52,103,106], the u * measured at height z is converted to the u *s expected at the sea-surface by u *s = u * −0.0007z [52]. From then on, the former u * is replaced by the u *s value.
The Baltic Sea experiment provided a total of 1304 valid {uw,vw} observations, which, when applied to these equations, generally yielded resembling values of u * (Figure 1a,b). In only a few situations, mainly under low winds, did both formulations provide conspicuously different estimates. In only rare occasions, Equation (2a) yielded complex numbers with negligible imaginary components. For comparison among formulations, these u * were used together with the u 10n to obtain the observed drag coefficients (C D = u * 2 /u 10 2 ). The C D , often increasing with the lowest wind speeds (ranging from 1 m·s −1 to 3 m·s −1 ) (Figure 1c), matched previous observations done worldwide [35,46,[62][63][64][65]. This dynamic has been advanced as a consequence of capillary-gravity waves replacing peak gravity waves as roughness elements on the sea-surface [35,62], and/or the occurrence of cross or counter swell [46,50,51,54], and led to u * peaking at u 10n ≈ 2 m·s −1 .
The Baltic Sea experiment provided a total of 1304 valid {uw,vw} observations, which, when applied to these equations, generally yielded resembling values of u* (Figure 1a,b). In only a few situations, mainly under low winds, did both formulations provide conspicuously different estimates. In only rare occasions, Equation (2a) yielded complex numbers with negligible imaginary components. For comparison among formulations, these u* were used together with the u10n to obtain the observed drag coefficients (CD = u* 2 /u10 2 ). The CD, often increasing with the lowest wind speeds (ranging from 1 m•s −1 to 3 m•s −1 ) (Figure 1c), matched previous observations done worldwide [35,46,[62][63][64][65]. This dynamic has been advanced as a consequence of capillary-gravity waves replacing peak gravity waves as roughness elements on the sea-surface [35,62], and/or the occurrence of cross or counter swell [46,50,51,54], and led to u* peaking at u10n ≈ 2 m•s −1 .  Alternatively, the FuGas allows the u* to be estimated from the Drag Coefficient. In this case, the u* is estimated from u10, the definition of drag coefficient (also frequently mentioned in the literature as 'wind stress coefficient') reversed to u* = u10  CD 0.5 , and one among the several empirical formulations for the estimation of CD from u10, namely Equation (3a) [108], Equation (3b) [92], Equation (3c) [109], Equation (3d) [46], or Equation (3e) [110]. In the case of Equation (3b), its application clearly departed from the observations as well as all other formulations (Figure 1a,c), and  Alternatively, the FuGas allows the u * to be estimated from the Drag Coefficient. In this case, the u * is estimated from u 10 , the definition of drag coefficient (also frequently mentioned in the literature as 'wind stress coefficient') reversed to u * = u 10 × C D 0.5 , and one among the several empirical formulations for the estimation of C D from u 10 , namely Equation (3a) [108], Equation (3b) [92], Equation (3c) [109], Equation (3d) [46], or Equation (3e) [110]. In the case of Equation (3b), its application clearly departed from the observations as well as all other formulations (Figure 1a,c), and thus should be rejected. Mackay and Yeun [92] developed this C D parameterization simultaneously with a formulation that also widely overestimated transfer velocities (debated in Section 3.6). The C D and u * forecasted by the remaining formulations, although reasonably close to each other, (i) tended to slightly overestimate relative to the bulk of the observations, and (ii) were unable to account for the large C D and u * episodically occurring under low winds.
The FuGas also allows u * to be estimated from the wave field following the formulations developed by Gao et al. [111] for the coastal ocean [Equation (4a)] or offshore conditions [Equation (4b)]. These formulations use the wind measured at 10 m under atmospherically neutral conditions (u 10n ), and the celerity of peak gravity waves (c p ). Both formulations predict u * reasonably close to that predicted by the former C D -based formulations ( Figure 1d).

Friction Velocities Estimated from the WLLP (under Rough Air-Flow)
In a FuGas exclusive [73][74][75], u * was estimated from the wind measured at height z (u z ) by solving the Wind Log-Linear Profile (WLLP) backwards [Equation (5)]. This solution congregates the effects of wind speed, wave state and atmospheric stability, and thus is central to the most comprehensive physically-based numerical scheme provided in FuGas. The WLLP is the air-side complementary of the Velocity Defect Law for the aqueous boundary layer [26,27], the boundary condition being that at height z = 0 the u z −u s = 0. Thus, theoretically, the collinear surface velocity u s should not be neglected [63,106]. These laws report to turbulent flow. In the WLLP's case, that is the rough air-flow regime.
The estimation of u s from the sea-state is presented in Section 3.3. The sea-state effect over the aerodynamic roughness length (z 0 ) is presented in Section 3.4. The atmospheric stability function (ψ m ) was presented in previous FuGas applications [73][74][75], but new parameterizations derived from data specific to this Baltic Sea site [45] were added to FuGas 2.5. The von Kármán constant (κ) is 0.4, although 0.41 is also frequently used. As the WLLP must be iteratively solved for convergence [63,64,74,75], a preliminary guess for u * was required. This could be taken from one of the former solutions: E-C data, C D -based formulations or wave-based formulations. Only afterwards could the definite u * be estimated from Equation (5).

Friction Velocities upon Smooth and Transient Air-Flow Regimes
Given that under smooth-flow, u * and z 0 are negatively correlated, the u * may be estimated as u * = R r ·ν a /z 0 . This is the principle also applied to the z 0 estimation under smooth-flow [52,[63][64][65]112], being usually applied the same R r of the smooth-flow threshold (see Section 3.1). Under transient-flow, u * is a compromise between rough-flow and smooth-flow formulations. Because the u * estimation under smooth-flow requires other properties (such as z 0 ) that need be estimated ahead, it will be addressed later in Section 3.5.

Surface Velocity
The wind at 0 heights (u 0 ≈ 0) over the sea-surface has a movement relative to the moving surface water [26,27,52,62,63,106,[113][114][115][116]. The surface water velocity (u s ) was shown to affect the wind stress over the sea-surface, even in the open ocean, most often reducing it, and with an impact on the transfer velocities [27,[114][115][116]. Hence, FuGas provides four options for its estimation: (i) Disregarding the velocity of the sea-surface, in which case u s = 0. (ii) Using the collinear component of the water current; this being most suited for estuaries and tidal lagoon systems. In this case u s = w × cosθ, where w is the current speed and θ is the angle between the wind and current directions. (iii) Using the collinear component of the sea-surface drift currents induced by waves, estimated adapting Stokes transport law [117]. In this case u s = 4π 2 ·f p ·H s 2 /L p ·cosθ. Note that this is the average horizontal velocity, i.e., disregarding both the vertical movement and the back-and-forth horizontal movement of particles on the sea-surface as waves pass by.
During the Baltic Sea experiment the surface water velocity was always very low, and only on rare occasions the estimates from available formulations slightly exceeded 10% of the wind velocity ( Figure 2a). This is in accordance with the bulk of the observations reported in previous works. According to Wu's formulation, the wind-induced drift currents are a fundamental component of u s (Figure 2b). Nevertheless, u s has little relevance for the estimation of u * from Equation (5), and is particularly irrelevant under low winds (Figure 2c). The wind-and wave-induced drift currents estimated from Wu [62,114] are related with the wave age (cp/u10), where cp is the celerity of peak waves. (c) The resulting u* is compared with the estimated neglecting the us.

Roughness Length
Both the Law of the Wall (LOW) and the Velocity Defect Law assume that in contact with the surface there is a sub-layer with laminar flow. Within this sublayer, stress is predominately carried by viscous transfer and, although there may be an effective wind velocity (i.e., uz > 0 with z ≈ 0), the turbulent velocity tends to be zero (i.e., u* = 0). The depth of this sub-layer corresponds to the roughness length z0, and may tend to be infinitesimal depending on air-flow regime and roughness elements over the sea-surface. Upon rough-flow the roughness elements protrude the laminar flow sub-layer, i.e., z0 is smaller than the height of the roughness elements, causing the air-flow to separate from the sea-surface. In this case, z0 increases with increasing u*. Upon transient-flow the roughness elements have approximately the same height as in the laminar flow sub-layer. Upon smooth-flow the roughness elements stay within the laminar flow sub-layer, i.e., z0 is larger than the height of the roughness elements, causing the air-flow to not separate from the sea-surface. In this case, z0 increases with decreasing u*. Either peak gravity waves or shorter waves, even on the capillary-gravity range, can be the roughness elements on the sea-surface. Wu [100,113] proposed that the roughness elements are the waves travelling sufficiently slower than the wind for wakes to generate in their front, and the flow to separate.

Roughness Length under Rough Air-Flow
Under rough (turbulent) air-flow, the roughness length (z0) scales with turbulence, and thus with u*. The bulk of the empirical formulations assume that the roughness length depends exclusively on gravity waves as the roughness elements on the sea-surface. The simplest of them, by Charnock [118], proposes a z0 (units in m) dependency on u* [Equation (7) and Figure 3]. Here, g is the gravitational acceleration constant, and Charnock's coefficient is 0.01 < αCh < 0.02 over water surfaces [52,53,64,103,112]. The wind-and wave-induced drift currents estimated from Wu [62,114] are related with the wave age (c p /u 10 ), where c p is the celerity of peak waves. (c) The resulting u * is compared with the estimated neglecting the u s .

Roughness Length
Both the Law of the Wall (LOW) and the Velocity Defect Law assume that in contact with the surface there is a sub-layer with laminar flow. Within this sublayer, stress is predominately carried by viscous transfer and, although there may be an effective wind velocity (i.e., u z > 0 with z ≈ 0), the turbulent velocity tends to be zero (i.e., u * = 0). The depth of this sub-layer corresponds to the roughness length z 0 , and may tend to be infinitesimal depending on air-flow regime and roughness elements over the sea-surface. Upon rough-flow the roughness elements protrude the laminar flow sub-layer, i.e., z 0 is smaller than the height of the roughness elements, causing the air-flow to separate from the sea-surface. In this case, z 0 increases with increasing u * . Upon transient-flow the roughness elements have approximately the same height as in the laminar flow sub-layer. Upon smooth-flow the roughness elements stay within the laminar flow sub-layer, i.e., z 0 is larger than the height of the roughness elements, causing the air-flow to not separate from the sea-surface. In this case, z 0 increases with decreasing u * . Either peak gravity waves or shorter waves, even on the capillary-gravity range, can be the roughness elements on the sea-surface. Wu [100,113] proposed that the roughness elements are the waves travelling sufficiently slower than the wind for wakes to generate in their front, and the flow to separate.

Roughness Length under Rough Air-Flow
Under rough (turbulent) air-flow, the roughness length (z 0 ) scales with turbulence, and thus with u * . The bulk of the empirical formulations assume that the roughness length depends exclusively on gravity waves as the roughness elements on the sea-surface. The simplest of them, by Charnock [118], proposes a z 0 (units in m) dependency on u * [Equation (7) and Figure 3]. Here, g is the gravitational acceleration constant, and Charnock's coefficient is 0.01 < α Ch < 0.02 over water surfaces [52,53,64,103,112]. Because Charnock's [118] formulation could not reflect the influence of the wave field, several authors proposed updating this classical formulation with a αCh dependency on the wave age (cp/u*), as given by Equation (8) [97][98][99]107,[119][120][121]. The wave age is given by the celerity of peak waves (cp) relative to u*. Usually, 1 < cp/u* < 1000. The larger the ratio, the older the waves. Different ACh and BCh constants were proposed by the respective authors. Their values usually mean z0 increasing with younger waves [98,99,107,[119][120][121] Representing wave age by cp/u10 instead facilitates the identification of swell, as in fully developed (mature) wind-seas both velocities should be approximate [97]. Swell corresponds to older waves travelling faster than the wind, and in this case aged > 1, obviously aside from waves traveling in other directions. Gao et al. [111] developed a z0 parameterization relying exclusively on this wave age, with a calibration for the coastal ocean [Equation (9a)] and another for offshore conditions [Equation (9b)]. Both formulations showed a z0 generally scaling with u*. Nevertheless, matching the smooth flow theory, under the lowest u* and z0, typical of smooth air-flows, the inverse scaling took place, i.e., z0 increased with decreasing u* and Rr ( Figure 3). The off-shore formulation fit particularly well the Baltic Sea data: under the lowest u* it predicted z0 matching exactly that expected for smooth air-flows. The larger z0 variability under similar u*, but given different wave ages, demonstrated an enhanced ability of this formulation to adapt to local conditions. Because Charnock's [118] formulation could not reflect the influence of the wave field, several authors proposed updating this classical formulation with a α Ch dependency on the wave age (c p /u * ), as given by Equation (8) [97][98][99]107,[119][120][121]. The wave age is given by the celerity of peak waves (c p ) relative to u * . Usually, 1 < c p /u * < 1000. The larger the ratio, the older the waves. Different A Ch and B Ch constants were proposed by the respective authors. Their values usually mean z 0 increasing with younger waves [98,99,107,[119][120][121] Representing wave age by c p /u 10 instead facilitates the identification of swell, as in fully developed (mature) wind-seas both velocities should be approximate [97]. Swell corresponds to older waves travelling faster than the wind, and in this case aged > 1, obviously aside from waves traveling in other directions. Gao et al. [111] developed a z 0 parameterization relying exclusively on this wave age, with a calibration for the coastal ocean [Equation (9a)] and another for offshore conditions [Equation (9b)]. Both formulations showed a z 0 generally scaling with u * . Nevertheless, matching the smooth flow theory, under the lowest u * and z 0 , typical of smooth air-flows, the inverse scaling took place, i.e., z 0 increased with decreasing u * and R r (Figure 3). The off-shore formulation fit particularly well the Baltic Sea data: under the lowest u * it predicted z 0 matching exactly that expected for smooth air-flows. The larger z 0 variability under similar u * , but given different wave ages, demonstrated an enhanced ability of this formulation to adapt to local conditions.
Anctil and Donelan [106] added extra factors to the z 0 dependency on wave age. First, they found that the dimensionless roughness length scaled with the inverse wave age. Solved for z 0 , it leads to Equation (10a). This formulation also showed a z 0 generally scaling with u * , but inverting the scaling under conditions corresponding to smooth air-flows ( Figure 3). This formulation also fit well the Baltic Sea data sampled under the lowest u * , predicting z 0 matching the expected for smooth air-flows. This formulation also showed an enhanced ability to adapt to local conditions, demonstrated by the larger variability of z 0 estimated under similar u * , but given different wave ages. Anctil and Donelan [106] found that the dimensionless roughness length scaled equally well with the mean squared slope of the wave field (θ S , but also named <S 2 > [20,27]). Solved for z 0 , it leads to Equation (10b). Finally, they merged both factors in a multiple regression, which, solved for z 0 , leads to Equation (10c). Our data did not include measurements of θ S necessary for the implementation of Equation (10b,c).
Taylor and Yelland [110] (following Hsu [122]) proposed a z 0 dependency on the wave slope as estimated from the peak wave length (L p ) and significant wave height (H s ) [Equation (11a)], and calibrated by A w = 1200 and B w = 4.5. This formulation also showed a z 0 generally scaling with u * , but inverting the scaling under conditions corresponding to smooth air-flows ( Figure 3); it also fit well the Baltic Sea data sampled under the lowest u * , predicting z 0 matching that expected for smooth air-flows. This formulation also showed an enhanced ability to adapt to local conditions, demonstrated by the larger variability of z 0 estimated under similar u * , but given different wave ages. Taylor and Yelland [110] (following Donelan [103] and Anctil and Donelan [106]) presented an alternative formulation for a z 0 dependency on, simultaneously, the wave age and slope [Equation (11b)]. However, this formulation was subsidiary, with its authors not attributing a value for C w . Preliminary tests showed that the B w and C w terms could easily become antagonistic. We tested the calibration A w = 1200, B w = 1.5 and C w = 3.5. This calibration lost the good fit to the transient-smooth air-flow conditions, bringing back the behavior more characteristic of the Charnock-based formulations ( Figure 3).
Pan et al. [123] went back to the z 0 dependency on the inverse wave age and size, proposing a new three-parameter formulation where the dimensionless roughness length scales with the u * -based inverse wave age, i.e., z 0 /H s u * /c p . We reformulated it into Equation (12a). In a secondary formulation, they also revived the α Ch u * /c p . We reformulated it into Equation (12b), which is basically Equation (8) with new coefficients. These formulations brought back the behavior more characteristic of the Charnock-based formulations (Figure 3).
Hwang [124] proposed that the vertical scales of properties of the air-water interface, like the drag coefficient or z 0 , should depend on the (horizontal) wave length. Consequently, both the drag coefficient and z 0 increase non-linearly with wave age for very young waves, while decreasing non-linearly with wave age for very old waves. Hwang [124] tuned his formulation for deep water conditions. We adapted it for any water depth [Equation (13)]. In this case, ω * is the inverse wave age. This formulation uses the same parameters (wave age and length) as COARE [64] z 0 formulation adapted from Oost et al. [125].
The Baltic Sea data only comprehended the range of older waves, for which Hwang [122] predicted that z 0 decreases with wave age. Our data did not comprehend the extremely young waves for which Hwang [122] proposed a z 0 increasing with wave age. For waves with different L p to be similarly old, the faster travelling, longer waves occurred under higher winds, and attained higher heights. Therefore, although here was used Lp, other metrics such as Hs or the peak wave frequency (fp) worked similarly. Due to this similitude, within the range of older waves, our adaptation of the formulation by Hwang [122] behaved reasonably similarly to our adaptation of the formulation by Pan et al. [123] (Figure 3), and thus this formulation also showed a behavior more characteristic of the Charnock-based formulations.

Roughness Length under Smooth Air-Flows
Contrary to the rough-flow regime, upon smooth-flow, the roughness length (z 0 ) increases with decreasing u * , due to the lack of flow separation [35,[62][63][64][65]112]. The most usual parametrization is z 0 = R r × ν a /u * , with the R r value given by the smooth-flow threshold, namely R r = 0.11 [52,64,103,112] or R r = 0.135 [65] (Figure 3). An alternative formulation with unconstrained R r was advanced by Wu [35,62]. It proposes that at u 10 ≈ 5 m·s −1 actually occurs the smoothest air-flows over the sea-surface, below which the sea-surface gets rougher with decreasing winds, because capillary waves (whose restoring force is the surface tension σ) replace the mild gravity waves (whose restoring force is the acceleration constant g) as roughness elements. This shift leads to both drag and z 0 increasing with decreasing u * . However, once the shift is done, the drag increase can be so intense as to turn the air-flow rough again, even upon light breezes ( Figure 3). Wu's algorithm [Equation (14a)] has the water surface tension (σ) governing the wind-wave interactions. The surface tension (N·m −1 ) of pure water (σ 0 ) was estimated from Equation (14b), accounting for the water temperature (T w in K) and considering a critical water temperature T c = 647.1K [126]. Rectification for salt water surface tension (σ w ) was done by Equation (14c), using the surface salinity (S in ppt) [127].

Definite Estimation of Roughness Length
Theory says that under typical smooth flows, z 0 should be exclusively set by the smooth-flow term, under typical rough flows z 0 should be exclusively set by the rough-flow term, and only under the transient regime should it be a function of both. Yet COARE [64] always adds both terms. Since when one term is maximized the other is minimized-tending to infinitesimal-always adding them yielded results almost identical to using exclusively the term adequate for each regime (Figure 4a-c). Besides, it is computationally simpler and more efficient. Choosing the maximum between each term also yielded almost identical results. Another option was averaging between both terms. Using the arithmetic mean yielded slightly different results (Figure 4d-f). The fundamental difference was that z 0 was estimated slightly smaller. However, using the geometric mean yielded widely different results as the rough-flow infinitesimal terms became preponderant. Consequently, with low u * the definite z 0 also becomes infinitesimal (Figure 4g-i), contradicting general theory on laminar flows. Geometric means never yielding infinitesimal z 0 could only be obtained if the rough-flow term also never yield infinitesimal z 0 , which was only satisfied by the formulations by Actil and Donelan [106] [Equation (10a)], Taylor and Yelland [110] [Equation (11b)] and Gao et al. [111] [Equation (9a)]. Using an analogous to the double resistance estimator (i.e., 1/z 0 = 1/z rough + 1/z smooth ) yielded results similar to those when using the geometric mean. The smooth-flow term can be estimated from a fixed-R r formulation or from the formulation by Wu [35,62]. Either application yielded widely different results ( Figure 3). Since z 0 cannot be directly measured, it was impossible to determine which numerical scheme provided the best estimates.   (7) and (8) Alternatively, Wu [100] also advanced a general formulation for z0 under any air-flow regime by combining the formulations related with rough and smooth air-flows [Equation (15)]. The υw is the kinematic viscosity of water and β is a constant. Wu proposed αCh = 0.0185 and 2 < β< 2.5. These β values implied that gravity waves dominate the determination of z0. Oddly, the smooth-flow component of Equation (15) was inversed relative to the previous smooth-flow formulations. The resulting z0 was far off the predicted values given by other algorithms (Figure 3), and followed a narrow logarithmic trend that was basically dependent on u*, while the kinematic viscosity and surface tension of water were irrelevant.   (7) and (8) [35,62]. Roughness Reynold number estimated from R r = z 0 × u * /υ a .
Alternatively, Wu [100] also advanced a general formulation for z 0 under any air-flow regime by combining the formulations related with rough and smooth air-flows [Equation (15)]. The υ w is the kinematic viscosity of water and β is a constant. Wu proposed α Ch = 0.0185 and 2 < β< 2.5. These β values implied that gravity waves dominate the determination of z 0 . Oddly, the smooth-flow component of Equation (15) was inversed relative to the previous smooth-flow formulations. The resulting z 0 was far off the predicted values given by other algorithms (Figure 3), and followed a narrow logarithmic trend that was basically dependent on u * , while the kinematic viscosity and surface tension of water were irrelevant.

Converging Friction Velocities
One of the purposes of FuGas is to be applied to atmospheric and oceanic L4 data with large spatial and temporal coverage. In such cases, since direct observations are not available, u * must alternatively be estimated from empirical formulations. The FuGas user may opt for one simple parameterization, or its comprehensive iterative estimation involving the wind log-linear profile: first, a preliminary u * is estimated, and then z 0 and u * must be solved for convergence. In both COARE [63,64] and previous FuGas algorithms [74][75][76], convergence was generally attained after three iterations. Although the new FuGas schemes could be substantially more complex, convergence was still attained within three iterations, at most, and often with only one or two. Once convergence was achieved, we compared the forecasted u * with the u * estimated from direct observation obtained through the EC method. Yet, there was no guarantee that the EC-based u * reflected the true u * , particularly under the low winds, when turbulence is reduced and EC method is prone to failure, as previously demonstrated for this Baltic Sea station [7].
The estimation of the preliminary u * from a C D -based formulation was an important step. Following the Smith [108] formulation in Equation (3a) leads to the best fits (Root Mean Square Deviation RMSD ≈ 0.0348 m·s −1 ). The worst estimation was from the formulation by Mackay and Yeun [92] given in Equation (3b) (RMSD ≈ 0.155 m·s −1 ). The z 0 estimation was also determinant. Its smooth-flow term was better estimated from a fixed-R r formulation (RMSD ≈ 0.0343 m·s −1 ) than from Wu's [35,62] smooth-flow formulation given in Equation (14) (RMSD ≈ 0.0376 m·s −1 ) ( Figure 5). The z 0' s rough-flow term was best estimated from the Pan et al. [123] formulation in Equation (12) or the Gao et al. [111] formulation for the coastal ocean in Equation (9a) (RMSD ≈ 0.0343 m·s −1 ), and was worst estimated from the Toba et al. [97] formulation in Equations (7) and (8), or Gao et al. [111] formulation for the open ocean in Equation (9b) (RMSD ≈ 0.0355 m·s −1 ) ( Figure 5). The rough-flow and smooth-flow terms need to be merged under transient air-flows. The tests were not conclusive about which z 0 integration method should be preferred. Yet, estimating their geometric average generally provided the best fits (RMSD ≈ 0.0343 m·s −1 ). On the opposite extreme, integrating z 0 rough-flow and smooth-flow terms in Equation (15)  for the open ocean in Equation (9b) (RMSD ≈ 0.0355 m•s ) ( Figure 5). The rough-flow and smooth-flow terms need to be merged under transient air-flows. The tests were not conclusive about which z0 integration method should be preferred. Yet, estimating their geometric average generally provided the best fits (RMSD ≈ 0.0343 m•s −1 ). On the opposite extreme, integrating z0 rough-flow and smooth-flow terms in Equation (15) [Wu 100] led to largely overestimated u* (RMSD ≈ 0.141 m•s −1 ). The choice of critical Rr for smooth-flow or rough-flow thresholds was of little relevance, and so was the choice of surface velocity from drift currents. Colored marks represent the roughness length under rough air-flow estimated following, Charnock [118] [Chk55, Equation (7)], Gao et al. [111] formulations for the coastal ocean [G09Co, Equation (9a)] or off-shore [G09Of, Equation (9b)], Pan et al. [123] [PanZC, Equation (12b)] or Toba et al. [97] [Tob90, Equations (7) and (8]). Roughness length under smooth air-flow was estimated following Wu [35,62] [Equation (14)] or formulations relying on a fixed Rr. Smooth and rough air-flow terms under transient air-flow regimes were added, arithmetically averaged or geometrically averaged.
When the numerical scheme was fine-tuned, the empirical estimation of u * was generally close to its estimation from direct observation obtained through the EC method ( Figure 5), although often slightly overestimated. Most of this overestimation occurred under lower wind speeds, when the air-flow should be transient or smooth, and the EC-based estimations may be biased. Melville [129] proposed that the shift from rough to transient air-flows occurs at u * from 0.15 to 0.3 m·s −1 . Wu [35] proposed this threshold at u 10 ≈ 5 m·s −1 , which in our data corresponded to u * predominantly between 0.05 and 0.2 m·s −1 , and median at ≈ 0.13 m·s −1 . Below this transient region, the empirically estimated u * was larger than that estimated from direct observation ( Figure 5). The lack of precision of EC methods under such low winds was already debated. The present observations match previous observations and models where, under low winds, the u * decreases less sharply associated with decreasing wind velocity [58], while drag and transfer coefficients increase [35,[62][63][64][65]98]. To justify such an effect, Wu [35] proposed that the sea surface is aerodynamically rough even under light winds, on the account of capillary-gravity waves replacing gravity waves as roughness elements on the sea-surface.

The k wind Term
The experiment provided a total of 163 observations with valid transfer velocities, of which only 94 included the wave properties necessary for simulations using the most advanced algorithms available in the FuGas. Still, starting with the simplest, i.e., the empirical u 10 -based formulations, the formulations by Jacobs et al. [83] and Wanninkhof and McGillis [84] estimated respectively the largest and the lowest transfer velocities under the observed wind range (Figure 6a). The formulation by Wanninkhof [18] is the standard in Earth-System Models, and roughly represents the average behavior of u 10 -based formulations. Their application should use the wind speed at a 10 m height under atmospherically neutral conditions (u 10 n). If that is not the case during field experiment, the equivalent u 10 n must be estimated from the data. These standard transformations, although derived from the wind log-linear profile, neglect the surface velocity, and that z 0 is variable. However, the FuGas numerical scheme allows accounting for them. Wanninkhof's [18] formulation applied to this u 10 n yields RMSD = 9.42 cm·h −1 (Figure 6a). All u 10 -based formulations largely misfit the observed transfer velocities under the lowest winds and swell. This had already been demonstrated for the Östergarnsholm setup [7]. Several explanations have been advanced for the swell interference in the transfer velocities. The older and larger waves have been demonstrated as having a direct interaction with the atmosphere, and as interfering in the properties of the shorter wind-waves that often regulate the turbulence and air-flow over the sea-surface [77][78][79]. The interference of swell was demonstrated to be stronger with lower winds, steeper waves and/or when opposing the wind direction [50,51,54], which was precisely our case. Disregarding the wave-state has already been shown to introduce systematic bias to the calibration of bulk flux formulations for the transfer of heat and moisture from direct observation obtained through EC method, particularly under near-collapsed turbulence [48]. and relation with (u10n) wind speed (cp/u10n) and wave age. Waves aged cp/u10n > 1 correspond to overdeveloped sea. The EC transfer velocities were compared with (a) those estimated from the u10based formulations by Jacobs et al. [83] (Jac99), Wanninkhof [18] (Wan92) and Wanninkhof and McGillis [84] (WMG99), and (b-f) with comprehensive physically-based algorithms (u*-based) by Esters et al. [30]. [Equations (17) and (18)]. (a,b) k600 estimated from Esters et al. [30] algorithm and its relation with u*. (d-f) validation of k600 estimated from from Esters et al. [30] algorithm against k600 estimated from direct observation obtained through EC method.
Turbulence-driven transfer velocities should be more accurately inferred from the wind effectively exerting drag on the sea-surface, i.e., u*. The most frequently used u*-based formulation was developed by Jähne et al. [94] [Equation (16a)] from wave tank experiments. In this formulation the transfer velocities were attributed to surface renewal by micro-scale wave breaking. Hence, this formulation is adequate for comprehensive physically-based algorithm splitting the kwind and kbubble components of transfer velocity [37][38][39], being usually conjugated with the kbubble formulation by Woolf [38]. This specific ensemble yields RMSD = 9.33 cm•h −1 . Mackay and Yeun [92] performed wave tank experiments to measure the transfer velocities of 11 organic solutes with solubilities ranging from extremely low to extremely high. Their formulation [Equation (16b)], representing the average behavior of these organic molecules, yields extreme transfer velocities about 10 times higher than everything else ever reported. Zhao et al. [40] used data from the open ocean, coastal ocean and lakes, with wind speeds ranging from 1 m•s −1 to 34 m•s −1 . They calibrated a generalist formulation [Equation (16c)] with the effect of wave breaking implicit, and thus unsuited for comprehensive physicallybased algorithms. Under the stronger winds, this formulation estimated too low transfer velocities, resulting in RMSD = 12.5 cm•h −1 . Landwehr et al. [31] developed two formulations with similar structures but calibrated to different data sets [Equation (16d,e)]. Their negative intercepts led to kwind Figure 6. Estimation of k wind term. (a) The transfer velocity standardized for CO 2 in fresh-water at 20 • C (k 600 ) estimated from (EC data) direct observation obtained through EC method for 2014 and 2015, and relation with (u 10n ) wind speed (c p /u 10n ) and wave age. Waves aged c p /u 10n > 1 correspond to overdeveloped sea. The EC transfer velocities were compared with (a) those estimated from the u 10 -based formulations by Jacobs et al. [83] (Jac99), Wanninkhof [18] (Wan92) and Wanninkhof and McGillis [84] (WMG99), and (b-f) with comprehensive physically-based algorithms (u * -based) by Esters et al. [30]. [Equations (17) and (18)]. (a,b) k 600 estimated from Esters et al. [30] algorithm and its relation with u*. (d-f) validation of k 600 estimated from from Esters et al. [30] algorithm against k 600 estimated from direct observation obtained through EC method.
Turbulence-driven transfer velocities should be more accurately inferred from the wind effectively exerting drag on the sea-surface, i.e., u * . The most frequently used u * -based formulation was developed by Jähne et al. [94] [Equation (16a)] from wave tank experiments. In this formulation the transfer velocities were attributed to surface renewal by micro-scale wave breaking. Hence, this formulation is adequate for comprehensive physically-based algorithm splitting the k wind and k bubble components of transfer velocity [37][38][39], being usually conjugated with the k bubble formulation by Woolf [38]. This specific ensemble yields RMSD = 9.33 cm·h −1 . Mackay and Yeun [92] performed wave tank experiments to measure the transfer velocities of 11 organic solutes with solubilities ranging from extremely low to extremely high. Their formulation [Equation (16b)], representing the average behavior of these organic molecules, yields extreme transfer velocities about 10 times higher than everything else ever reported.
Zhao et al. [40] used data from the open ocean, coastal ocean and lakes, with wind speeds ranging from 1 m·s −1 to 34 m·s −1 . They calibrated a generalist formulation [Equation (16c)] with the effect of wave breaking implicit, and thus unsuited for comprehensive physically-based algorithms. Under the stronger winds, this formulation estimated too low transfer velocities, resulting in RMSD = 12.5 cm·h −1 . Landwehr et al. [31] developed two formulations with similar structures but calibrated to different data sets [Equation (16d,e)]. Their negative intercepts led to k wind estimates lower than those by the previous formulations, and to the estimation of k wind = 0 for u 10 = 2.5 m·s −1 and u 10 = 2 m·s −1 , respectively. This opposes the general acceptance of effective transfer velocities even under zero winds. Thus, despite RMSD of 8.7 cm·h −1 and 8.24 cm·h −1 , respectively, they are unsuited for application at the coastal ocean featuring frequent low winds. Other formulations use Equation (16f) to relate the transfer velocity to the water-side friction velocity (u *w ) [20,93,94]. In this case, the β s and n s parameters reflect the effects of capillary-gravity waves, surface tension and surfactants. Sc w is the Schmidt number of water, dependent on temperature and salinity. We could not test these formulations due to the lack of data on shorter waves. All these equations yield k wind in units of m·s −1 .
Evidence points to k wind scaling preferably with the turbulent kinetic energy dissipation rate ε [19,[22][23][24]27,28,30,31,39,53,82] than with u * . The major difficulty with this approach is to determine the vertical profile of ε and its integration depth. In their recent algorithm, Esters et al. [30] start from Equation (17a), where A is a scaling constant and υ w the kinematic viscosity of water. Following the Law of the Wall, ε at the surface (ε s ) was given by Equation (17b) [21,27,30,39,53]. However, this formulation only accounts for shear-induced turbulence, and not for turbulence from wave breaking. The integration depth z corresponds to the thickness of the viscous sub-layer, being inversely proportional to the wind speed [Equation (17c)]. Esters et al. [30] observed an off-set (here named δ) ≈ 0.394 between model forecasts and their observations. Equation (17a-c) were integrated into Equation (17d). However, its incorporation in FuGas 2.5 conjugated with Equation (1) required a constant with the value 660.
A can either be a proportionality constant or scale allometrically with ε as A = βε γ , β and γ being fixed coefficients. Esters et al. [30] propose specific A for CO 2 and DMS, and list the coefficients for their alternative calibrations. Beware the typo in Esters et al. [30] publication, with β erroneously having negative signs. The CO 2 coefficients are also provided in FuGas, whereas for other gases these must be updated from the literature. This formulation calibrated with n = 1/2, constant A, z 0 estimated from the wave field, and conjugated with the k bubble formulation by Woolf [38], showed a mild variability in the gas transfer velocities forecasted under similar wind speeds, since they now depended on ε instead ( Figure 6). However, they did not differ much from the previous formulations, and still largely misfitted the observed transfer velocities ( Figure 6). The best fit was provided when A = 0.25 and δ = 1, yielding RMSD = 7.62 cm·h −1 . Given the large variability in the observed transfer velocities, the uncertainty about their accuracy, and their systematic bias under swell, no solution can be considered overall better. The model implementation where A scaled with ε did not improve the fit, yielding RMSD = 10.28 cm·h −1 (Figure 6).
The k wind effect on transfer velocity is regulated by the Schimdt number, i.e., the ratio between the kinematic viscosity of water and the molecular diffusivity of the trace gas in water [see Equation (1)]. Its exponent n was demonstrated to scale with the wave field, being usually considered 2/3 for a smooth solid wall model to 1/2 for a free wavy surface [20,21,30,49,93,94]. Esters et al. [30] also provide the possibility of n depending on u *w instead of the wave field [Equation (18)]. However, when included in FuGas 2.5, this represents an implicit dependency of n from the wave field. n = −0.22 · log 10 (u * w ) + 0.13 (18) Applied to the Baltic Sea data, the variable n exponent was generally 0.58 < n < 0.71 with just four exceptions going beyond to a maximum of 0.86. This update did not improve the model fit to our data ( Figure 6), yielding a RMSD = 13.9 cm·h −1 . This was not surprising, since the Esters et al. [30] LOW formulation is only valid in situations in which no waves are present, so that the prevailing turbulence is purely shear induced. It was conjugated with the k bubble formulation by Woolf [38], which considers only the whitecap fraction relation with bubble-mediated gas transfer, and not with turbulence induced by breaking waves, as some other k bubble formulations do. Therefore, this ensemble neglects the gas transfer mediated by turbulence from wave breaking.
Terray et al. [53] estimated the depth profile of ε accounting for breaking waves and for a transition layer dependent on the wave size. From it, Esters et al. [30] developed an alternative formulation that we adapted [Equation (19a)]. In this case,ċ is the effective wave velocity, related with the energy flux from wind to waves, and given by Equation (19b) after reworking the original model. However, the original fit by Terray et al. [83], based exclusively on young waves, applied to old waves, led to negativeċ/c p and consequently negativeċ, both being illogical. This happened in most of the observations in our survey. We re-fit Terray's data to a Weibull curve (Figure 7a) and reworked to estimateċ [Equation (19c)]. Under low and moderate winds, the waves were generally old (u * /c p < 0.05) and the k wind was particularly sensitive to the x-axis displacement of the Weibull (= 18) and to its increment rate (= 2) determining the shape of this side of the function. Under higher winds the waves were younger (u * /c p > 0.05) and the k wind became also sensitive to the asymptotic y (= 0.5).
This formulation showed an improved ability to adapt the transfer velocity forecasts to the local conditions (Figure 7a,b). However, different calibrations adapted better to different groups of observations, and no specific calibration could adjust better to the whole set of observations. Consequently, the new transfer velocity fits to the transfer velocities, estimated from direct observation obtained through the EC method, were not better than the fits provided by the previous algorithms (Figure 7d, RMSD = 8.4 cm•h −1 ).

Conclusions
Transfer velocities obtained through the EC method are inherently subject to large uncertainty, making model calibration, validation and comparisons among formulations difficult. Nevertheless, our simulations showed that comprehensive physical-based algorithms are more sensitive, and better adaptable to local conditions, than simpler u10-based empirical formulations (Figure 8a). The more comprehensive algorithms in FuGas 2.5 included factors besides the wind speed that have been proved fundamental drivers of transfer velocities, such as sea-surface roughness and atmospheric stability. In the coastal ocean, and particularly under low winds, these factors lead to gas transfer velocity estimates systematically larger than the estimates by empirical u10-based formulations. Nevertheless, many other factors are still to be added, and those included can still be improved. One such important update is the effect of shorter (capillary) waves, when these act as roughness elements on the sea-surface [130][131][132], leading the air-flow to separate. Other potentially important updates are the saturation of the turbulent transfer velocities [49,[57][58][59][60][61] and its dependency on the relative air and water motions [49]. Empirical u10-based formulations can see their predictive ability improved if Figure 7. Estimation of k wind term adapting the algorithm by Esters et al. [30] based on the formulation by Terray et al. [53] for the ε depth profile [Equation (19)]. (a) Terray et al. [53] function and FuGas re-fit; (b) and (c) Esters et al. [30] calibration; (d) Esters et al. [30] validation against EC-derived transfer velocities. (T96 data) Terray et al. [53] data, (T96 fit) Terray et al. [53] fit, and (FuGas fit) FuGas re-fit of effective wave velocity. (EC data) transfer velocities estimated from the Eddy-Covariance fluxes.
This formulation showed an improved ability to adapt the transfer velocity forecasts to the local conditions (Figure 7a,b). However, different calibrations adapted better to different groups of observations, and no specific calibration could adjust better to the whole set of observations. Consequently, the new transfer velocity fits to the transfer velocities estimated from direct observation obtained through the EC method were not better than the fits provided by the previous algorithms ( Figure 7d, RMSD = 8.4 cm·h −1 ).

Conclusions
Transfer velocities obtained through the EC method are inherently subject to large uncertainty, making model calibration, validation and comparisons among formulations difficult. Nevertheless, our simulations showed that comprehensive physical-based algorithms are more sensitive, and better adaptable to local conditions, than simpler u 10 -based empirical formulations (Figure 8a). The more comprehensive algorithms in FuGas 2.5 included factors besides the wind speed that have been proved fundamental drivers of transfer velocities, such as sea-surface roughness and atmospheric stability. In the coastal ocean, and particularly under low winds, these factors lead to gas transfer velocity estimates systematically larger than the estimates by empirical u 10 -based formulations. Nevertheless, many other factors are still to be added, and those included can still be improved. One such important update is the effect of shorter (capillary) waves, when these act as roughness elements on the sea-surface [130][131][132], leading the air-flow to separate. Other potentially important updates are the saturation of the turbulent transfer velocities [49,[57][58][59][60][61] and its dependency on the relative air and water motions [49]. Empirical u 10 -based formulations can see their predictive ability improved if using the wind speed under atmospherically neutral conditions (u 10n ) (Figure 8b). In fact, after wind speed, atmospheric stability is probably the most important factor governing gas transfer velocity across the atmosphere-ocean interface.
using the wind speed under atmospherically neutral conditions (u10n) (Figure 8b). In fact, after wind speed, atmospheric stability is probably the most important factor governing gas transfer velocity across the atmosphere-ocean interface. Figure 8. Comparing transfer velocities predicted by comprehensive physical-based algorithms and empirical u10-based algorithms. Comprehensive algorithm included kwind term by Terray et al. [53] and kbubble term by Woolf [38]. Empirical u10-based algorithms by Jacobs et al. (1999) [83], Wanninkhof (1992) [18] and Wanninkhof and McGillis (1999) [84]. Please confirm if the (a) and (b) should be added in this figure caption.

Code and Data Availability
The

Conflicts of Interest:
The authors declare no conflict of interest. Figure 8. Comparing transfer velocities predicted by comprehensive physical-based algorithms and empirical u 10 -based algorithms. Comprehensive algorithm included k wind term by Terray et al. [53] and k bubble term by Woolf [38]. Empirical u 10 -based algorithms by Jacobs et al. (1999) [83], Wanninkhof (1992) [18] and Wanninkhof and McGillis (1999) [84]. Please confirm if the (a) and (b) should be added in this figure caption.

Code and Data Availability
The FuGas was first published by Vasco Vieira in 2012. The latest FuGas 2.5, data and videos related with previous publications are available at http://www.maretec.org/en/models/fugas.