The FuGas 2.3 Framework for Atmosphere–Ocean Coupling: Comparing Algorithms for the Estimation of Solubilities and Gas Fluxes

Accurate estimates of the atmosphere–ocean fluxes of greenhouse gases and dimethyl sulphide (DMS) have great importance in climate change models. A significant part of these fluxes occur at the coastal ocean which, although much smaller than the open ocean, have more heterogeneous conditions. Hence, Earth System Modelling (ESM) requires representing the oceans at finer resolutions which, in turn, requires better descriptions of the chemical, physical and biological processes. The standard formulations for the solubilities and gas transfer velocities across air–water surfaces are 36 and 24 years old, and new alternatives have emerged. We have developed a framework combining the related geophysical processes and choosing from alternative formulations with different degrees of complexity. The framework was tested with fine resolution data from the European coastal ocean. Although the benchmark and alternative solubility formulations generally agreed well, their minor divergences yielded differences of up to 5.8% for CH4 dissolved at the ocean surface. The transfer velocities differ strongly (often more than 100%), a consequence of the benchmark empirical wind-based formulation disregarding significant factors that were included in the alternatives. We conclude that ESM requires more comprehensive simulations of atmosphere–ocean interactions, and that further calibration and validation is needed for the formulations to be able to reproduce it. We propose this framework as a basis to update with formulations for processes specific to the air–water boundary, such as the presence of surfactants, rain, the hydration reaction or biological activity.


