The FuGas 2 . 1 framework for atmosphere-ocean coupling in geoscientific models : improving estimates of the solubilities and fluxes of greenhouse gases and aerosols

Accurate estimates of the atmosphere-ocean balances and fluxes of greenhouse gases and aerosols are fundamental for 20 geoscientific models dealing with climate change. A significant part of these fluxes occur at the coastal ocean which, although much smaller than the open ocean, is also much more heterogenic. The scientific community is becoming increasingly aware of the necessity to model the Earth at finer spatial and temporal resolutions, which also requires better descriptions of the chemical, physical and biological processes involved. The standard formulations for the gas transfer velocities and solubilities are 24 and 36 years old, respectively, and recently, new alternatives have emerged. 25 We developed a framework congregating the geophysical processes involved which are customizable with alternative formulations with different degrees of complexity and/or different theoretical backgrounds. We propose this framework as basis for novel couplers of atmospheric and oceanographic model components. We tested it with fine resolution data from the European coastal ocean. Although the benchmark and alternative solubility formulations agreed well, their minor divergences yielded differences of many tons of greenhouse gases dissolved at the ocean surface. The transfer 30 velocities largely mismatched their estimates, a consequence of the benchmark formulation not considering factors that were proved determinant at the coastal ocean. Climate Change research requires more comprehensive simulations of atmosphere-ocean interactions but the formulations able to do it require further calibration and validation.

by the flow.The traditional formulation estimates the flux from F=k w •k H cp•Δp gas , in units of mol•m -2 •s -1 .The Δp gas is the difference between air and water gas partial pressures (atm).The k H cp is the Henry's constant for the gas solubility in its C w /p a form (mol•m -3 •atm -1 ), where p a is its air partial pressure (atm) and C w its concentration in the water (mol•m -3 ).
The k w is the transfer velocity of gases across the sub-millimetrically thick water surface layer in m•s -1 although usually 80 plotted in cm•h -1 .The alternative double layer model (Liss and Slater, 1974) estimates the flux taking into consideration both the water-side and air-side sub-millimetrically thick surface layers and thus, F = K w (C a /k H -C w ) = K a (C a -C w •k H ).
The C a and C w are the concentrations of the gas in air and water given in mol•m -3 and the k H is Henry's constant in its equivalent dimensionless quantity (C a /C w ).The transfer velocity is averaged over both layers from K w =(1/k w +1/(k H •k a )) - 1 or its equivalent K a =(k H /k w +1/k a ) -1 .85

Solubility
Sarmiento and Gruber (2013) compiled the algorithm for the k H cp dependence on temperature and salinity provided by Weiss (1974) and Weiss and Price (1980).We converted it to its corresponding dimensionless k H preserving the constants required to estimate Bunsen's solubility coefficient β.This formulation accounted for fugacity (f) of non-ideal 90 gases (Eq. 1) and corrected the gas partial pressure for moisture effects from the expression p moist =(1-p H2O /P)p dry considering water vapour saturation over the sea-surface (Eq.2).P is air pressure (atm), T w is water temperature (K), S is salinity (‰), p is the gas partial pressure (atm), R is the ideal gas law constant (Pa•m 3 • mol -1 •K -1 ), V m is the molar volume of the specific gas (22.3 for CO 2 and CH 4 , and 22.2432 for N 2 O) and V ideal =22.4136 mol•L -1 is the molar volume of ideal gases.Solubility coefficients were estimated from the Virial expansion (Eq.3), where Β was β or β/V m , 95 depending on which gas it was applied to (Table 3.2.2 in Sarmiento and Gruber (2013)).Our software automatically detected the gas from the a i coefficient.When Β=β the k H was estimated from Eq. (4).When Β=β/V m the k H was estimated from Eq. ( 5).(6) .33 3 -.
log H .