Introduction
The dynamics of the atmosphere-ocean gas exchanges plays a central role in the Earth's climate, with the ocean acting simultaneously as a sink and source of greenhouse gases and a source of dimethyl sulfide (DMS) to the atmosphere.However, and despite the observed seasonal, inter-annual and regional variability [1][2][3], the open ocean is generally assumed to act as a sink for atmospheric CO 2 .In the sub-polar regions, for instance, the solubility pump retrieves large amounts of greenhouse gases from the atmosphere, which are then advected to the deep ocean [1,2].On the other hand, at the surface of coastal waters the balances and fluxes of CO 2 , CH 4 , N 2 O and DMS are very heterogenic, mostly due to factors like upwelling, plankton productivity, and land-originated loads of both natural and anthropogenic cause.Consequently, besides DMS [3,4] the coastal ocean can be, at least seasonally, also a source of CO 2 [1,2,[5][6][7][8][9][10][11], CH 4 [2,10,12], and N 2 O [2,10,[13][14][15] to the atmosphere.Thus, the uncertainty on the atmosphere-coastal ocean exchange of CO 2 begs the question about the overall budget of the global ocean [1,3,6,16].
The atmosphere-ocean flux of gases is controlled by the unbalance between the air and water concentrations (the latter mediated by solubility), and by the transfer velocity [17][18][19][20][21][22].The solubility of gases at the surface can be estimated by algorithms with different physical and chemical backgrounds [20][21][22][23][24][25][26][27].The transfer velocity depends on numerous processes occurring at the sea surface [17][18][19][20][21][22]28].However, large-scale simulations of relevant processes for the air-water gas exchange at the coastal ocean and their inherent variability poses serious difficulties for Earth System Models (ESM) and marine biogeochemical models.Constraints related to computational demands explain why such models usually simulate the biosphere at decadal or centennial time scales with daily intervals and spatial resolutions of hundreds to 1000 km.Such are the cases of ESM applications by the Intergovernmental Panel on Climate Change (IPCC), Max Planck Institute (MPI) or Centro Euro-Mediterraneo sui Cambiamenti Climatici (CMCC).Hand in hand with the low resolution for space and time, the atmosphere-ocean gas fluxes are estimated with simple algorithms that disregard the complexity of processes recently unveiled at the coastal ocean.Hence, the standard in current ESM is the parameterization by Wanninkhof [29], relying on wind speed at 10 m elevation under atmospherically neutral conditions (u 10N ) as the sole driver of transfer velocity.ESM relies on the dominance of strong winds over the global ocean to support the coarse resolution of geophysical process.Contradicting this rationale, Wanninkhof et al. [19] demonstrated that the largest amount of atmosphere-ocean CO 2 exchange occurs under low-to-moderate winds.Concomitantly, there emerged empirical formulae best fitting low wind data collected from estuaries.Such are the formulations proposed by Carini et al. [30] or Raymond and Cole [31], accounting only for u 10 , or by Borges et al. [32], adding current drag with the bottom.Meanwhile, new evidence emerged for the effects of the sea surface cool-skin and warm-layer [33][34][35][36][37][38], atmospheric stability [39,40], wave-breaking [41][42][43][44][45][46][47][48][49][50], water surface renewal [28,[51][52][53][54][55][56][57], surfactants [28,[58][59][60] and rain [61][62][63][64], just to mention a few examples from an extensive list of published works.In addition to their single influence, these factors also interact.For instance, while surfactants increase surface tension, although decreasing the water surface renewal and micro-scale wave breaking [28,65], there is no report of them affecting whitecap from breaking waves.Surfactants are expected to be more relevant on the coastal ocean due to their supply from land and their effect being more pronounced under low-moderate winds [60,66].
There are many simplistic formulae estimating the gas transfer velocity from the few factors that were determinant for the specific set of data used in their calibration.However, modelling the coastal oceans with fine resolution requires algorithms providing adequate representation of the many processes involved.Developing such an algorithm demands a modular framework able to be updated with the best formulation for each of the mediator processes, and to test and compare different algorithms as well as different degrees of complexity.Hence, we developed the FuGas (Flux of Gases) framework [21,22], currently in its version 2.3.This simulation tool includes 67 alternative formulations to account for factors such as solubility, wind or current-mediated turbulence, atmospheric stability, sea-surface roughness, breaking waves, air and water viscosities, temperature and salinity.This framework is also able to automatically select simpler algorithms in particular locations with insufficient data for more comprehensive solutions.The framework can be used to estimate the air-water fluxes of any atmospheric gas by selecting the respective constants.In this work we rely on the FuGas 2.3 to compare estimates of CO 2 transfer velocities obtained by simpler empirical formulations with some more comprehensive alternatives focusing on the sea state.We tested the effect of wave-state with field data from the Baltic Sea.Also, we coupled the Weather Research and Forecasting (WRF) atmospheric model to the WaveWatch III (WW3)-NEMO oceanographic model using simulated data from the European coastal ocean.These new tests are particularly relevant given recent findings that the sea state in fetch-limited locations (as is the coastal ocean) has an influence on the transfer velocity, by contrast with that from the fetch-unlimited open ocean [50,55,57].

Methods
Air-water gas fluxes result from the interaction of two factors: (i) the imbalance between the gas concentrations in the air and in the water sets the strength and direction of the flux; and (ii) the resistance the medium does for being crossed by the flow [17][18][19][20][21][22].Most formulae estimate the flux from F = k w •k H cp•∆p gas .The k H cp is the Henry's constant for gas solubility in its C w /p a form, where p a is the gas partial pressure in the air and C w its concentration in the water.The ∆p gas is the difference between the partial pressures of the gas in the air (p a ) and in the water (p w ) i.e., ∆p gas = p w − p a .In this form, a positive flux represents uptake by the atmosphere while a negative flux represents uptake by the ocean.The p a must be corrected for the partial pressure of water vapour considering saturation over the sea surface: p moist = (1 − p H2O /P) p dry .This conversion is detailed in Section 2.1 below.The k w is the transfer velocity of gases across the sub-millimetre thick water surface layer.For moderately soluble gases the airside resistance is not negligible.In these cases, the double layer model [17][18][19][20][21][22] estimates the flux taking into consideration both the water-side and air-side surface layers and thus, . The C a and C w are the concentrations of the gas in air and water and the k H is Henry's constant in its equivalent dimensionless quantity (C a /C w ).The transfer velocity is estimated from both layers as S1 provides a list of the variables and constants presented in this article.Other variables and constants used in subsidiary calculations are not included.The link to software and data is provided in the Supplementary Materials.

Solubility
Sarmiento and Gruber [2] compiled the algorithm for the k H cp dependence on temperature and salinity provided by Weiss [23] and Weiss and Price [24].We converted it to its corresponding dimensionless k H preserving the constants required to estimate Bunsen's solubility coefficient (β) (Equation (1) or Equation ( 2)).Solubility coefficients were estimated from the Virial expansion (Equation (3)), where B was β or β/V m , depending on which gas it was applied to (Table 3.2.2 in Sarmiento and Gruber [2]).Our software automatically detected the gas from the a i coefficient.When B = β the k H was estimated from Equation (1), whereas when B = β/V m the k H was estimated from Equation (2).This formulation accounted for fugacity (f) of non-ideal gases (Equation ( 4)) and corrected the gas partial pressure for moisture effects from p moist = (1 − p H2O /P) p dry considering water vapour saturation over the sea-surface (Equation ( 5)).Calculations required the air pressure (P), water temperature (T w ), salinity (S), ideal gas law constant (R), molar volume of the specific gas (V m ) and the molar volume of ideal gases (V ideal ).
Johnson [20] developed an algorithm from another approach that corrected for the effects of temperature and salinity on the molecular and thermodynamic properties of the water, its solutes and the specified gas, based on a ideal gas behaviour.His formulae were developed from the compilation by Sander [27] (available on the web since 1999) of the k H cp for nearly all gases in the atmosphere at 25 • C (298.15 K) and 0 ppt.Then, Equation ( 6) converted the k H cp to k H at a given temperature and 0 ppt salinity.The term −∆ soln H/R reflected the temperature dependence of solubility, having a value of 2400 for CO 2 , 1700 for CH 4 and 2600 for N 2 O.The correction to a given salinity (Equation ( 7)) relied on the empirical Setschenow constants (K S = θ•logV b ) reporting the effect of electrolytes salting-out gases proportionally to their liquid molar volume at boiling point (V b ).The V b was estimated using the additive Schroeder method, whereas θ was estimated from Equation (8) using a provisional k H# = 0.0409/k H cp. k H,0 = 12.1866

Transfer Velocity
The simpler transfer velocity formulations rely on the wind velocity 10 m above the sea surface (u 10 ) under atmospherically neutral conditions.Among then, the formulation by Wanninkhof [29] (henceforth also referred to as 'Wan92') became the standard used in ESM and satellite data processing (Equation (9a,b)).It includes the Schmidt number of the water (Sc w ), related to viscosity, raised to an exponent sensitive to the assumed model of surface renewal.For sparingly soluble gases or for CO 2 in alkaline waters, and under low winds, the transfer velocity is enhanced by the hydration reaction (b T ), which is sensitive to temperature.Other simple empirical formulations based only on u 10 [30,31], or also accounting for bottom-generated turbulence [32], used data collected in estuaries under low wind conditions.However, modelling the coastal ocean at finer resolutions requires an enhanced representation of the multitude of processes involved.Hence, we updated the framework by Vieira et al. [21], with the k w being decomposed into its shear produced turbulence (k wind + k current ) and bubbles from whitecapping (k bublle ) forcings [32,41,42,44].Sc w was determined from temperature and salinity following Johnson [20]: The formulation by Zhao et al. [43], merged k wind within the k bubble (Equation (11a)), being this term dependent on the wave breaking parameter (R B given by Equation (11b)).The u * was the friction velocity i.e., the velocity of wind dragging on the sea surface, and f p the peak angular frequency of the wind-waves (in rad•s −1 ).We estimated the kinematic viscosity of air (υ a ) from Johnson [20].
Later, Zhao and Xie [47] developed a more elaborate algorithm.While the output was still given by Equation (11a), R b was now estimated from Equation (12a).R HU was a turbulent Reynolds number (Equation (12b)) accounting for the sea state by using the significant wave height (H s ).The wave age (β) was estimated from Equation (12c) using the gravitational acceleration constant (g).While Zhao and Xie [47] tested their algorithm with the drag coefficient (C D ) estimated from several empirical formulations, here we estimated C D from its definition (C D = u * 2 /u 10 2 ).This was only possible because u * was estimated from u 10 using the wind log-linear profile as presented below.
The formulation by Woolf [44] split the estimation of k wind and k bubble .The k wind , taken from Jähne et al. [66], related to the transfer mediated by the turbulence generated by wind drag (Equation ( 13)).The k bubble = 850 W related to the transfer mediated by the bubbles generated by breaking waves.The whitecap cover was estimated from W = 4.02 The friction velocity (u * ) was estimated from the wind log-linear profile (WLLP: Equation ( 14)) accounting for wind speed at height z (u z ), atmospheric stability of the surface boundary layer (through ψ m ) and sea-surface roughness (through the roughness length z 0 ).Here, κ is von Kármàn's constant.Historically, the WLLP originated from Monin-Obukhov similarity theory [67].
Roughness length (z 0 ) is the theoretical minimal height (most often sub-millimetric) at which wind speed averages zero.It is dependent on surface roughness and often used as its index.It is more difficult to determine over water than over land as there is a strong bidirectional interaction between wind and sea surface roughness.Taylor and Yelland [68] proposed a dimensionless z 0 dependency from the wave field, increasing with the wave slope (Equation ( 15)) as estimated from the significant wave height (H s ) and peak wave period (L p ). Due to the bidirectional nature of the z 0 and u * relation, we also tested an iterative solution (iWLP) where Equation ( 15) was used as a first guess for the z 0 and Equation ( 14) for its subsequent u * .A second iteration re-estimated z 0 from the COARE 3.0 [69] adaptation of the Taylor and Yelland [68] formulation, whereby a term for smooth flow was added (Equation ( 16)), and u * was estimated back again from Equation (14).Applying four iterations were enough for an excellent convergence of the full data array.The coefficients proposed by Taylor and Yelland [68] yielded implausibly high z 0 when applied to steep young waves originated by a storm over the central Mediterranean Sea and by strong Mistral blowing off-shore of southern France, corresponding to 2.85% of our observations.These extreme z 0 led to absurdly high u * and k w , and thus we imposed a maximum roughness length of 0.01 m.This value roughly matches the maximum expected from the application of alternative formulations for the estimation of z 0 from the wave field [70][71][72][73][74].
Atmospheric stability characterized the tendency of the surface boundary layer (SBL) to be well mixed (unstable SBL with ψ m < 0) or stratified (stable SBL with ψ m > 0).The ψ m was inferred from the bulk Richardson number (Ri b ), weighting the air vertical heat gradient and kinetic energy (Equation ( 17)).Hence, it required the air virtual potential temperature (T v ) and wind speed at two heights (z) within the SBL.Alternatively, the Ri b was directly estimated from the air potential temperature (T p ) i.e., replacing T v and, thus, neglecting humidity [75].
T v was estimated from T p and humidity, following Equation ( 18) [39] or Equation ( 19) [67].Each equation requires different humidity inputs, namely the observed specific humidity (q), relative humidity (h r ), liquid water mixing ratio (r L ) and water vapour mixing ratio at saturation (r sat ): T p and r sat were estimated from the air temperature (T ka ) and pressure (P) [67]: ln(e sat ) = ln 0.61078 At surface height z ≈ 0 m, the wind velocity must be set to the theoretical u 0 = 0, temperature given by the SST [34,39,69] with rectification for cool-skin and warm-layer effects [34][35][36][37][38], and humidity set to the saturation level at P 0 and T 0 using q = r sat /(1 + r sat ) [69].The Ri b was used to estimate the length L from the Monin-Obukhov's similarity theory, a discontinuous exponential function tending to ±∞ when Ri b tends to ±0 and tending to ±0 when Ri b tends to ±∞.Ri b and L were used to estimate ψ m following Stull [67] or Lee [75] algorithms.The transfer velocity of moderately soluble gases, besides taking into consideration the molecular crossing of the water-side surface layer, should also take into consideration the molecular crossing of the air-side surface layer [17][18][19][20][21][22].We compared values from a single layer model and two-layer "thin film" model [17], the later estimating the air-side transfer velocity (k a ) from the COARE formulation as in Equation ( 23) [46].C D is the drag coefficient and Sc a the Schmidt number of air, which were determined for a given temperature and salinity following Johnson [20].

Validation with Field Data
The algorithms presented in Section 2.2 were applied to observed oceanographic and atmospheric variables (Level 2 data).These are derived geophysical variables at the same resolution and location as the instrument source data (Level 1).For these estimates we used the FuGas 2.2 single processing ensemble.The forecasted CO 2 transfer velocities were compared to those inferred from k w = F/∆p gas , where F were the CO 2 vertical fluxes measured by eddy-covariance, and the ∆p gas were the differences between measured atmospheric and oceanic CO 2 concentrations.The field sampling occurred from 22 May 2014 to 26 May 2014 using the atmospheric tower at Östergarnsholm in the Baltic Sea (57 • 27 N, 18 • 59 E), the Submersible Autonomous Moored Instrument (SAMI-CO 2 ) 1 km away and the Directional Waverider (DWR) 3.5 km away, both south-eastward of the tower [9,76].The atmospheric measurements were done at 12 m heights.The DWR measured temperatures at 0.5 m depth, taken as representative for the sea-surface.Salinity was obtained from the Asko mooring data provided by the Baltic In-Situ Near-Real-Time Observations available in Copernicus Marine catalogue.
The eddy-covariance (E-C) measurements were subject to quality control.We only used the fluxes for which the wind direction set the SAMI-CO 2 and DWR in the footprint of the atmospheric tower (100 • < wind direction < 160 • ).We rejected periods for which there were abrupt changes in wind speed or direction as these represented strong violations of the required homogeneous conditions [77].An extra filter was applied to remove measurements in which wave age denoted large dominance of swell conditions (β > 2), since they have been identified as biasing measurements of transfer velocity [47,48].When the ratio is >1 the waves move faster than the wind, indicating the presence of swell [47,48].Using only the observations passing the quality control check, the air-water CO 2 fluxes measured by eddy-covariance were smoothed over 30 min bins and corrected according to the Webb-Pearman-Leuning (WPL) method [78].

Application within Geoscientific Models
The transfer velocities of CO 2 , CH 4 and N 2 O were estimated from atmospheric and oceanographic re-analysed (Level 4) data retrieved for the European shores (Mediterranean and North Atlantic to the tip of Scotland) from the 24 May 2014 at 06 h to the 27 May 2014 at 00 h.All variables were interpolated to the same 0.09 • grid (roughly 11 km at Europe's latitudes) and 1 h time steps.This resulted in a data set with 17 variables × 41,776 locations × 66 time instances.For a faster processing of the ≈2 Gb of Level 4 data, we used the multiple processing ensemble of the FuGas 2.3.
The atmospheric variables were retrieved from the standard operational application of the WRF by Meteodata.cz,with model output at 9 km horizontal and 1 h time resolutions.The vertical thickness of the WRF layers varied with space and time.The atmospheric variables where retrieved from the two lowest levels within the SBL.These levels refer to the horizontal layers in the atmospheric model, and are not to be mistaken by the data processing level.Over the ocean, these levels occurred roughly at 0 m and 12 m heights.The second level provided the wind velocity, temperature, pressure and humidity at z meters above sea surface.At surface height z ≈ 0 m, the wind velocity was set to the theoretical u 0 = 0, and temperature was given by the SST without rectification for cool-skin and warm-layer effects due to the lack of some required variables.Yet, these effects tend to compensate each other [34,37,38].Air pressure was given by the WRF at the lower level (≈0 m) and humidity was set to the saturation level using P 0 and T 0 .
The sea surface temperature (SST, here also named T w ) and salinity (S) were estimated by the NEMO modelling system provided in the MyOcean catalogue with 1/12 • and 1 day resolutions.The WW3 wave field data for the Mediterranean Sea was supplied by the Istituto Nazionale di Geofisica e Vulcanologia (INGV) using the WW3-NEMO modelling system at 0.0625 • and 1 h resolutions [79], and for the North Atlantic by Windguru at roughly 0.5 • and 3 h resolutions.The variables included significant wave height (H s ) and peak frequency (f p ) for wind sea i.e., disregarding swell.A few aspects did not correspond to the ideal data format for atmosphere-ocean coupling and required further calculations: (i) The peak wave length (L p ) was estimated from the peak frequency assuming the deep-water approximation: L p = 2πg/f p 2 ; (ii) the Windguru data did not provide a wind sea component where (and when) the wind was too low.For these missing cases, the lowest H s were attributed and L p was simulated everywhere else; and (iii) the Windguru and the INGV data overlapped along the Iberian shores, in which case the INGV was given a 2:1 weight over the Windguru data.SST and S were used to estimate how the solubility algorithms could lead to different estimates of greenhouse gases dissolved below the ocean surface.First, the solubility was converted to the C w /C a form corresponding to 1/k H .Then, the differences in solubility between the two algorithms, as estimated from Equation (24), led to the differences in the estimates of greenhouse gases dissolved per cubic meter immediately below the ocean surface (Equation ( 25)).
For these calculations, time variability was disregarded assuming the average ∆s over the 66 h time interval.The ∆C w was converted to mass units using the molar weights of the specific gas (M gas ), respectively 44.01 for CO 2 , 16.043 for CH 4 and 44.013 for N 2 O.The C a was estimated from the atmospheric pressure (P in atm) and the partial pressure (ρ gas ) of CO 2 , CH 4 or N 2 O, 390 ppm, 1.75 ppm and 0.325 ppm respectively [80], assuming that they were approximately uniform all over the atmospheric surface boundary layer (SBL).The ρ gas was divided by 10 6 to re-scale from ppm to dimensionless.Using the ideal gas law, we divided by R and T a (assumed T a = T w at the sea-surface), and multiplied by 101325 to convert atm back to Pascal.

Solubility Estimates
The solubility formulae were compared for the range of environmental conditions commonly found in nature: water temperature (T w ) ranged from −2 • C to 30 • C at 1 • C intervals and a salinity (S) ranged from 0 ppt to 38 ppt at 1 ppt intervals.In general, both formulations estimated similar solubilities despite their distinct chemistry backgrounds, as demonstrated by their ratio k H [20] /k H [2] close to 1 (Figure 1a-c).Nevertheless, differences from 10% to −5% for CO 2 , 11% to 5% for CH 4 , and 7% to −5% for N 2 O were observed in freshwaters.Determining which formulation provides more accurate estimates of greenhouse gas (GHG) solubilities in freshwater is fundamental for Earth System Modelling, because rivers and freshwater reservoirs are important sources of GHG to the atmosphere, even if only seasonably [7][8][9][10][11][12][13][14]16].In marine waters, our simulations found differences of up to 12% for CO 2 , 9% for CH 4 and 7% for N 2 O, suggesting that for marine waters is also important to determine which formulation provides more accurate estimates of GHG solubilities.Afterwards, both formulae were applied to the conditions observed at the European coastal ocean during the experiment (the Level 4 SST, S and P).The surface water temperature changed significantly and there were large fresh water inflows from the Black Sea and the Baltic Sea (Figure 2 and Supplementary Video S1).Solubility estimates were averaged over a 66 h time interval and compared at each location by their ratio (Figure 1d-f).Higher divergences reached +4.5% for CO 2 solubility in cooler marine waters, +5.65% for CH 4 solubility in warmer saltier waters, and −2.1% for N 2 O solubility in brackish waters.Hence, the algorithms diverge by as much as 0.045 mol•mol −1 of CO 2 , 0.0015 mol•mol −1 of CH 4 and 0.012 mol•mol −1 of N 2 O (i.e., mol of gas in the ocean surface per mol of gas in the atmosphere) in some of the most sensitive situations for Earth System Modelling and satellite data processing, namely: (i) the cooler marine sub-polar waters where the solubility pump traps greenhouse gases and carries them to the deep ocean [2], and (ii) the warmer waters at the coastal ocean and seas, which have regularly been observed having greenhouse gases and DMS dissolved in concentrations highly unbalanced with those of the atmosphere [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16].
The discrepancies between solubilities estimated by either algorithm led to differences in the estimates of GHG dissolved per cubic meter immediately below the ocean surface (Figure 3).The major differences occurred in the regions with the cooler waters, the warmer waters and the brackish waters.These biases in the estimated total amount of greenhouse gases for the European coastal ocean's surface during late May 2014, may be an indicator of higher global biases.Hence, similar comparisons with data from other locations and/or periods are required.Performing these comparisons is facilitated the by the FuGas framework.The discrepancies between solubilities estimated by either algorithm led to differences in the estimates of GHG dissolved per cubic meter immediately below the ocean surface (Figure 3).The major differences occurred in the regions with the cooler waters, the warmer waters and the brackish waters.These biases in the estimated total amount of greenhouse gases for the European coastal ocean's surface during late May 2014, may be an indicator of higher global biases.Hence, similar comparisons with data from other locations and/or periods are required.Performing these comparisons is facilitated the by the FuGas framework.The discrepancies between solubilities estimated by either algorithm led to differences in the estimates of GHG dissolved per cubic meter immediately below the ocean surface (Figure 3).The major differences occurred in the regions with the cooler waters, the warmer waters and the brackish waters.These biases in the estimated total amount of greenhouse gases for the European coastal ocean's surface during late May 2014, may be an indicator of higher global biases.Hence, similar comparisons with data from other locations and/or periods are required.Performing these comparisons is facilitated the by the FuGas framework.
The SBL was stable i.e., 0.016 < Rib < 0.093.The sea surface was slightly to moderately rough (z0 < 0.49 mm).Under these conditions, the kw estimated from the E-C measurements were approximate to the kw estimated by both empirical u10-based parameterizations and extended physically based algorithms (Figures 4 and 5).The small kw fluctuations estimated by empirical u10-based parameterizations were the sole consequence of changes in water viscosity (as estimated by the Scw) driven by changes in water temperature.On the other hand, the physically based algorithms, by accounting for the effects of sea state, breaking waves and atmospheric stability, showed a greater adaptability to local conditions (Figure 5).As an example, depending on the sea state, the formulations by Woolf [44] or Zhao and Xie [47] estimated kw similar to that estimated by the Carini

Transfer Velocity Estimates from Field Data
The transfer velocity (k w ) scaled with u 10 (r = 0.91) and with the inverse of the wave age (r = −0.88)i.e., the younger the waves, the larger the k w (Figure 4).This fits the general theory that under low-moderate wind speeds, the k w is set by surface renewal and micro-scale wave breaking (i.e., by k wind ), increasing with steeper younger waves [28,51,65], and only under higher wind speeds is k w set by bubbles from breaking waves (i.e., k bubble ) [41][42][43][44][45][46][47][48][49][50].The preponderance of surface renewal and micro-scale wave breaking at the coastal ocean turn it more susceptible to additional factors driving k w .As an example, under the low-moderate wind speeds often verified at the coastal ocean, the hydration reaction enhances k w [19] and surfactants attenuate surface renewal [28,58,59,65,66].[20] to compilation by Sarmiento and Gruber [2].Colourscale: Δg m −3 .
The SBL was stable i.e., 0.016 < Rib < 0.093.The sea surface was slightly to moderately rough (z0 < 0.49 mm).Under these conditions, the kw estimated from the E-C measurements were approximate to the kw estimated by both empirical u10-based parameterizations and extended physically based algorithms (Figures 4 and 5).The small kw fluctuations estimated by empirical u10-based parameterizations were the sole consequence of changes in water viscosity (as estimated by the Scw) driven by changes in water temperature.On the other hand, the physically based algorithms, by accounting for the effects of sea state, breaking waves and atmospheric stability, showed a greater adaptability to local conditions (Figure 5).As an example, depending on the sea state, the formulations by Woolf [44] or Zhao and Xie [47] estimated kw similar to that estimated by the Carini  (E-C).The Jac99 (Jacobs et al. [81]) and WMG99 (Wanninkhof and McGillis [82]) are the empirical u 10 -based parameterizations estimating higher and lower k w within this wind range.
The SBL was stable i.e., 0.016 < Ri b < 0.093.The sea surface was slightly to moderately rough (z 0 < 0.49 mm).Under these conditions, the k w estimated from the E-C measurements were approximate to the k w estimated by both empirical u 10 -based parameterizations and extended physically based algorithms (Figures 4 and 5).The small k w fluctuations estimated by empirical u 10 -based parameterizations were the sole consequence of changes in water viscosity (as estimated by the Sc w ) driven by changes in water temperature.On the other hand, the physically based algorithms, by accounting for the effects of sea state, breaking waves and atmospheric stability, showed a greater adaptability to local conditions (Figure 5).As an example, depending on the sea state, the formulations by Woolf [44] or Zhao and Xie [47] estimated k w similar to that estimated by the Carini et al. [30] parameterization for the Parker river estuary, or similar to that estimated by the Wanninkhof [29] parameterization for the global ocean.Consequently, the root mean square deviation (RMSD) between the Woolf [44] formulation (henceforth 'W05va') and the E-C transfer velocities was 6.34 cm•h −1 , whereas between the Wan92 and the E-C transfer velocities it was 8.4 cm•h −1 .On relative terms, the physically based algorithm showed 75% the deviation of the empirical u 10 -based algorithm.Nevertheless, the physically based algorithms were still unable to accurately fit E-C observations (Figures 4 and 5).
Atmosphere 2018, 9, x FOR PEER REVIEW 11 of 18 et al. [30] parameterization for the Parker river estuary, or similar to that estimated by the Wanninkhof [29] parameterization for the global ocean.Consequently, the root mean square deviation (RMSD) between the Woolf [44] formulation (henceforth 'W05va') and the E-C transfer velocities was 6.34 cm•h −1 , whereas between the Wan92 and the E-C transfer velocities it was 8.4 cm•h −1 .On relative terms, the physically based algorithm showed 75% the deviation of the empirical u10-based algorithm.Nevertheless, the physically based algorithms were still unable to accurately fit E-C observations (Figures 4 and 5).With the optional coupling of the sea state and wind velocity forcings (through the iWLP joint estimate of z0 and u*), the FuGas enables an enhanced representation of local conditions that was demonstrated to be fundamental for the estimation of transfer velocities in fetch-limited situations [50,55,57].This coupling agrees with observations that kw scales with the turbulent kinetic energy dissipation at the sea surface (ε), and that ε is better reflected by the sea state [54,56].Accordingly, the COARE 3.0 included the wave state in the estimation of the roughness parameters essential for the transfer of mass, heat and momentum [69].We highlight the fact that the roughness Reynolds number is given by R = u*•z0/υa i.e., relying precisely on the u* and z0 that are jointly estimated by the FuGas from the SBL and wave field.
The estimation of kw as wind asymptotically tends to zero is generally considered uninteresting for open ocean formulations, although the equatorial global ocean consistently presents such low winds.The bulk of the empirical u10-based parameterizations undervalues the accuracy of the kw estimated under the lowest winds and make it asymptotically tending to zero as well.The parameterizations by Mackay and Yeun [87], McGillis et al. [84], Wanninkhof et al. [19] and Johnson [20] propose eliminating this misrepresentation of reality by making kw asymptotically tend to a constant.Arguably, the accurate estimation of kw under low wind speeds becomes more relevant when studying the coastal ocean and inland waters, and extended physically-based algorithms are better suited for the task.The COARE algorithm added a gustiness term stabilizing the kw in finite velocities under lighter winds [69].Clayson et al. [51], on the other hand, replaced the gustiness term by a capillary wave parameterization.The FuGas achieved this through the iterative estimation of With the optional coupling of the sea state and wind velocity forcings (through the iWLP joint estimate of z 0 and u * ), the FuGas enables an enhanced representation of local conditions that was demonstrated to be fundamental for the estimation of transfer velocities in fetch-limited situations [50,55,57].This coupling agrees with observations that k w scales with the turbulent kinetic energy dissipation at the sea surface (ε), and that ε is better reflected by the sea state [54,56].Accordingly, the COARE 3.0 included the wave state in the estimation of the roughness parameters essential for the transfer of mass, heat and momentum [69].We highlight the fact that the roughness Reynolds number is given by R = u * •z 0 /υ a i.e., relying precisely on the u * and z 0 that are jointly estimated by the FuGas from the SBL and wave field.
The estimation of k w as wind asymptotically tends to zero is generally considered uninteresting for open ocean formulations, although the equatorial global ocean consistently presents such low winds.The bulk of the empirical u 10 -based parameterizations undervalues the accuracy of the k w estimated under the lowest winds and make it asymptotically tending to zero as well.The parameterizations by Mackay and Yeun [87], McGillis et al. [84], Wanninkhof et al. [19] and Johnson [20] propose eliminating this misrepresentation of reality by making k w asymptotically tend to a constant.Arguably, the accurate estimation of k w under low wind speeds becomes more relevant when studying the coastal ocean and inland waters, and extended physically-based algorithms are better suited for the task.The COARE algorithm added a gustiness term stabilizing the k w in finite velocities under lighter winds [69].Clayson et al. [51], on the other hand, replaced the gustiness term by a capillary wave parameterization.The FuGas achieved this through the iterative estimation of the wind log-linear profile (the iWLP: Equations ( 15) and ( 16)).Under the lowest winds but as long as there were waves, the smooth flow term in the iWLP imposed finite u * , and consequently, finite transfer velocities (Figure 5).

Transfer Velocity Estimates from Level 4 Data
We performed simulations for the European costal ocean to compare between the ESM standard (the 'Wan92') and our extended physically-based alternatives.We show the comparison with the iWLP implemented with the k wind formulation by Jähne et al. [66] and the k bubble formulation by Woolf [44] using the kinematic viscosity of air (the W05va).Their k w estimates diverged mostly under unstable SBL, very rough or very smooth sea-surfaces, or higher friction velocities (Figure 6).The details of the simulations and the differences between k w estimates are presented hereafter.
Atmosphere 2018, 9, x FOR PEER REVIEW 12 of 18 the wind log-linear profile (the iWLP: Equations ( 15) and ( 16)).Under the lowest but as long as there were waves, the smooth flow term in the iWLP imposed finite u*, and consequently, finite transfer velocities (Figure 5).

Transfer Velocity Estimates from Level 4 Data
We performed simulations for the European costal ocean to compare between the ESM standard (the 'Wan92') and our extended physically-based alternatives.We show the comparison with the iWLP implemented with the kwind formulation by Jähne et al. [66] and the kbubble formulation by Woolf [44] using the kinematic viscosity of air (the W05va).Their kw estimates diverged mostly under unstable SBL, very rough or very smooth sea-surfaces, or higher friction velocities (Figure 6).The details of the simulations and the differences between kw estimates are presented hereafter.

Atmospheric Stability
Strong winds occurred along the European shores from the 24-26 May 2014.Besides, the air was unusually cold for the season and colder than the sea surface (Supplementary Video S1).Where the wind blew lighter, the rise of the warmer air heated by the sea surface generated turbulent eddies that enhanced mixing within the SBL and, consequently, also enhanced u* and kw.These unstable conditions, identified by Rib < 0, L tending to −0 and ψm < 0 (Supplementary Video S2), occurred more frequently and intensively nearby land masses and often associated to cooler continental breezes blowing off-shore.Its correct simulation required the estimation of the Rib, L and ψm from the algorithms by Grachev and Fairall [39] and Stull [67] that account for humidity considering saturation at 0 m heights.The Rib estimates neglecting humidity, following Lee [75], yielded biased estimates of the SBL conditions as a consequence of biased estimates of the virtual potential temperature.Stull [67] (page 9), Grachev and Fairall [39] and Fairall et al. [69] had already highlighted the importance of accounting for humidity.

Roughness Length
The sea surface agitation was very heterogeneous, particularly at the coastal ocean where it attained both the highest and the lowest estimated roughness lengths (the z0 in Supplementary Video S3).There, the steeper waves, a consequence of shorter fetches, should extract more momentum from the atmosphere under similar u10 conditions [54][55][56]58].Thus, the rougher coastal ocean surfaces were expected to have relatively more turbulent layers through which gases were transferred at higher rates.The comprehensive formulations simulated this by increasing u* (and consequently kwind) with z0 under similar uz i.e., similar winds generate more drag when blowing over rougher sea surfaces.Closer to the shores, and very often within the Mediterranean, lighter wind blew over smoother sea surfaces.Under these conditions, the iWLP estimated much higher z0 than the WLLP (Supplementary Video S4), demonstrating that the smooth flow was a fundamental driver for the z0 under calmer weather.The increase in z0 led to significantly higher u*, often 1.5 times higher and sometimes more,

Atmospheric Stability
Strong winds occurred along the European shores from the 24-26 May 2014.Besides, the air was unusually cold for the season and colder than the sea surface (Supplementary Video S1).Where the wind blew lighter, the rise of the warmer air heated by the sea surface generated turbulent eddies that enhanced mixing within the SBL and, consequently, also enhanced u * and k w .These unstable conditions, identified by Ri b < 0, L tending to −0 and ψ m < 0 (Supplementary Video S2), occurred more frequently and intensively nearby land masses and often associated to cooler continental breezes blowing off-shore.Its correct simulation required the estimation of the Ri b , L and ψ m from the algorithms by Grachev and Fairall [39] and Stull [67] that account for humidity considering saturation at 0 m heights.The Ri b estimates neglecting humidity, following Lee [75], yielded biased estimates of the SBL conditions as a consequence of biased estimates of the virtual potential temperature.Stull [67] (page 9), Grachev and Fairall [39] and Fairall et al. [69] had already highlighted the importance of accounting for humidity.

Roughness Length
The sea surface agitation was very heterogeneous, particularly at the coastal ocean where it attained both the highest and the lowest estimated roughness lengths (the z 0 in Supplementary Video S3).There, the steeper waves, a consequence of shorter fetches, should extract more momentum from the atmosphere under similar u 10 conditions [54][55][56]58].Thus, the rougher coastal ocean surfaces were expected to have relatively more turbulent layers through which gases were transferred at higher rates.The comprehensive formulations simulated this by increasing u * (and consequently k wind ) with z 0 under similar u z i.e., similar winds generate more drag when blowing over rougher sea surfaces.Closer to the shores, and very often within the Mediterranean, lighter wind blew over smoother sea surfaces.Under these conditions, the iWLP estimated much higher z 0 than the WLLP (Supplementary Video S4), demonstrating that the smooth flow was a fundamental driver for the z 0 under calmer weather.The increase in z 0 led to significantly higher u * , often 1.5 times higher and sometimes more, anticipating a significant impact on the local k wind estimates.These results demonstrated how the iWLP stabilized the k w in finite velocities as wind speed approached zero, and how this becomes a relevant matter when simulating the coastal ocean.

Transfer Velocities
In some rare situations the algorithms estimated unreasonably high k w despite the z 0 bounds imposed in the software.To avoid this bias, a 200 cm•h −1 ceiling was imposed on all k w estimates.The higher transfer velocities were associated to storms, both on the Atlantic side and on the Mediterranean side (Figure 7 and Supplementary Video S5).The root mean square deviation (RMSD) was estimated between the empirical u 10 -based Wan92 formula and the physically based W05va formula, averaged over time for each location.Results showed that both formulations diverged the most, in absolute terms, under storms and under strong winds near shore (Figure 7).On the other hand, the normalized root mean square deviation (NRMSD)-i.e., the deviation standardized to the Wan92 prediction-showed that, in relative terms, these formulations diverged the most under low and moderate winds and sea states near shore (Figure 7).The physically-based W05va formulation was also used to compare between the k wind and k bubble components of the CO 2 transfer velocity.The k bubble term was always lower than the k wind term, and at best 90% of k wind in the fetch-unlimited Atlantic or in the Mediterranean storms (Supplementary Video S5).Over the full domain, the k bubble term was on average 12% and on median 15.9% of the k wind term.Whatever the greenhouse gas, the differences were negligible between estimating k w using the single-layer or the double-layer schemes (Supplementary Video S5).Nevertheless, it is worth noting that the bigger differences were found under Mediterranean storm conditions.
Atmosphere 2018, 9, x FOR PEER REVIEW 13 of 18 anticipating a significant impact on the local kwind estimates.These results demonstrated how the iWLP stabilized the kw in finite velocities as wind speed approached zero, and how this becomes a relevant matter when simulating the coastal ocean.