Transfer velocity
The available algorithms consider that the rate at which gases cross the sea-surface is basically set by the turbulence upon it.E.g. wind drag, wave breaking, currents and rain promote turbulence.The water viscosity, set by temperature and salinity and enhanced by the presence of surfactants, antagonizes turbulence.With all these forcings, it becomes difficult to develop an algorithm that estimates the transfer velocity accurately.The literature has many of them, either 125 fitted to specific surface conditions or rougher generalizations, focusing on different factors and relying in different theoretical backgrounds.The simpler ones rely on the wind velocity 10m above the sea-surface (u 10 ).Among then, the formulation by Wanninkhof (1992)  Other simple empirical formulations based only on u 10 (Carini et al., 1996;Raymond and Cole, 2001), or also accounting for current drag with the bottom (Borges et al., 2004), 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. (2013), with the k w being decomposed into its shear produced turbulence (k wind ) and bubbles from whitecapping (k bublle ) forcings (Asher and 140 Farley, 1995;Borges et al, 2004;Woolf, 2005;Zhang et al., 2006).The effect of currents was disregarded at this stage (Eq.10).Sc w was determined from temperature and salinity following Johnson (2010).The formulation by Zhao et al. (2003), merged k wind into k bubble (Eq.11a) using the wave breaking parameter (R B given by Eq. 11b).The u * is the friction velocity i.e, the velocity of wind dragging on the sea-surface, and f p is the peak angular frequency of the wind-waves.The kinematic viscosity of air (υ a ) was estimated from Johnson (2010).This solution used the wave field as a proxy for whitecapping that increased transfer velocity with wind-wave age.However, it simultaneously used the wave field as a proxy for the sea-surface roughness that increased transfer velocity from wind-drag over steeper younger waves (through the WLLP estimation of u * explained in a section below).bubble .
A more comprehensive solution split the two drives of transfer velocity (Woolf, 2005;Zhang et al., 2006): k wind for the transfer mediated by the turbulence generated by wind drag (Eq.12) (Jähne et al., 1987) and k bubble for the transfer mediated by the bubbles generated by breaking waves (Eq.13) (Zhang et al., 2006).Β is Bunsen's solubility coefficient estimated for the local sea-surface conditions.W=3.88×10 -7 R B 1.09 is the whitecap cover requiring the R B estimated from (Eq. 11b), V=4900, e=14 and n=1.2.
These formulations required friction velocity (u * ), which was estimated from the Wind Log-Linear Profile (WLLP: Eq. 14) accounting for wind speed at height z (u z ), atmospheric stability of the surface boundary layer (through ψ m ) and seasurface roughness (through the roughness length z 0 ).The κ is von Kármàn's constant.
ln ln (14) Roughness length (z 0 ) is the theoretical minimal height (most often sub-millimetrical) 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 (2001) proposed a dimensionless z 0 dependency from the wave field, increasing with the wave slope (Eq.15).Due to the bidirectional nature of the z 0 and u * relation, we also tested an iterative solution (iWLP) where Eq.15 was used as a first guess for the z 0 and Eq.14 for its subsequent u * .A second iteration re-estimated z 0 from the COARE 3.0 (Fairall et al.;2003) adaptation of the Taylor and Yelland (2001) formulation, which added a term for smooth flow (Eq.16), and u * again from Eq.14.Applying four iterations were enough for an excellent convergence of the full data array.in its turn estimated from air temperature, air pressure and specific humidity (Grachev and Fairall, 1997) or from the liquid water mixing ratio (Stull, 1988).Alternatively, the use of the air potential temperature neglected humidity (Lee, 1997).The wind velocity (u z ), temperature (T z ), pressure (P z ) and humidity (q z ) z meters above sea-surface were given by the WRF second level.The wind velocity at z 0 (u 0 ) was set to the theoretical u 0 =0.Temperature at the height of 0 m (T 0 ) was given by the SST (Grachev and Fairall, 1997;Fairal et al., 2003;Brunke et al., 2008) without rectification for 190 cool-skin and warm-layer effects due to the lack of some required variables.Yet, these effects tend to compensate each other (Brunke et al., 2008;Fairall et al., 1996;Zeng and Beljars, 2005).Air pressure at 0 m (P 0 ) was given by the WRF at the lower first level (at roughly 0 m).Humidity at 0 m (q 0 ) was set to the saturation level at P 0 and T 0 (Grachev and Fairall, 1997).The Ri b was used to estimate the length L from, Monin-Obukhov's similarity theory, a discontinuous e ponential function tending to ± ∞ when Ri b tends to ±0 and tending to ±0 when Ri b tends to ±∞.Ri b and L were used 195 to estimate ψ m following Stull (1988) or Lee (1997) algorithms.