Transfer Velocities
In some rare situations the algorithms estimated unreasonably high kw despite the z0 bounds imposed in the software.To avoid this bias, a 200 cm•h −1 ceiling was imposed on all kw estimates.The higher transfer velocities were associated to storms, both on the Atlantic side and on the Mediterranean side (Figure 7 and Supplementary Video S5).The root mean square deviation (RMSD) was estimated between the empirical u10-based Wan92 formula and the physically based W05va formula, averaged over time for each location.Results showed that both formulations diverged the most, in absolute terms, under storms and under strong winds near shore (Figure 7).On the other hand, the normalized root mean square deviation (NRMSD)-i.e., the deviation standardized to the Wan92 prediction-showed that, in relative terms, these formulations diverged the most under low and moderate winds and sea states near shore (Figure 7).The physically-based W05va formulation was also used to compare between the kwind and kbubble components of the CO2 transfer velocity.The kbubble term was always lower than the kwind term, and at best 90% of kwind in the fetch-unlimited Atlantic or in the Mediterranean storms (Supplementary Video S5).Over the full domain, the kbubble term was on average 12% and on median 15.9% of the kwind term.Whatever the greenhouse gas, the differences were negligible between estimating kw using the single-layer or the double-layer schemes (Supplementary Video S5).Nevertheless, it is worth noting that the bigger differences were found under Mediterranean storm conditions.Comparing transfer velocity (kw) algorithms using modelled data: (Wan92) transfer velocity of CO2 estimated from Wanninkhof [29] and averaged over the 66 h; (RMSD) Root mean square deviation between estimating the transfer velocity following Wanninkhof [29] or Woolf [44] conjugated with the iterative wind log-linear profile.(NRMSD) Root mean square deviation normalized to the Wan92 value at each time and location.

Conclusions
The sea state and atmospheric stability were determinant factors for the gas transfer velocities at the coastal ocean, irrespective of current algorithms still needing calibration and validation.Our findings corroborate recent conclusions that the wave state is fundamental for accurate estimates of gas transfer velocities at the fetch-limited coastal ocean [50,57].Furthermore, our conclusions also agree with Jackson et al. [88] and Shuiqing and Dongliang [49] when inferring about the open ocean.In fact, the COAREG (the gas-dedicated version of the COARE) by Jackson et al. [88] was applied to data from the GasEx cruises on the Northern Atlantic, Equatorial Pacific and Southern Ocean, with their results also suggesting that the sea state and atmospheric stability are important additional Comparing transfer velocity (k w ) algorithms using modelled data: (Wan92) transfer velocity of CO 2 estimated from Wanninkhof [29] and averaged over the 66 h; (RMSD) Root mean square deviation between estimating the transfer velocity following Wanninkhof [29] or Woolf [44] conjugated with the iterative wind log-linear profile.(NRMSD) Root mean square deviation normalized to the Wan92 value at each time and location.

Conclusions
The sea state and atmospheric stability were determinant factors for the gas transfer velocities at the coastal ocean, irrespective of current algorithms still needing calibration and validation.Our findings corroborate recent conclusions that the wave state is fundamental for accurate estimates of gas transfer velocities at the fetch-limited coastal ocean [50,57].Furthermore, our conclusions also agree with Jackson et al. [88] and Shuiqing and Dongliang [49] when inferring about the open ocean.In fact, the COAREG (the gas-dedicated version of the COARE) by Jackson et al. [88] was applied to data from the GasEx cruises on the Northern Atlantic, Equatorial Pacific and Southern Ocean, with their results also suggesting that the sea state and atmospheric stability are important additional factors driving k w in particular situations.Our application of the FuGas to the coastal ocean and the Jackson et al. ([88] Figure 11) application of the COAREG to the global ocean even agree in the magnitude of the enhancement of the transfer velocities, when compared to u 10 -based empirical parameterizations.Regarding the coarse resolution of Earth System Models, our FuGas 2.2 application demonstrated that physically-based algorithms diverge even more from empirical u 10 -based parameterizations when applied over the coastal ocean with finer spatial and temporal resolutions.Jackson et al. [88] also implemented a temporal fine resolution (3 h), having observed a wide discrepancy with the k w estimated from monthly means, and concluded that the k w estimated from monthly averaged factors is highly biased due to linearization of highly non-linear processes.The chemical enhancement of the transfer velocity is another factor that may become important under the low wind speeds often verified at the coastal ocean.Yet, it still needs to be implemented in the FuGas and tested regarding its functional form, the choice of gases, and the required environmental conditions [28,89,90].