Ri (17)
CO 2 is mildly soluble with a K H =1.17 for pure water at 25 ºC.Its transfer velocity is limited by the molecular crossing 200 of the water-side surface layer.CH 4 is much less soluble with a K H =31.5 for pure water at 25 ºC.Its transfer velocity should also take into consideration the molecular crossing of the air-side surface layer (Johnson, 2010).We compared between the use of the traditional single layer and the double layer "thin film" model (Liss and Slater, 1974;Johnson, 2010;Vieira et al, 2013), the later requiring the air-side transfer velocity (k a ) estimated from the COARE formulation as in Eq. 18 (Jeffrey et al., 2010).CD is the drag coefficient and Sc a the Schmidt number of air, which were determined for 205 a given temperature and salinity following Johnson (2010).

Validation with field data 210
The field sampling occurred from the 22 nd of May 2014 to the 26 th of May 2014 using the atmospheric tower at Östergarnsholm in the Baltic Sea (57° ′ N, 18° ′ E), the Submersible Autonomous Moored Instrument (SAMI-CO 2 ) 1 km away and the Directional Waverider (DWR) 3.5 km away, both south-eastward from the tower (see e.g.Högström et al. (2008) and Rutgersson et al. (2008) for detailed description of the sites).The air-water CO 2 fluxes measured by eddy-covariance were smoothed over 30 min bins and corrected according to the Webb-Pearman-Leuning 215 (WPL) method (Webb et al., 1980).We used only the fluxes for which the wind direction set the SAMI-CO 2 and DWR in the footprint of the atmospheric tower (90º < wind direction < 180º).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.We applied this data set to the single processing software ensemble of the FuGas 2.1 in order to test which algorithms provide better approximations to 220 reality.