Figure 1 .
Figure 1.Match-mismatch between estimated solubilities of greenhouse gases.Colour scale is quotient between kH estimated from Johnson [20] and from Sarmiento and Gruber [2].Solubilities tested for the range of values observable at the Earth's surface (a-c) and at the European coastal ocean from the 24 May 2014 at 06 h to the 27 May 2014 at 00 h (d-f) with each pixel corresponding to a square sided 11 km and averaged over the 66 h survey period.

Figure 2 .
Figure 2. Sea surface temperature (SST), salinity (S) and air pressure (P) at the European coastal ocean from the 24 May 2014 at 06 h to the 27 May 2014 at 00 h.Each pixel corresponds to a square sided 11 km and averaged over the 66 h survey period.

Figure 1 .
Figure 1.Match-mismatch between estimated solubilities of greenhouse gases.Colour scale is quotient between k H estimated from Johnson [20] and from Sarmiento and Gruber [2].Solubilities tested for the range of values observable at the Earth's surface (a-c) and at the European coastal ocean from the 24 May 2014 at 06 h to the 27 May 2014 at 00 h (d-f) with each pixel corresponding to a square sided 11 km and averaged over the 66 h survey period.

Figure 1 .
Figure 1.Match-mismatch between estimated solubilities of greenhouse gases.Colour scale is quotient between kH estimated from Johnson [20] and from Sarmiento and Gruber [2].Solubilities tested for the range of values observable at the Earth's surface (a-c) and at the European coastal ocean from the 24 May 2014 at 06 h to the 27 May 2014 at 00 h (d-f) with each pixel corresponding to a square sided 11 km and averaged over the 66 h survey period.

Figure 2 .
Figure 2. Sea surface temperature (SST), salinity (S) and air pressure (P) at the European coastal ocean from the 24 May 2014 at 06 h to the 27 May 2014 at 00 h.Each pixel corresponds to a square sided 11 km and averaged over the 66 h survey period.

Figure 2 .
Figure 2. Sea surface temperature (SST), salinity (S) and air pressure (P) at the European coastal ocean from the 24 May 2014 at 06 h to the 27 May 2014 at 00 h.Each pixel corresponds to a square sided 11 km and averaged over the 66 h survey period.

Figure 4 .
Figure 4. Transfer velocities (k w ) estimated for the Baltic experiment from eddy-covariance observations (E-C).The Jac99 (Jacobs et al.[81]) and WMG99(Wanninkhof and McGillis [82]) are the empirical u 10 -based parameterizations estimating higher and lower k w within this wind range.

Figure 6 .
Figure 6.Comparing the formulation by Wanninkhof [29] (Wan92) to the formulation by Woolf [44] (W05va) using modelled data about the European coastal ocean.Differences in the estimates of transfer velocity (kw) as functions of the bulk Richardson number (Rib), roughness length (z0) and friction velocity (u*).The 1:1 line (in black) represent equal kw estimates.

Figure 6 .
Figure 6.Comparing the formulation by Wanninkhof [29] (Wan92) to the formulation by Woolf [44] (W05va) using modelled data about the European coastal ocean.Differences in the estimates of transfer velocity (k w ) as functions of the bulk Richardson number (Ri b ), roughness length (z 0 ) and friction velocity (u * ).The 1:1 line (in black) represent equal k w estimates.

Figure 7 .
Figure 7.Comparing transfer velocity (kw) algorithms using modelled data: (Wan92) transfer velocity of CO2 estimated from Wanninkhof[29] and averaged over the 66 h; (RMSD) Root mean square deviation between estimating the transfer velocity following Wanninkhof[29] or Woolf[44] conjugated with the iterative wind log-linear profile.(NRMSD) Root mean square deviation normalized to the Wan92 value at each time and location.

Figure 7 .
Figure 7.Comparing transfer velocity (k w ) algorithms using modelled data: (Wan92) transfer velocity of CO 2 estimated from Wanninkhof[29] and averaged over the 66 h; (RMSD) Root mean square deviation between estimating the transfer velocity following Wanninkhof[29] or Woolf[44] conjugated with the iterative wind log-linear profile.(NRMSD) Root mean square deviation normalized to the Wan92 value at each time and location.