2.4
Atmosphere-ocean coupler The atmospheric model was the standard operational application of the WRF by Meteodata.cz, with 9 km and 1 h resolutions.Air temperature 'T' (º ), pressure 'P' (atm), U and V components of wind velocity (m•s -1 ), water vapour mi ing ratio 'Q' (scalar) and height 'h' (m), where retrieved at the two lowest levels within the atmospheric surface boundary layer (SBL).The vertical thickness of the WRF horizontal layers varied with space and time.Over the ocean, the two lowest levels occurred roughly at 0 m and 12 m heights.The WRF output decomposes height, temperature and pressure into their base level plus perturbation values.
Sea-surface temperature (SST) 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 INGV using the WW3-NEMO modelling system at 0.0625 º and 1 h resolutions (Clementi, 2013), and for the North Atlantic by Windguru at roughly .º and 3 h resolutions.The variables included significant wave height 'H s ' (m) and peak frequency 'f p ' (rad•s -1 ) 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 ' (m) was estimated from the peak frequency assuming the deep-water approximation: L p = πg/f p 2 , where g is the gravitational acceleration constant; (ii) the Windguru data did not provide wind sea component (where and) when the wind was too low.For these missing cases were attributed the lowest H s and L p simulated everywhere else; (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.
The WRF and WW3-Nemo outputs were retrieved for the European shores from the 24 th of May 2014 at 06h to the 27 th of May 2014 at 00 h.All variables were interpolated to the same .º grid (roughly km at Europe's latitudes) and 1 h time steps.This resulted in a data set with 17 variables × 41776 locations × 66 time instances, that occupied nearly 1Gb ram memory (with another 1Gb taken by the software).To optimize the computations, the calculus was first vectorized and then parallelized using the Single Program Multiple Data (spmd) programming strategy.Hence, in the FuGas 2.1 multiple processing software ensemble, the variables were first organized in matrices with locations along the 1 st dimension and time along the 2 nd .Running the calculus applying matrix algebra to the whole data set, by itself represented an improved speed of several orders of magnitude.Furthermore, the spmd replicated the data, split the replicates into n approximately equal-sized arrays, and distributed their calculus among the n available cpu cores, which represented an extra improvement of computational speed.However, it also bared computational costs: (i) invoking the parallel processing toolbox was time consuming, (ii) replicating 1Gb ram was time consuming, (iii) once running the calculus, the 4Gb ram memory was soon exhausted, which stall the calculus, (iv) to avoid it, the spmd were split into several sequential code blocks and in-between the variables no longer necessary were deleted.This spmd fragmentation was time consuming.In conclusion, there is no perfect solution for calculus parallelization, and although spmd is the best strategy available for this task, its application needs to be carefully programmed according to the data and hardware characteristics.

Results
Both solubility formulations were tested simulating T w from 4ºC to 30ºC at 1ºC intervals and S from 0 ppt to 36 ppt at 1 ppt intervals.The metric k H , Joh10 /k H , Sar13 showed better how much their estimates could diverge (Fig. 1).Afterwards, both formulations were applied to the data from the European coastal ocean.Their estimates were compared applying the previous metric averaged over the 66h time interval using the geometric mean (Fig. 1).From the 24 th to the 26 th May the water temperature at the ocean surface changed significantly and there were large fresh water inputs from the Black Sea and the Baltic Sea (Video 1).The widest divergences were up to 4.5% in the CO 2 solubility estimates associated to cooler waters, 5.8% in the CH 4 solubility estimates associated to both temperature extremes, and 2.1% in the N 2 O solubility estimates associated to cooler and less saline waters (Fig. 1).These mismatches lead to large differences in the estimates of greenhouse gases dissolved in the first meter below the ocean surface (Fig. 2).These The k w estimated from the E-C measurements presented a systematic bias.To detect its source, the difference (∆k w ) between the k w estimated from the E-C measurements and the one estimated from the Wan92 formulation was compared to the potential sources of bias.Besides well correlated with u 10 (r=0.55), the ∆k w was also well correlated with the relative humidity (r=-0.7)and with the first (r=0.49),second (r=0.47) and third (r=0.67)terms of the WPL correction.The distortion of the E-C flux estimates by cross-sensitivity to humidity is a common problem with openpath IRGA, raising substantially their detection limit.The observed differences between the concentrations of CO 2 in the air and in the water during our survey varied within 120 and 270∆ppm, well below the limit for a 25% error in the flux estimates as reported by Blomquist et al. (2014) for our IRGA model, the LI-COR LI-7500.We hypothesize whether the E-C data lacked quality to calibrate and validate the formulations.However, our formulations were close matches to the estimates by widely used transfer velocity formulations subject to thorough calibration and validation, which proved them reasonable estimators of the central tendency (Fig. 3).Hence, we were confident about the potential of our newly proposed formulations to replicate the central tendency similarly well while improving the accuracy of the estimates for each particular location.
During this Baltic Sea sampling at the Östergarnsholm site, the SBL was generally stable (0<Ri b <0.5 with a few exceptions) and the sea-surface was little to moderately rough (z 0 <0.49mm).These conditions were used as reference to estimate the elasticity of k w to its forcing functions (Fig. 4).The variables related with the SBL stability, namely the u 10 , temperature, pressure and humidity, were the variables able to induce larger changes in k w .Several renowned u 10 -based formulations for the estimation of k w were used and compared with the most comprehensive alternatives provided in our software and framework (Fig. 3).Although their estimates were close matches, there were a few fundamental differences: the comprehensive algorithms split the data points into two distinct scatter lines, the upper line for k w obtained under rougher sea-surfaces and the lower line for k w obtained under smother ones.The red markers representing the ZRb03 iWLP give the best example.The u 10 -based formulations were unable to perform this adjustment to the local wave state.Their small k w fluctuations were a sole consequence of changes in water viscosity (as estimated by the Sc w ) driven by changes in water temperature.These results highlight the potential of the SBL stability and the sea-surface agitation as additional k w mediators.It is curious that the wave variables were the responsible for the big differences between k w estimates (as shown in Fig. 3) although these were the variables to which the k w was least elastic (as shown in Fig. 4).It demonstrates that more important than model sensitivity (or elasticity) is how much the There is yet the interesting detail of how the WLLP and the iWLP diverged under smoother sea-surfaces, supporting the solution suggested in the COARE 3.0 (Fairall et al., 2003) for the iterative estimation of u * and z 0 .
Complementary to the analysis above, we also used the simulations of the European costal oceans to compare between the ESM standard (the Wan92) and one of our comprehensive alternatives (the iWLP-ZRb03), chosen on the basis of two factors: it was both the most elastic formulation and the one providing the closest estimates to the Wan92 (recall Fig. 3).Since the Wan92 often represented the central tendency of the iWLP-ZRb03, this choice provided the best probability that the differences between the k w estimates were due to the enhanced representation of the environmental processes involved and not to systematic biases associated to uncertainty in the parameter estimation.
Both k w estimates diverged under two particular situations (Fig. 5): (i) under low winds and unstable SBL, and (ii) under high winds and rougher sea-surfaces.
Strong winds occurred along the European shores from the 24 th to the 26 th of May of 2014.Besides, the air was unusually cold for the season and colder than the sea-surface (Video 1).The upward advection of the warmer air, heated by the sea-surface, generated turbulent eddies that enhanced mixing within the SBL.These unstable conditions were identified by Ri b <0, L tending to -and ψ m <0 (Video 2).The mixing of the SBL enhanced u * and k w everywhere the wind blew lighter.This situation 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 (1997) and Stull (1988) that account for humidity considering saturation at 0 m heights.The Ri b estimates neglecting humidity (Lee, 1997) often yielded neutral conditions (i.e, with Ri b ≈ ) or unreasonably stable SBL (i.e, with Ri b >0).
The sea-surface agitation was very heterogenic, particularly at the coastal ocean where it attained both the highest and the lowest estimated roughness lengths (the z 0 in Video 3).There, the steeper waves as a consequence of shorter fetches, should extract more momentum from the atmosphere under similar u 10 conditions (Taylor and Yelland, 2001;Fairall et al, 2003).Thus, the rougher coastal ocean surfaces were expected to possess 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 harsher seasurfaces.Aside the rougher weather, whenever lighter wind blew over smoother sea-surfaces, the iWLP estimated much higher z 0 than the WLLP (video 4), demonstrating that the smooth flow was a fundamental driver for the z 0 under calmer weather.This increase in z 0 lead to significantly higher u * , often 1.5 times higher and sometimes more, anticipating a significant impact on the k wind estimates.
The comprehensive formulation (i.e, ZRb03 iWLP) often estimated k w largely higher than the one estimated by the ESM standard formulation (i.e, Wan92), although it occasionally estimated lower k w (Video 5).Its largest estimates of k w were associated to unreasonably high estimates of z 0 that biased the subsequent results.These biased estimates of z 0 could either be due to a poor calibration of the Taylor and Yelland (2001) model estimating z 0 from the wave field or due to biased wave field provided by the WW3-NEMO.To avoid this bias, k w was imposed a cm•h -1 ceiling, corresponding to the maximum reported in the bulk literature associated to similar wind speeds.With this restriction, the difference in the CO 2 volume transferred by either formulation across the ≈ , , km 2 of ocean surface during the 66 h was of 12997 km 3 , corresponding to 33.7% of the 38551 km 3 of CO 2 total volume transferred using the ESM standard formulation (Fig. 6).These differences were higher at the coastal ocean, a consequence of the factors that were and 41158 Km 3 , respectively.The differences were negligible between using the single layer or the double layer scheme to estimate k w , even for a rather insoluble gas as is CH 4 (Fig. 6 and Video 5).Nevertheless, it is worth noting that it was again in the fetch limited coastal ocean where most of the bigger differences were found.

4 Discussion
The accurate estimation of the balances of greenhouse gases and aerosols in the atmosphere and in the oceans, as well as their fluxes across the surfaces of the coastal oceans, is an important issue for biogeosciences and Earth-System modelling (ESM).Previous estimates of CO 2 uptake by the global oceans done by coarse resolution implementations diverged in about 70 % depending on the transfer velocity formulations being used (Takahashi et al., 2002), whereas the 350 wide uncertainty in the ocean N 2 O source to the atmosphere mostly originated from the uncertainty in the air-water transfer velocities (Nevison et al., 1995).However, the knowledge on this subject is still limited, with plenty of room for improvement.As an example, the simpler formulations for the estimation of k w assume either a quadratic or cubic dependency from u 10 depending mostly on the sensing method, time scale and fetch at the particular location.
Furthermore, the simulation of atmosphere-ocean interactions by regional and Earth-system models, by still using these 355 simpler formulations, are decades behind the state-of-the-art.Our work proposes a framework to incorporate this stateof-the-art in an atmosphere-ocean coupler and demonstrates that this is fundamental for reliable simulations of coastal ocean systems.
Remarkably, both solubility formulations generally matched their estimates despite their distinct backgrounds.
Nevertheless, they did diverge in as much as 0.045 mol•mol -1 of CO 2 , .mol•mol -1 of CH 4 and .mol•mol -1 of 360 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: (i) the cooler marine waters occur closer to the poles, where the solubility pump traps greenhouse gases and carries them to the deep ocean (Sarmiento and Gruber, 2013), and (ii) the warmer and the less saline waters occurring at the coastal ocean and seas, which have regularly been observed having greenhouse gases and aerosols dissolved in concentrations highly unbalanced with those of the atmosphere (Nevison et 365 al., 2004;Borges et al, 2005;Barnes and Upstill-Goddard, 2011;Sarmiento and Gruber, 2013;Dutta et al., 2015;Gypens and Borges, 2015;Harley et al., 2015).Therefore, the biases in the estimated total amount of greenhouse gases in the first meter depth of the European coastal ocean during late May 2014 may be an indicator of higher global biases.This work showed that the accurate estimation of the transfer velocity of greenhouse gases and aerosols across the coastal oceans' surface requires taking into consideration at least the atmospheric stability of the SBL and the sea-370 surface roughness, as recently suggested by Jackson et al. (2012) and Shuiqing and Dongliang (2016).Our results show that, by neglecting these factors, the simpler u 10 -based formulations tend to provide lower estimates of the transfer velocity than the provided by comprehensive formulations.Similar conclusions were achieved by Jackson et al. (2012).
However, the more comprehensive formulations still need improvement and validation.It is imperative to calibrate and validate the estimation of transfer velocity from friction velocity and wind-wave breaking, and the roughness length 375 from the wave field.All the available formulations for these specific purposes lack robust parameter estimations.
Generally, there seems to be a great dependency of the available algorithms from the particular data sets that were used to calibrate them.Nevertheless, there is a general consensus that the k bubble term is fundamental under high wind speeds, with its estimate being central to current k w research.The latest developments have been on the dependency of k bubble Geosci. Model Dev. Discuss., doi:10.5194/gmd-2016-273, 2016 Manuscript under review for journal Geosci.Model Dev.Published: 23 November 2016 c Author(s) 2016.CC-BY 3.0 License.from the interactions among the wind, the wave state, the bubble plume and the properties of the gas being transferred (Woolf et al. 2007;Callaghan et al., 2008Callaghan et al., , 2014;;Goddijn-Murphy et al., 2011, Crosswell, 2015).The effect of sea-spray is the new buzz on this topic and only recently started emerging algorithms like the ones by Zhao et al. (2006) andWu et al. (2015).So far, these focused on the momentum transfer from wind to the ocean surface and the attenuation of the friction velocity.It should be interesting to understand how the intrusion of the sea-spray on the atmosphere affects the transfer velocity of gases, being anticipated a process symmetrical to that of the intrusion of bubbles on the ocean.The new algorithms for the effects of surfactants are particularly concerned with the variability of the coastal ocean (Pereira et al., 2016).These no longer associate the surfactants to the Schmidt number's e ponent but rather to a coefficient setting a proportional decay of k w .The effect of sea-ice must take into consideration its distortion of the ocean surface and its effect upon the SBL stability (Loose et al., 2014).Our coupling solution still needs to integrate the effects of the sea-surface cool-skin and warm-layer, surfactants, rain, sea-spray and sea-ice.From these, the cool-skin and warm-layer algorithms are the only with robust calibrations and validations, mostly done under the COARE (Fairall et al., 1996;Fairal et al., 2003;Zeng and Beljars, 2005;Brunke et al., 2008).The addition of complexity to any coupling solution must be carefully thought as these cannot become intricate to the point of calculus becoming unbearable for ESM application.In particular, any algorithm demanding for-loops is unviable as it disables calculus vectorization and its coordination with parallel processing.In our software, vectorization enabled improving calculus roughly 12× faster in a single core.
Besides finding the appropriate algorithms and parameter values to be used by the coupler, there is also the issue of accurately retrieving the variables mediating the gas transfer.The results showed that the k w was most elastic to the variables related with the SBL stability, namely the u 10 , temperature, pressure and humidity.Although these are provided by the oceanic and atmospheric model components at courser vertical resolutions, they need to be transposed to finer vertical resolutions taking into consideration the processes occurring at the sea-surface.While the u 10 is given by the atmospheric model, the water temperature needs to account for the cool-skin and warm-layer effects and the heat and humidity at the SBL need to account for their vertical fluxes over the sea-surface.The COARE algorithm is the state-of-the-art for these tasks.During most of its development it focused on E-C methods to estimate the fluxes of heat and humidity across the SBL using a framework with an intricate mathematical structure going deeper into the simulation of the geophysical process.Given its complexity, it must be quite a challenge to perform the calculus vectorization and parallelization required for the substantial improvement of computational speed and its application to ESM.Only in its latter developments did the COARE explicitly addressed the fluxes of gases and the importance of sea-surface roughness (Fairall et al, 2003;Jeffrey et al., 2010;Blomquist et al., 2006Blomquist et al., , 2014)).
developed an algorithm from an alternative chemistry background.It accounts for the effects of temperature and salinity taking into consideration the molecular and thermodynamic properties of the water, its solutes and the specified gas, but disregarding the non-ideal behaviour of the gases and moisture.His formulation was developed from the compilation by Sander (2015) (although available in 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 110 temperature and 0 ppt salinity.The term -Δ soln H/R reflected the temperature (in Kelvin) 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 (Eq.7) relied on the Geosci.ModelDev.Discuss., doi:10.5194/gmd-2016-273,2016   Manuscript under review for journal Geosci.Model Dev.Published: 23 November 2016 c Author(s) 2016.CC-BY 3.0 License.empirical Setschenow constants (K S =θ•logVb) reporting the effect of electrolytes salting-out gases proportionally to their liquid molar volume at boiling point (Vb).The Vb was estimated using the additive Schroeder method, whereas θ was estimated from Eq.8 using a provisional k H# =0.0409/k H cp.
(henceforth also mentioned as 'Wan ') became the standard used in ESM and satellite data processing.It further considers the Schmidt number of the water (Sc w ) related to viscosity and with its e ponent reflecting the surface layer's rate of turbulent renewal, and the temperature dependent chemical enhancement 130 due to CO 2 reaction with water (α Ch ): Geosci.ModelDev.Discuss., doi:10.5194/gmd-2016-273,2016   Manuscript under review for journal Geosci.Model Dev.Published: 23 November 2016 c Author(s) 2016.CC-BY 3.0 License.
differences were estimated from Δton•m -1 •121 km -2 = 11 2 •Δs•p gas •P• 3 •M a /(10 9 •R•T), where Δs was the difference in the solubility estimated by either algorithm in its C w /C a form at each 11 km wide cells and averaged over the 66 h time interval.M a =28.97 was the air molecular mass and p gas the atmospheric partial pressure of CO 2 , CH 4 or N 2 O, 390 ppm, 1.75 ppm and 0.325 ppm respectively(EPA, 2015), assuming that they were approximately uniform all over the atmospheric SBL.These differences summed to 3.86×10 6 ton of CO 2 , 880.7 ton of CH 4 and 401 ton of N 2 O.Because the bias of N 2 O changed from positive to negative with location, the overall bias was 163 ton.
Geosci.Model Dev.Discuss., doi:10.5194/gmd-2016-273,2016 Manuscript under review for journal Geosci.Model Dev.Published: 23 November 2016 c Author(s) 2016.CC-BY 3.0 License.not taken into consideration by the ESM standard.The total volumes of CH 4 and of N 2 O transferred were 41156 Km 3

FigureFigure 5 :
Figure 4: Elasticities of the transfer velocity 575 to its forcing functions.

Figure 5 :
Figure 5: Applying the modelled data about the European coastal ocean for a direct comparison between the kw 620 estimates by the ESM standard and our best performing comprehensive.